© ESO, 2012

1. Introduction

Massive stars play an important role in cycling and mixing material between stars and the interstellar medium in galaxies. However, our understanding of how high-mass stars (M ≥ 8   M_⊙) form is still far from complete (Zinnecker & Yorke 2007). Observationally, it is challenging to gain information from high-mass young stellar objects (YSOs) because high-mass stars usually form in groups at large distances (typically a few kpc), requiring high angular resolution to resolve individual sources. High-mass YSOs are embedded in regions with high extinctions (~10–100 A_v) for  ~15% of their lifetime, where complex processes, such as gravitational interaction, powerful outflows and stellar winds, high accretion rates and ionizing radiation fields make the interpretation of observations difficult. Moreover, high-mass YSOs evolve very fast (~10⁵ yrs) and reach the main sequence while still in the embedded phase. Structures such as circumstellar disks have an even shorter lifetime and are harder to detect.

Accretion via a circumstellar disk is suggested to be the most effective way to overcome radiative pressure and build up a high-mass star (e.g., Krumholz et al. 2009). The detection of massive bipolar outflows with high mass loss rates of  ~10^-4 M_⊙ yr^-1 around young high-mass stars implies the presence of circumstellar disks with comparable mass accretion rates (e.g., Henning et al. 2000; Beuther et al. 2002; Zhang et al. 2005). To date, only a dozen or so B-type young stars have been proposed to have circumstellar disks (e.g., Cesaroni et al. 1997, 1999, 2006; Shepherd & Kurtz 1999; Beuther et al. 2005), with disk masses less than or comparable to the masses of the host stars. In contrast, around young O-type stars, only large rotating toroids are found with masses greater than stellar masses (Beltrán et al. 2005; Furuya et al. 2008).

The presence of circumstellar disks around low-mass young stars is firmly established, and our current observational and theoretical understanding is reviewed by Williams & Cieza (2011) and Armitage (2011), respectively. In contrast, the applicability of a scaled-up version of low-mass star formation via disk accretion for massive stars has not been firmly established (Zinnecker & Yorke 2007), although promising cases have been reported (e.g., Keto & Zhang 2010). Clearly, more observational work is needed to establish the kinematics of material associated with high-mass YSOs.

In this paper, we present observations of massive disk candidate, AFGL 2591–VLA3, a young early B-type star in the Cygnus X region. The distance to this object is uncertain, with reported values ranging from 0.5 to 2 kpc (0.5–2 kpc as discussed in van der Tak et al. 1999) and recent measurements of water masers even suggesting 3.3 kpc (Rygl et al. 2012). In this work, we assume a distance of 1 kpc in order to aid comparison with previous publications and discuss how our conclusions change for a distance of 0.5 or 3 kpc. At 1 kpc, the total luminosity of AFGL 2591 is  ~2 × 10⁴ L_⊙ (Lada et al. 1984). Compact radio emission is associated with the source (named as VLA1, VLA2 and VLA3; Campbell 1984; Trinidad et al. 2003). Modeling of its free-free emission implies the stellar mass is about 16 M_⊙ (van der Tak & Menten 2005). An envelope mass of only 40 M_⊙ within 30 000 AU around the source is derived through comprehensive radiative transfer modeling of molecular lines by van der Tak et al. (1999).

The geometry of the system can be derived from bipolar ¹²CO outflows at arcminute scale, which extend in east (red)-west (blue) direction (position angle (PA)  ~  270°–325°) observed at resolutions of 14″ to 21″ by Mitchell et al. (1992) and Hasegawa & Mitchell (1995). Observations of ¹³CO at resolution of 14″ suggests that the PA of the outflow is about 280° (van der Tak et al. 1999), confirmed by K-band bispectrum speckle imaging of AFGL 2591 at a resolution of  ~0.2″ showing multiple loops on the western side of the source with PA  ~ 260° (Preibisch et al. 2003). Precession of the outflow axis is also implied. Dense gas in the outflow walls traced by CS and HCN shows a consistent morphology (Bruderer et al. 2009). Based on this orientation of the outflow, a north-south (PA  ~ 0°) orientation of the equatorial plane of AFGL 2591–VLA3 is inferred.

Observations of HDO 1_1,0–1_1,1, H₂¹⁸O 3_1,3–2_2,0 and SO₂ 12_0,12–11_1,11 with the Plateau de Bure interferometer (PdBI) of the Institut de Radio Astronomie Millimétrique (IRAM)¹ reveal a velocity gradient in the north(east)-south(west) direction across a barely resolved source with a diameter of  ~800 AU, suggestive of a circumstellar disk (van der Tak et al. 2006). However, the kinematics could not be studied due to the limited angular resolution (~1″). Here we present new PdBI observations toward AFGL 2591 at a higher resolution of  ~0.5″ (~500 AU at 1 kpc), which resolve the kinematics. This paper is organized as follows. Section 2 summarizes the PdBI observations and data reduction. Section 3 presents the results. In Sect. 4, we describe a simple model of the kinematics. Discussions and conclusions are given in Sects. 5 and 6, respectively.

2. Observations

The source AFGL 2591 was observed with PdBI in 2006–2007 in 4 tracks. The 3 mm receiver was tuned to cover HDO 1_1,0–1_1,1 at 80 578.3 MHz (hereafter HDO) and the 1.3 mm receiver was tuned to cover H₂¹⁸O 3_1,3–2_2,0 at 203 407.5 MHz (hereafter H₂¹⁸O) and SO₂ 12_0,12–11_1,11 at 203391.6 MHz (hereafter SO₂). Table 1 presents a log of the observations and Table 2 summarizes the observed lines. The combined projected baselines range from 26 kλ (~10″) to 206 kλ in length at 3 mm and from 48 kλ (~5″) to 525 kλ in length at 1.3 mm. The phase tracking center was chosen at α(2000) = 20^h29^m2487, δ(2000) = 40°11′1950. The nominal LSR velocity adopted for observations was  − 5.5 km s^-1. The quasars 2005+403 and 2013+370, which are about 5° away from AFGL 2591, were frequently observed for phase and gain calibrations in all four observations. 3C 273, 3C 345 and 3C 454.3 were used as the bandpass calibrators, and MWC349 with a total flux density of  ~ 1.3 Jy at 203.4 GHz during the time of the observations, as flux calibrator. The flux error is estimated to be within 20%. The spectral resolutions are 39.0625 kHz per channel or 0.145 km s^-1 per channel at 3 mm band, and 78.125 kHz per channel or 0.115 km s^-1 per channel at 1.3 mm band. Data reduction was conducted at the IRAM headquarters in Grenoble, using the GILDAS² software.

Fig. 1 Natural-weighted continuum images in the region of AFGL 2591. (Left) 80.6-GHz continuum. Contours are at 3, 5, 10, 20, 30... σ with 1-σ of 0.15 mJy beam-1. Sources within the field following Campbell (1984) and Trinidad et al. (2003) are marked on the map. (Right) 203.4-GHz continuum. Contours are at 3, 5, 10, 30, 50... σ with 1-σ of 0.4 mJy beam-1.

Fig. 2 Visibility analysis of the continuum at 203.4 GHz. a), c) Show the observed visibilities (filled circles) and model visibilities (grey area). The observed data are vector-averaged for the purpose of better visualization. Only the averaged structure in circular symmetry can be observed in the plot. Full model visibilities are plotted in order to show the non-circular symmetric structure of the source. b), d) Represent a graphical view of the visibility fits of VLA3 shown in a) and c), respectively. The red plus sign is the peak position of the emission. Blue ellipse and filled circle stand for the Gaussian component and the point source.

3. Results

3.1. Continuum emission

The natural-weighted continuum images of AFGL 2591–VLA3 at 80.6 GHz and 203.4 GHz are presented in Fig. 1. The angular resolutions are , PA 28° at 80.6 GHz and , PA 18° at 203.4 GHz. With an arcsecond beam, VLA3 is unresolved at 80.6 GHz. In contrast, at 203.4 GHz VLA3 is well resolved with a subarcsecond beam and shows a diamond shaped morphology in the north-south direction. In addition to VLA3, VLA1 and VLA2 are detected within the field of view at 80.6 GHz. At 203.4 GHz, VLA1 is likely resolved out due to missing short spacings. A previously unknown source (VLA3–NE) located at  ~2.5″ north-eastern to VLA3 is detected at 5-σ at 203.4 GHz.

3.1.1. Measurements of positions, sizes and total flux densities

The positions, sizes and total flux densities of the sources in the region of AFGL 2591 are derived in the visibility domain and summarized in Table 3. Judging from the images in Fig. 1, all sources are fit with point sources, except VLA1 at 80.6 GHz and VLA3 at 203.4 GHz, where Gaussians are used. The vector-averaged visibility plot in Fig. 2 shows a compact component on baselines longer than 200 kλ for VLA3 at 203.4 GHz. We therefore fit VLA3 with a 20 mJy point source plus a 65 mJy, 1.0″  ×  0.6″ (~800 AU diameter at 1 kpc) Gaussian source. Including the flux densities of VLA2 and VLA3–NE, the total flux density within the field of view at 203.4 GHz is about 93 mJy which is consistent with the visibility measurement (Fig. 2).

Compared to van der Tak et al. (2006), our 80.6-GHz observations measure 100%, 12% and 45% of the total flux densities of VLA1, VLA2 and VLA3, respectively, while the new 203.4-GHz observations detect 0%, 7% and 44% of the total flux densities of VLA1, VLA2 and VLA3, respectively. Considering the different array configurations (limited by the antenna shadowing effect, the shortest baselines of the observations made by van der Tak et al. 2006, are 8 kλ at 3 mm and 20 kλ at 1.3 mm) and the absolute uncertainty of the flux density, the new data are consistent with optically thin free-free emission, optically thick dust or free-free emission and optically thin dust emission for VLA1, VLA2 and VLA3, respectively, as discussed by Trinidad et al. (2003) and van der Tak et al. (1999, 2006). For the newly detected source VLA3–NE, more sensitive high-resolution observations at different frequencies are needed to derive its nature.

3.1.2. Estimation of densities and masses

For the “dust” source VLA3, we derive the beam averaged H₂ column density from the measured peak intensity of 30 mJy per beam at 203.4 GHz. The H₂ column density at the emission peak can be estimated as (1)where S_ν is the peak flux density, a is the gas-to-dust ratio (100), Ω_b is the beam solid angle, m_H is the mass of atomic hydrogen, κ_ν is the dust opacity per unit mass and B_ν(T_d) is the Planck function at dust temperature T_d. We apply the value of κ_ν(1.3   mm) ≈ 1.0 cm² g^-1 as suggested by Ossenkopf & Henning (1994) for gas densities of 10⁶ − 10⁸ cm^-3 and coagulated dust particles with thin ice mantles. The beam averaged H₂ column densities for T_d from 100 to 300 K (see Sect. 3.2.4 for the validity of this assumption) are 0.5 − 1.7    ×    10²⁴ cm^-2. The gas mass is estimated from the integrated continuum flux density of 85 mJy via (2)where F_ν is the total flux density of dust emission and d is the distance to the source. At a distance of 1 kpc and temperatures of 100–300 K, the mass is 0.1–0.3 M_⊙, significantly less than the stellar mass of VLA3 ( ≈ 10–16 M_⊙, Lada et al. 1984; van der Tak & Menten 2005). If the distance is 3 kpc, the gas mass is modified to 1–3 M_⊙, still less than the scaled stellar mass of  >20 M_⊙. Regardless of the distance, the inferred mass of the surrounding gas around VLA3 is much less than the stellar mass. With the further assumption of a spherical source of 0.8″ in diameter, as measured from the visibilities of VLA3 at 203.4 GHz, H₂ number densities of 7−22    ×    10⁷ cm^-3 are derived for a distance of 1 kpc, and of 2 − 7    ×    10⁷ cm^-3 for a distance of 3 kpc. These densities are lower limits if the source structure is flattened. The inferred gas number density is consistent with the envelope density profile within radii of 400 AU (10^7 − 9 cm^-3) as derived by van der Tak et al. (1999). We summarize these results in Table 4.

3.2. Line emission

Compact emission of HDO, HO and SO₂ is found toward VLA3 (Figs. 3a–c). With the improved angular resolutions of our new observations (natural weighting; 3 mm: , PA 26°; 1.3 mm: , PA 22°), all three molecular images are resolved and show a nearly round shape, unlike the diamond shape with a centrally peaked feature seen in the 203.4-GHz continuum. At 1σ noise level, no significant line emission is found toward other sources in AFGL 2591.

Fig. 3 a)–c) Zero moment maps of HDO, HO and SO2 around VLA3, respectively. The contours are at  − 3, 3, 6, 9...σ for all three maps, where 1σ noise levels are 10, 37 and 65 mJy km s-1 for HDO, HO and SO2, respectively. The green crosses mark the positions of emission peaks. The green ellipses represent FWHM sizes derived from visibility fit. The green filled circles indicate the positions of the additional point sources in the visibility fit of HO and SO2. The natural-weighted clean beam sizes, shown in blue filled ellipses, of HDO and HO/SO2 are , PA 26° and , PA 22°, respectively. d)–l) First moment maps of HDO, HO and SO2 around VLA3, respectively. For all nine maps, the color scales are set to be the same. Dashed lines represent the isovelocity contours with an interval of 0.25 km s-1. m)–o) Molecular emission integrated over the velocity ranges of  −5.35 < VLSR < 0.6 km s-1 (red contours) and  −10.1 < VLSR <  −5.35 km s-1 (blue contours) for HDO, HO and SO2, respectively. Contours are plotted with 10% interval of peak value of each velocity component. The inferred position angles of the velocity gradients across VLA3 are plotted as green dashed lines.

3.2.1. Measurements of positions and sizes

We derive the positions and the sizes of the line emission in the visibility domain (Fig. 4). In the vector-averaged visibility plots, the amplitude of the visibility decreases up to 200 kλ and is more or less flat onward. This feature is similar to what we observed in the continuum visibility at 203.4 GHz, which implies a point source surrounded by an extended source. As a result, we adopted an elliptical Gaussian plus a point source in the visibility fit of H₂¹⁸O and SO₂. For HDO, we used a Gaussian only in the fit since our baselines at 80.6 GHz only extend to 200 kλ. The source sizes derived in this approach are 0.4″–1.0″. The lower end of the result is estimated from HO, which could be due to the large dispersion of the visibilities or reflect a smaller source size traced by this higher excitation line (E_u = 204 K; 47 and 70 K for HDO and SO₂, respectively). The size measurements of the molecular lines are consistent with those derived from the continuum, indicating that the dust and the molecular gas have a similar spatial distribution. The locations of the emission peaks differ by much less than the synthesized beam and are negligible. We summarize the results of the visibility analysis in Table 5.

3.2.2. Velocity gradients across VLA3

The new dataset allows us to spatially resolve the velocity distribution across VLA3. In Fig. 3d–f, velocity gradients in the northeast-southwest direction can be seen with the red-shifted part toward the southwest side and the blue-shifted part toward the northeast side. Apart from the difference in angular resolution, HDO and HO show a similar distribution in isovelocity contours. Detailed investigation of H₂¹⁸O image reveals that the isovelocity contours near the systemic velocity (green part) exhibits a mirrored S shape, while SO₂ shows a normal S-shaped distribution of isovelocity contours. The difference implies that HDO/H₂¹⁸O and SO₂ trace different kinematics of VLA3, where SO₂ may be dominated by the east-west outflow. First moment maps plotted in the velocity ranges of  −5.35 < V_LSR < 0.6 km s^-1 (Fig. 3g–i) and  − 10.1 < V_LSR <  −5.35 km s^-1 (Fig. 3j–l) support our speculation that SO₂ emission is dominated by the outflow since the extreme velocities originate near the source center, while HDO and H₂¹⁸O are less affected by the outflow and presumably dominated by the equatorial motion. We derive the position angles of the velocity gradients across the source seen in Fig. 3d–f by integrating over the velocity ranges of  −5.35 < V_LSR < 0.6 km s^-1 and  − 10.1 < V_LSR <  −5.35 km s^-1 for each line and overplotting them in Fig. 3m–o, respectively. The position angles of velocity gradients are calculated from the emission peaks of the blue- and red-shifted components derived from Gaussian fit. For HDO, H₂¹⁸O and SO₂, the position angles are 54  ±  12°, 40  ±  7° and 33  ±  9°, respectively (Table 5).

Fig. 4 Observed visibilities of HDO, HO and SO2 (filled circles) overplot with model visibilities (grey area). The observed data are vector-averaged corresponding to the circularly symmetric structure. Full model visibilities are plotted in order to show the non-circular symmetric structure of the source.

The orientation of the observed velocity gradients (northeast-southwest) is not perpendicular to the large scale outflow (east-west) (Mitchell et al. 1992; Hasegawa & Mitchell 1995; van der Tak et al. 1999) with a difference of  ~45° in position angle. We rule out the possibility of contamination from the bipolar outflow since the red-shifted emission of the molecular emission lies to the southwest, which is inconsistent with the red lobe of the outflow in the east. High-resolution K-band imaging (Preibisch et al. 2003) shows that the orientation of the large-scale outflow remains unchanged on small scales, ruling out contamination by a small-scale outflow to explain the observed velocity pattern. A more likely interpretation for the observed velocity pattern is rotation in the equatorial plane around AFGL 2591–VLA3. For purely azimuthal motion we would expect the velocity gradient to be perpendicular to the bipolar outflow direction. Section 4 further explores the effect of radial motion in the equatorial plane that affects the orientation of the velocity gradient.

To further quantify the observed velocity gradients, we plot the position-velocity maps with the position angles defined in Fig. 3m–o (Fig. 5). Clearly, in all cases a linear velocity gradient is found, which is very different from the typical Keplerian rotation found in the protoplanetary disks around low-mass stars (e.g., Simon et al. 2000). We derive (projected) velocity gradients of 5.5    ±    0.4, 10.2    ±    1.3 and 12.4    ±    0.5 km s^-1 arcsec^-1 or km s^-1 per 1000 AU (at 1 kpc) for HDO, HO and SO₂, respectively, suggesting that SO₂ is more affected by the east-west outflow while HDO is less affected. To minimize the contamination of the outflow and highlight the motion in the equatorial plane of VLA3, we also derive the velocity gradients by taking the position-velocity cut with PA = 0°. A consistent strength of linear velocity gradient (~15 km s^-1 arcsec^-1) is derived from all three molecules, indicating that the effect of outflow is indeed minimized.

3.2.3. Molecular line profile

Spectra of HDO, HO and SO₂ centered at the mean emission peak (Δα =  −0.15″, Δδ = 0.02″) in a single synthesized beam (, PA 26°) are shown in Fig. 6. A triple-peaked line profile is observed for HDO, while the HO and SO₂ line profiles are smooth Gaussians with weak hints of secondary peaks coincident with those seen in HDO. We use three Gaussians to decompose the HDO line profile (Table 6) and find that Gaussians with FWHM of  ~ 1.20    ±    0.05 km s^-1 at V_LSR of  − 6.89    ±    0.03,  − 5.35    ±    0.34 and  − 3.92    ±    0.26 km s^-1 best reproduce the observed spectrum, although significant residuals remain near  − 8.0 and  − 3.0 km s^-1. The velocity differences between the blue-shifted and systemic components (1.54    ±    0.34 km s^-1), and, the red-shifted and systemic components (1.43    ±    0.43 km s^-1) are comparable, which suggest that the red- and blue-shifted components are likely from the same structure undergoing systematic motion, not just a spatial coincidence of three different gas clumps. To decompose the HO and SO₂ line profiles, we fit three Gaussians fixed at V_LSR of  − 6.89,  − 5.35 and  − 3.92 km s^-1 derived from HDO line profile (Table 6). For an FWHM line width of 1.45 km s^-1, this fit reproduces most of the line profile, with only significant residuals at extreme velocities (near  − 2.5 and  − 8.5 km s^-1 for H₂¹⁸O, and near  − 2.0 and  − 9.0 km s^-1 for SO₂). Based on the line profiles and the observed velocity gradient across the source, we suggest that the velocity structure of VLA3 consists of three components: two exhibiting motions (azimuthal and radial) in the equatorial plane and one static component.

Fig. 5 a)–c) Position-velocity maps of HDO (PA 54°; see Fig. 3m–o), HO (PA 40°) and SO2 (PA 33°). Linear velocity gradients are fitted to the maps and shown as dotted lines. d)–f) Same as a)–c) except the position angle of the position-velocity cut for all three molecules is 0°.

3.2.4. Molecular excitation

We estimate the excitation conditions of VLA3 based on the line ratios of the chemically related molecules HDO and HO. If we assume optically thin emission and local thermodynamic equilibrium (LTE), the total column density at a given excitation temperature can be derived via (3)where W is the integrated intensity in Jy beam^-1 km s^-1, θ_a and θ_b are the FWHM widths of the synthesized beam in arcsec, Q_tot is the partition function, E_u is the upper energy level in K, T_ex is the excitation temperature in K, g_I and g_K are the spin and K degeneracies, respectively, S is the line strength, μ is the dipole moment in Debye, and ν is the rest frequency in GHz. The source filling factor is assumed be unity and background emission is neglected. In other words, the excitation temperature can be expressed as a (nonlinear) function of the column density ratio of HDO and HO, if HDO and HO share the same excitation temperature. Once the excitation temperature is estimated, the column density can be derived via Eq. (3).

Fig. 6 Molecular spectra taken at the mean emission peak after convolving all data to the same angular resolution of , PA 26°. Fitted Gaussians are overplotted in color. Residuals are shown in dots.

The column density ratio N_tot[HO]/N_tot[HDO] can be expressed as (4)For the column density ratio of H₂O and HDO, we adopt the abundance ratio [H₂O]/[HDO] of 1400–2000 for the inner region of the envelope around VLA3 (assuming this ratio can be applied to the observed rotating structure), which is derived based on the “jump” model of 1D radiative transfer calculations (see Table 11 of van der Tak et al. 2006). We assume the column density ratio of HO and H₂O to be 0.002 based on the [¹⁸O]/[¹⁶O] ratio (Wilson & Rood 1994). Therefore, the column density ratio N_tot[HO]/N_tot[HDO] is 2.8–4.0. The total column densities of HDO and HO as a function of T_ex are shown in the upper three panels of Fig. 7, while the column density ratio N_tot[HO]/N_tot[HDO] as a function of T_ex can be seen in the bottom three panels of Fig. 7. From the column density ratio plots, a range of excitation temperatures for a given velocity component can be estimated by finding the intersections of the derived curves of the column density ratio from the new observations (in purple) and the horizontal lines of the adopted column density ratios from Eq. (4) (in green). We summarize the molecular excitation condition of VLA3 in Table 7.

We estimate the excitation temperatures for the blue-shifted, systemic and red-shifted components to be 120–165 K, 155–240 K and 120–165 K, respectively, or 90–290 K, 110–540 K and 90–290 K, respectively, if the error in absolute flux calibration is considered. We find that the blue- and red-shifted components have similar excitation temperatures, implying again that both components are from a single source structure such as the observed equatorial gas component around VLA3. There is also an indication that the systemic component might be warmer. The overall column densities counting all three velocity components (without absolute error in flux density) for HDO, H₂¹⁸O and SO₂ are 1    ×    10¹⁶, 4    ×    10¹⁶ and 1    ×    10¹⁶ cm^-2, respectively. The fractional abundances are estimated to be 5    ×    10^-8, 2    ×    10^-7 and 6    ×    10^-8 for HDO, H₂¹⁸O and SO₂, respectively, adopting the 203.4-GHz continuum peak flux at the resolution of , PA 26° and using T_d = 200 K in the calculation of Eq. (1).

Fig. 7 HO and HDO column density and column density ratio plots. The left, middle and right columns represent the plots for the blue-shifted, systemic and red-shifted components, respectively. In the first row, total column densities at different excitation temperatures are plotted in blue (HO) and red (HDO) thick curves. Their column density density ratios, shown in the second row, are plotted in thick purple curves. The thin curves represent the errors. In addition, the HO over HDO column density ratio, derived from van der Tak et al. (2006) and Wilson & Rood (1994), are plotted as green horizontal lines. The excitation temperatures for a given component can be read out by finding the intersection of the green lines and thick purple lines. The thick vertical dashed lines mark the derived temperatures for ratio 2.8 and 4.0. Here we derive the excitation temperature without considering the error of absolute flux density. Derived temperatures and column densities are summarized in Table 7.

Fig. 8 HDO rotation diagram plotted with the data published by van der Tak et al. (2006). Red open circles are the data observed with JCMT and IRAM 30 m. The green crosses represent the best-fit model from population diagram analysis (Goldsmith & Langer 1999).

We note that the line ratio analysis of HDO and HO discussed above is sensitive to the assumed abundance ratios of HDO/H₂O and ¹⁶O/¹⁸O. For example, if the assumed N_tot[HO]/N_tot[HDO] is smaller (larger) by a factor of two, the corresponding temperatures would be less than 100 K (greater than 250 K) for all three velocity components. Alternatively, if we adopt a temperature of 100 K (the ice evaporation temperature), the derived column density ratio N_tot[HO]/N_tot[HDO] would be less than 7, consistent with van der Tak et al. (2006). Only observations of multiple transitions of HDO and HO can further constrain the excitation conditions of the inner region of VLA3.

Alternatively, the excitation conditions on large scales near VLA3 can be estimated from the HDO single-dish data published by van der Tak et al. (2006) with the population diagram method (Goldsmith & Langer 1999). Including the effects of beam dilution and line opacity, the excitation temperature, total column density and source size are solved self-consistently. The observed data (red open circles) and best fit model (green crosses) are plotted in Fig. 8. The excitation temperature is estimated to be 130  ±  30 K with a column density of  ~(1.0  ±  0.1)  ×  10¹⁴ cm^-2 distributed uniformly in a source with angular size of 17  ±  2″. Line opacities of all the transitions used are estimated to be less than 0.006. The excitation temperatures derived from the single-dish data are consistent with the ones derived from our new interferometric data. The total HDO column density including all three velocity components is about 3  ×  10¹³ cm^-2 if we observe the same structure (roughly 0.8″ − 1.0″) seen by the interferometer with a 17″ beam which is consistent with the value derived from single-dish data if missing flux is considered. As a result, both single-dish and interferometric data imply that HDO is present in both the outer envelope and the inner dense region toward VLA3 with a temperature greater than  ~100 K.

4. Kinematics of the equatorial region of VLA3

4.1. Evidence of rotation and outward radial motions

The new (sub)arcsecond resolution PdBI observations resolve the velocity field in the inner region of AFGL 2591–VLA3. Clear linear velocity gradients in the northeast-southwest direction are observed in HDO, HO and SO₂ lines. Comparison of the distributions of isovelocity contours of these lines reveals that HDO and HO best reflect the kinematics of the equatorial region of VLA3, while SO₂ is more affected by the east-west outflow. In addition, a triple-peaked line profile is observed in HDO. Given the bipolar outflow geometry in the east-west direction (Mitchell et al. 1992; Hasegawa & Mitchell 1995; Preibisch et al. 2003), which defines a north-south equatorial plane, we propose that the observed lines trace rotation in the inner envelope around VLA3. We suggest that radial motions alter the velocity field’s position angle from the expected north-south direction to the observed northeast-southwest direction. Consider an inclined disk-like structure at PA 0°, with rotation and radial motion but without vertical motion (as expected for a model in vertical hydrostatic equilibrium), four possible velocity patterns exist (as shown schematically in Fig. 9), depending on whether the rotation appears clockwise or counterclockwise as seen by the observer, and whether the radial motion are inward or outward. Only the combination of counterclockwise rotation and outward radial motions reproduce the observed velocity pattern (cf. Figs. 3 and 9).

4.2. A toy model

To explore if a kinematical model as suggested in the previous section can reproduce the observed velocity pattern and molecular line spectra, we constructed a toy model and calculate the resulting emission using the axisymmetric radiative transfer code RATRAN (Hogerheijde & van der Tak 2000). Our model includes a flattened “disk-like” structure surrounded by a static spherical envelope and a bipolar outflow cavity (Fig. 10). We adopt a parametrized description of the density, temperature, and velocity field, and assume LTE and optically thin conditions for the line emission. Given the critical densities of the observed transitions (7    ×    10⁴, 3    ×    10⁵ and 5    ×    10⁶ cm^-3 for HDO, HO and SO₂ at 200 K, respectively), and the densities derived in Sect. 3.1.2, the excitation is likely in LTE. We scale the model intensity to match the observed line strengths, which means that we can derive conclusions about the velocity field but not the density or temperature. Our only aim is to show that our simple toy model can reproduce the observed velocity pattern and spectral line shapes. We do not perform any global model optimization in density, temperature and velocity profiles. Instead, for fixed sets of density and temperature profiles, we conduct a simplified optimization to the parametrized velocity field. Our derived parameters are therefore not best-fit solutions but rather indicative values.

Fig. 9 Schematic view of the velocity fields combining radial motion and rotation.

Fig. 10 Gas density (left) and gas temperature (right) distributions of the source in the radiative transfer models, including a disk-like structure, a spherical turbulent envelope, and an outflow cavity.

4.3. Details of the model

The density profiles for the disk-like structure (for simplicity, we use “disk” hereafter in this section) and the envelope are assumed to follow power-law distributions, respectively, as n_disk(r,θ) = An_e0(r/r_d) ^− αsin^f(θ) and n_env(r) = n_e0(r/r_d) ^− α, where n_e0 is the density at radius r_d (the outer radius of the disk), A is the disk-to-envelope density ratio at r_d in the mid-plane, θ is the angle measured from the polar axis, f is the flattening parameter and α is the power-law index. In the disk density profile, regions with densities lower than An_e0 are set to have zero density, resulting in a dumpling-like, flattened structure. The overall density in the model is the sum of n_disk and n_env. For the temperature profiles in the disk and the envelope, we assume power-law distributions, respectively, as T_disk(R) = BT_e0(R/r_d) ^− β and T_env(r) = T_e0(r/r_d) ^− β (i.e., we adopt a vertically isothermal model for the disk) where R denotes the radial direction in cylindrical coordinates, T_e0 is the temperature at r_d, B is the disk-to-envelope temperature ratio, and β is the power-law index. The gas temperature for a given position is the weighted mean of the disk and envelope values: T = (n_diskT_disk + n_envT_env)/(n_disk + n_env). The distance to the source is assumed to be 1 kpc. The inner radius of the model is set to 12.5 AU (half the size of the HII region; van der Tak & Menten 2005) and the outer radius to 400 AU (=r_d, based on the size derived from our 203.4 GHz continuum). In addition, a bipolar outflow cone with an opening angle of 60° is included in the model by setting the density inside this cone to be zero. We describe the velocity field in the disk by three orthogonal components as radial (v_R), vertical (v_z) and azimuthal (v_a). The velocity field in the static envelope is set to be purely turbulent, inspired by the observed spectra showing emission peaks near systemic velocity. The gas velocity for a given position is the weighted mean of the disk and envelope values: v = (n_envv_env + n_diskv_disk)/(n_disk + n_env). The inclination of the source is set to be 30° (van der Tak et al. 2006) with the blue-shifted outflow cone pointed to the west. Model parameters are chosen to ensure that the LTE and optically thin assumptions are valid. We adopt A = 10, n_e0 = 1    ×    10⁷ cm^-3, α = 1.5, f = 5, B = 1, T_e0 = 100 K, β = 1, a molecular fractional abundance of 1 × 10^-13 and a turbulent line width of 0.5 km s^-1. A detailed description of the parametrized velocity field is given in the next section. We scale the line intensity to the observed values, and only analyze the velocity signatures (moment 1 maps and position-velocity maps) and spectral line profiles. We note that this scaling process makes the exact profiles of density and temperature irrelevant.

Fig. 11 Model moment 1 maps of the combined velocity profiles in azimuthal and radial directions. The model position-velocity maps are included in each panel (bottom-left: PA 0°, bottom-right: PA 45°). In all cases, linear velocity gradient can be inferred. However, only Keplerian-like rotation (~r-0.5) plus Hubble-law like radial expansion (~r + 1) can reproduce the observed mirrored S-shape isovelocity contours seen in HO (cf. Fig. 3e).

4.4. Evidence of sub-Keplerian rotation and Hubble-law like radial expansion traced by HDO and HO

The analysis of the position-velocity maps of the observed lines toward VLA3 reveals that the velocity field is characterized by a linear velocity gradient (Fig. 5), which suggests that the equatorial region is undergoing solid-body rotation. However, the analysis demonstrated in Fig. 9 implies that an extra outward radial component is required to reproduce the observed northeast-southwest velocity pattern (Fig. 3d–f). Therefore, a direct link between the observed linear velocity gradient and the inferred solid-body rotation may not be trivial. In our toy model, we consider two types of velocity profiles in the azimuthal and radial directions for the disk-like structure: one is proportional to r^-0.5 and the other follows r ^+ 1. Although all the combinations of the velocity field can reproduce the triple-peaked line signature seen in HDO and the linear velocity gradient seen in the position-velocity maps, only Keplerian-like rotation (v_a ~ r^-0.5) plus Hubble-law like radial expansion (v_R ~ r ^+ 1) can reproduce the mirrored S-shape isovelocity contours observed in HO, which is less affected by the outflow than SO₂ (Fig. 11; cf., Fig. 3e). An additional velocity component in the z direction may redistribute the isovelocity contours to match the observed SO₂ moment 1 map. In our toy model, the outflow component is not included due to the complexity of outflow morphology and velocity field. Therefore, we use HDO and H₂¹⁸O data to constrain the kinematics of the equatorial region of VLA3. To summarize, we adopt the parametrized velocity for the disk-like structure as: (5)where v_R0 is an analog of the Hubble constant, ξ is a constant between 0 and 1, G is the gravitational constant and M _⋆  = 16   M_⊙ (van der Tak & Menten 2005) is the stellar mass.

Fig. 12 Simplified parameter optimization of our kinematic model. The grey scales represent the similarity of model and observed HDO triple-peaked line profile (lower values for better match). The solid curves with values in unit of km s-1 arcsec-1 indicate the measured strength of linear velocity gradient at PA 0° in the grid models. The open circle highlights our best matching parameter set (vR0 = 10 km s-1 arcsec-1, ξ = 0.35).

Fig. 13 (Left) Best matching model spectrum a), moment 1 maps b)–d) and position-velocity map f) for HDO. The color scales and contour levels for moment 1 maps are identical to Fig. 3d–l. The dashed line in the position-velocity map represents a linear velocity gradient with strength of 15 km s-1 arcsec-1. (Right) Best matching model for HO.

To estimate the strength of the two free parameters v_R0 and ξ in our kinematic model, we perform a simplified parameter optimization which utilizes the peak velocities of the HDO triple-peaked line profile (−6.7,  −5.4 and  −3.9 km s^-1 for the blue-shifted, systemic and red-shifted components, respectively), and the measured linear velocity gradient of  ~15 km s^-1 arcsec^-1 at P.A 0° as the constraints. We constructed a grid of models for HDO with v_R0 between 4 and 20 km s^-1 per 1000 AU and ξ between 0.2 and 0.5. For each set of parameters (v_R0, ξ), we compute Δ² = Σ(observed  peak  velocities − model  peak  velocities)², which is smaller for a better match. Figure 12 shows that the models with the lowest Δ² values span a diagonal in (v_R0, ξ)-space, indicating that the model is degenerate. To break the degeneracy, we compute the strength of the linear velocity gradient at PA 0° in our grid model and plot it as the solid curves in Fig. 12. By matching the observed value, we conclude that based on our toy model, the azimuthal motion of the disk-like structure traced by HDO can be characterized by sub-Keplerian rotation (about 1/3 the strength of Keplerian rotation with a stellar mass of 16 M_⊙) and the radial motion is the accelerated expansion with the analogue Hubble constant of 10 km s^-1 per 1000 AU at 1 kpc. The best matching model spectrum, moment 1 maps and position-velocity map for HDO can be found in Fig. 13a–e.

To model the H₂¹⁸O data, we use the best matching solution for HDO (v_R0 = 10 km s^-1 arcsec^-1, ξ = 0.35) as the starting point. Certainly, the turbulent line width of 0.5 km s^-1 is too narrow for H₂¹⁸O since its line profile is Gaussian-like with hints of secondary velocity components. We therefore adopt a turbulent line width of 1.0 km s^-1 in the model. The best matching model for HO can be seen in Fig. 13f–j. Though not perfectly matching the observations, our toy model of kinematics is able to reproduce both HDO and H₂¹⁸O data qualitatively. In our toy model, we found that the kinematics of the proposed disk-like structure around VLA3 is best described by two velocity components: sub-Keplerian rotation and Hubble-law like expansion. We also found that the conclusion of rigid-body rotation across the source judging directly from the linear velocity gradient observed in position-velocity maps may not be made if the observed velocity gradient does not align well with the equatorial plane defined by the outflow.

5. Discussion

5.1. An embedded disk-like structure in AFGL 2591–VLA3?

Our (sub)arcsecond PdBI images of HDO, H₂¹⁸O and SO₂ reveal a clear velocity gradient in the northeast-southwest direction across the source AFGL 2591–VLA3 (Sect. 3.2.2). Together with the shape of the HDO line profile (Sect. 3.2.3), and the symmetric temperature distribution of the blue- and red-velocity components (Sect. 3.2.4), we propose that the source AFGL 2591–VLA3 is surrounded by a disk-like structure undergoing rotation and expansion. The kinematics of the equatorial region is further characterized as sub-Keplerian rotation plus Hubble-law like expansion as shown by radiative transfer modeling of a toy model (Sect. 4). In the models for HDO, a microturbulent line width of  ~0.5 km s^-1 matches the triple-peaked line profile, which appears small for the violent environment of massive star-forming regions (Zinnecker & Yorke 2007). On the other hand, a microturbulent line width of 1.0 km s^-1 is needed to match the almost single-peaked line profiles seen in H₂¹⁸O. We expect that an even larger turbulent line width is required to model SO₂, judging from the observed SO₂ spectrum. This suggests that these molecules trace different regions within the disk-like structure.

We hypothesize that the disk-like structure has a layered distribution with HDO close to the mid-plane, HO in the warm mid-layer and SO₂ in the upper disk layers (Fig. 14). In addition, significant emission from all three molecules originates in a static turbulent envelope surrounding the disk-like structure, especially for H₂¹⁸O and SO₂ where it dominates the line profiles (Table 7). In the disk-like structure, we find temperatures of about 120–165 K, well above the evaporation temperature of water ice (Fraser et al. 2001). The HDO molecules may be formed on grain surfaces (Tielens 1983; Parise et al. 2005), and released into the gas phase when the disk-like structure is heated by shocks, stellar radiation and winds. Compared to HDO, H₂¹⁸O originates from a more turbulent region at larger height with increased turbulence in the disk-like structure. In this picture, the presence of HO but absence of HDO in the mid-layer of the disk-like structure could be due to increased local water formation because of warm gas chemistry as opposed to evaporation of deuterated ices formed at low temperature. SO₂ is a good shock/outflow tracer in massive star-forming regions (Schilke et al. 1997). Judging from the facts that the observed SO₂ line profile has significant emission from line wings, and that a relatively large velocity gradient is measured in SO₂ compared to HDO/H₂¹⁸O, we suggest that SO₂ traces a part of the disk-like structure where the outflow interacts with the upper layers. The concentric velocity patterns seen both in the blue- and red-shifted part of moment 1 maps of SO₂ (Fig. 3i and l) likely support the idea since the most extreme velocities occur in the region close to the star (unlike HDO and H₂¹⁸O) as expected for an outflow origin.

Fig. 14 A cartoon showing the idea of a layered disk-like structure around AFGL 2591–VLA3. HDO is close to the cooler, quiescent mid-plane. HO is distributed in the warm mid-layer and SO2 originates in the upper disk layers.

5.2. Kinematics of the proposed disk-like structure and comparison to other disk sources

Our toy model suggests that the proposed disk-like structure around VLA3 undergoes sub-Keplerian rotation which cannot be directly inferred from the observed position-velocity maps. We speculate that a magnetic field slows the rotation down below pure Keplerian, which may help to transport the angular momentum outward in VLA3. Similar to the formation of low-mass stars, the magnetic field has been suggested to play an important role in the formation of high-mass stars as well (Girart et al. 2009; Beuther et al. 2010). In strongly magnetized clouds during collapse, magnetic braking is thought to regulate the formation of disks around stars (e.g., Basu & Mouschovias 1994; Allen et al. 2003; Mellon & Li 2008). Certainly in our case, the magnetic field does not fully dominate the kinematics of VLA3, otherwise a solid-body like velocity field would be expected in the inner region of the disk-like structure (cf. Galli et al. 2006, Fig. 2). From OH maser observations, Hutawarakorn & Cohen (2005) measure a maser disk with radius about 750 AU and the magnetic field strength of a few mG toward AFGL 2591–VLA3, which may be sufficient to slow down the rotation. Subarcsecond dust polarization measurements of the magnetic field toward VLA3 and hence deriving the mass-to-flux ratio would be helpful to test our hypothesis if the magnetic field indeed regulates the kinematics of the equatorial regions of young high-mass stars.

In addition to sub-Keplerian rotation, we found that Hubble-law like expansion occurs in the equatorial region of VLA3. The overall kinematics is similar to the SiO maser disk around Orion–IRc2 in which Keplerian rotation and decelerated expansion characterize the kinematics (Plambeck et al. 1990). The expanding motions in the disk could be due to the interaction between the stellar radiation and winds with the disk-like structure, providing the forces to push the gas outward (e.g., Hollenbach et al. 1994; Yorke & Welz 1996; Richling & Yorke 1997). On the other hand, the presence of an expanding velocity component in the disk-like structure could be the result of accretion through the mid-plane of the very inner region so that the expansion helps to conserve the angular momentum in the equatorial region. The infrared speckle imaging of AFGL 2591 by Preibisch et al. (2003) suggests that the outflow or wind is variable in time thus the accretion in the disk may not be continuous but episodic. In addition, little free-free emission is detected toward the source, indicating an evolutionary stage before the development of an ultracompact HII region (van der Tak & Menten 2005). As a result, the main accretion phase might have ceased although models show that massive young stars can continue to gain mass through the disk even after an ultracompact HII region has been developed (e.g., Keto 2007). Future high-resolution observations of multiple HDO lines and other dense gas tracers, such as N₂D⁺ and CH₃CN, would be useful to constrain the physical condition and the kinematics of the disk-like structure around VLA3 better and understand its evolutionary stage of high-mass star formation.

Although massive disks showing pure Keplerian rotation traced by various molecules are reported such as the cases in IRAS 20126+4104 (C³⁴S; Cesaroni et al. 2005), AFGL 490 (C¹⁷O; Schreyer et al. 2006) and NGC 7538S (H¹³CN; Sandell et al. 2003), cases showing non-Keplerian rotation were observed as well such as AFGL 2591 (HDO, H₂¹⁸O; this work) and IRAS 18089 − 1732 (NH₃; Beuther & Walsh 2008). In addition, massive disk candidates have a wide range of mass and size distributions derived via various tracers as summarized by Cesaroni et al. (2007). Despite the fact that the detections of massive disks to date are still limited, Keplerian rotation in the disks is typically found toward young massive stars with masses about 10 M_⊙ such as the cases quoted above. For stellar masses about 16 M_⊙, massive disks start to show non-Keplerian motion such as rigid-body rotation reported in the literature. For even higher stellar masses, the rotating structures (toroids) typically show solid-body rotation and are gravitationally unstable (Beltrán et al. 2005; Furuya et al. 2008). The sub-Keplerian rotation discovered in our work suggests that AFGL 2591–VLA3 may be a special case linking transition of velocity field of massive disks from pure Keplerian rotation to solid-body rotation though definitely more new detections of circumstellar disks around high-mass YSOs are required to examine this hypothesis.

5.3. Effect of the uncertainty in distance

The distance to AFGL 2591 is uncertain as discussed by van der Tak et al. (1999) and Rygl et al. (2012). In this paper, we assume a distance of 1 kpc with total luminosity of 20 000 L_⊙. From the radiative transfer modeling of HDO/H₂¹⁸O with our toy model, we found that the azimuthal motion of the proposed disk-like structure is dominated by sub-Keplerian rotation which is about 1/3 in strength of pure Keplerian rotation with a stellar mass of 16 M_⊙, and the radial motion is consistent with Hubble-law like expansion with the analogue Hubble constant of 10 km s^-1 arcsec^-1. For other distances, our results remain unchanged except for the relative strengths of rotation and expansion. In particular, v_R0 scales as 1/d and ξ scales as , so that for d = 0.5 and 3.0 kpc, v_R0 is changed to 20 and 3.3 km s^-1 arcsec^-1, while ξ becomes 0.5 (effective Keplerian mass = 4 M_⊙) and 0.2 (effective Keplerian mass = 0.7 M_⊙). We suggest that the disk-like structure cannot undergo pure Keplerian rotation at any of these distances since the effective Keplerian masses are not consistent with the scaled total luminosities.

6. Conclusions

Our (sub)arcsecond resolution PdBI observations in millimeter continuum and in lines of HDO 1_1,0–1_1,1, H₂¹⁸O 3_1,3–2_2,0 and SO₂ 12_0,12–11_1,11 successfully resolve the inner thousand AU region of the young high-mass star AFGL 2591–VLA3. We summarize our findings below.

• 1.

  From the visibility analysis of millimeter continuum and molecular lines, a compact source with diameter of  ~ 800 AU at 1 kpc is found toward VLA3. Depending on the assumptions, the H₂ column density, number density and mass are estimated to be 1 × 10²⁴ cm^-2, 1 × 10⁸ cm^-3 and 0.2 M_⊙, respectively.

• 2.

  A linear velocity gradient in the northeast-southwest direction is found across the source, which has an offset in PA by  ~ 40° to the equatorial region (PA 0°) defined by the east-west outflow. We interpret this kinematical signature as a combination of rotation and radial expansion in the equatorial region of VLA3.

• 3.

  By introducing a toy model with simplified parameter optimization, we further constrain the rotation to be sub-Keplerian and the radial expansion to be Hubble-law like. This result is independent to the distance of the source.

• 4.

  We speculate that sub-Keplerian rotation in the equatorial region could be effect of magnetic field which turns the velocity profile from pure Keplerian to sub-Keplerian. The expansion could be due to the interaction between the stellar radiation/winds and the gas in the equatorial region of VLA3.

• 5.

  Based on the molecular line profiles, the chemical structure of the proposed disk appears layered, with HDO near the mid-plane, H₂¹⁸O in the mid-layer and SO₂ in the upper layers, judging from their turbulent line widths.

Acknowledgments

The authors appreciate the anonymous referee for the critical comments which help the preparation of the manuscript. The authors are also grateful to the staff of the Plateau de Bure telescope for assisting with the observations, especially Jan Martin Winters at IRAM Grenoble. The research of K.-S.W. at Leiden Observatory is supported through a Nederlandse Onderzoekschool Voor Astronomie (NOVA) Ph.D. grant.

References

All Tables

All Figures

	Fig. 1 Natural-weighted continuum images in the region of AFGL 2591. (Left) 80.6-GHz continuum. Contours are at 3, 5, 10, 20, 30... σ with 1-σ of 0.15 mJy beam-1. Sources within the field following Campbell (1984) and Trinidad et al. (2003) are marked on the map. (Right) 203.4-GHz continuum. Contours are at 3, 5, 10, 30, 50... σ with 1-σ of 0.4 mJy beam-1.
In the text

Fig. 2 Visibility analysis of the continuum at 203.4 GHz. a), c) Show the observed visibilities (filled circles) and model visibilities (grey area). The observed data are vector-averaged for the purpose of better visualization. Only the averaged structure in circular symmetry can be observed in the plot. Full model visibilities are plotted in order to show the non-circular symmetric structure of the source. b), d) Represent a graphical view of the visibility fits of VLA3 shown in a) and c), respectively. The red plus sign is the peak position of the emission. Blue ellipse and filled circle stand for the Gaussian component and the point source.

In the text

Fig. 3 a)–c) Zero moment maps of HDO, HO and SO2 around VLA3, respectively. The contours are at  − 3, 3, 6, 9...σ for all three maps, where 1σ noise levels are 10, 37 and 65 mJy km s-1 for HDO, HO and SO2, respectively. The green crosses mark the positions of emission peaks. The green ellipses represent FWHM sizes derived from visibility fit. The green filled circles indicate the positions of the additional point sources in the visibility fit of HO and SO2. The natural-weighted clean beam sizes, shown in blue filled ellipses, of HDO and HO/SO2 are , PA 26° and , PA 22°, respectively. d)–l) First moment maps of HDO, HO and SO2 around VLA3, respectively. For all nine maps, the color scales are set to be the same. Dashed lines represent the isovelocity contours with an interval of 0.25 km s-1. m)–o) Molecular emission integrated over the velocity ranges of  −5.35 < VLSR < 0.6 km s-1 (red contours) and  −10.1 < VLSR <  −5.35 km s-1 (blue contours) for HDO, HO and SO2, respectively. Contours are plotted with 10% interval of peak value of each velocity component. The inferred position angles of the velocity gradients across VLA3 are plotted as green dashed lines.

In the text

	Fig. 4 Observed visibilities of HDO, HO and SO2 (filled circles) overplot with model visibilities (grey area). The observed data are vector-averaged corresponding to the circularly symmetric structure. Full model visibilities are plotted in order to show the non-circular symmetric structure of the source.
In the text

	Fig. 5 a)–c) Position-velocity maps of HDO (PA 54°; see Fig. 3m–o), HO (PA 40°) and SO2 (PA 33°). Linear velocity gradients are fitted to the maps and shown as dotted lines. d)–f) Same as a)–c) except the position angle of the position-velocity cut for all three molecules is 0°.
In the text

	Fig. 6 Molecular spectra taken at the mean emission peak after convolving all data to the same angular resolution of , PA 26°. Fitted Gaussians are overplotted in color. Residuals are shown in dots.
In the text

Fig. 7 HO and HDO column density and column density ratio plots. The left, middle and right columns represent the plots for the blue-shifted, systemic and red-shifted components, respectively. In the first row, total column densities at different excitation temperatures are plotted in blue (HO) and red (HDO) thick curves. Their column density density ratios, shown in the second row, are plotted in thick purple curves. The thin curves represent the errors. In addition, the HO over HDO column density ratio, derived from van der Tak et al. (2006) and Wilson & Rood (1994), are plotted as green horizontal lines. The excitation temperatures for a given component can be read out by finding the intersection of the green lines and thick purple lines. The thick vertical dashed lines mark the derived temperatures for ratio 2.8 and 4.0. Here we derive the excitation temperature without considering the error of absolute flux density. Derived temperatures and column densities are summarized in Table 7.

In the text

	Fig. 8 HDO rotation diagram plotted with the data published by van der Tak et al. (2006). Red open circles are the data observed with JCMT and IRAM 30 m. The green crosses represent the best-fit model from population diagram analysis (Goldsmith & Langer 1999).
In the text

	Fig. 9 Schematic view of the velocity fields combining radial motion and rotation.
In the text

	Fig. 10 Gas density (left) and gas temperature (right) distributions of the source in the radiative transfer models, including a disk-like structure, a spherical turbulent envelope, and an outflow cavity.
In the text

Fig. 11 Model moment 1 maps of the combined velocity profiles in azimuthal and radial directions. The model position-velocity maps are included in each panel (bottom-left: PA 0°, bottom-right: PA 45°). In all cases, linear velocity gradient can be inferred. However, only Keplerian-like rotation (~r-0.5) plus Hubble-law like radial expansion (~r + 1) can reproduce the observed mirrored S-shape isovelocity contours seen in HO (cf. Fig. 3e).

In the text

Fig. 12 Simplified parameter optimization of our kinematic model. The grey scales represent the similarity of model and observed HDO triple-peaked line profile (lower values for better match). The solid curves with values in unit of km s-1 arcsec-1 indicate the measured strength of linear velocity gradient at PA 0° in the grid models. The open circle highlights our best matching parameter set (vR0 = 10 km s-1 arcsec-1, ξ = 0.35).

In the text

	Fig. 13 (Left) Best matching model spectrum a), moment 1 maps b)–d) and position-velocity map f) for HDO. The color scales and contour levels for moment 1 maps are identical to Fig. 3d–l. The dashed line in the position-velocity map represents a linear velocity gradient with strength of 15 km s-1 arcsec-1. (Right) Best matching model for HO.
In the text

	Fig. 14 A cartoon showing the idea of a layered disk-like structure around AFGL 2591–VLA3. HDO is close to the cooler, quiescent mid-plane. HO is distributed in the warm mid-layer and SO2 originates in the upper disk layers.
In the text