Multi-line Herschel/HIFI observations of water reveal infall motions and chemical segregation around high-mass protostars

TL;DR

This study uses Herschel/HIFI observations to analyze water emission and absorption lines around high-mass protostars, revealing infall motions, chemical segregation, and water abundance variations within the envelopes.

Contribution

It provides new measurements of infall velocities, mass inflow rates, and chemical differentiation in high-mass protostellar envelopes using multi-line Herschel/HIFI data.

Findings

Infall velocities range from 0.6 to 3.2 km/s.

Mass inflow rates are between 7e-5 and 2e-2 solar masses per year.

Water abundance varies with radius, being low in outer regions and very high in outflows.

Abstract

(Abridged) We use HIFI maps of the 987 GHz H2O 2(02)-1(11) emission to measure the sizes and shapes of 19 high-mass protostellar envelopes. To identify infall, we use HIFI spectra of the optically thin C18O 9-8 and H2O-18 1(11)-0(00) lines. The high-J C18O line traces the warm central material and redshifted H2O-18 1(11)-0(00) absorption indicates material falling onto the warm core. We probe small-scale chemical differentiation by comparing H2O 752 and 987 GHz spectra with those of H2O-18. Our measured radii of the central part of the H2O 2(02)-1(11) emission are 30-40% larger than the predictions from spherical envelope models, and axis ratios are <2, which we consider good agreement. For 11 of the 19 sources, we find a significant redshift of the H2O-18 1(11)-0(00) line relative to C18O 9-8. The inferred infall velocities are 0.6-3.2 km/s, and estimated mass inflow rates range from 7e-5 to 2e-2 M0/yr, with the highest mass inflow rates occurring toward the sources with the highest masses, and possibly the youngest ages. The other sources show either expanding motions or H2O-18 lines in emission. The H2O-18 1(11)-0(00) line profiles are remarkably similar to the "differences" between the H2O 2(02)-1(11) and 2(11)-2(02) profiles, suggesting that the H2O-18 line and the H2O 2(02)-1(11) absorption originate just inside the radius where water evaporates from grains, typically 1000-5000 au from the center. In some sources, the H2O-18 line is detectable in the outflow, where no C18O emission is seen. Together, the H2O-18 absorption and C18O emission profiles show that the water abundance around high-mass protostars has at least three levels: low in the cool outer envelope, high within the 100 K radius, and very high in the outflowing gas. Thus, despite the small regions, the combination of lines presented here reveals systematic inflows and chemical information about the outflows.

Tables

Table 1. Table 1: Source sample.

Sourcea	RA (J2000.0)	Dec	$L_{\mathrm{bol}}$	$d$	$M_{\mathrm{env}}$	Distance
	hh mm ss.s	0 ′ ′′	$L_{\odot }$	kpc	$M_{\odot }$	Reference
mid-IR-quiet HMPOsb
IRAS 05358+3543	05 39 13.1	+35 45 50	6.3$\times$103	1.8	142	(1)
IRAS 16272–4837	16 30 58.7	$-$48 43 55	2.4$\times$104	3.4	2170	(1)
NGC 6334I(N)	17 20 55.2	$-$35 45 04	1.1$\times$103	1.3	2237	(5)
W43 MM1	18 47 47.0	$-$01 54 28	1.8$\times$104	5.5	5992	(2)
DR21(OH)	20 39 00.8	+42 22 48	1.3$\times$104	1.5	472	(1)
mid-IR-bright HMPOs
W3 IRS5	02 25 40.6	+62 05 51	1.7$\times$105	2.0	424	(1)
IRAS 18089–1732	18 11 51.5	$-$17 31 29	1.3$\times$104	2.3	172	(1)
W33A	18 14 39.1	$-$17 52 07	4.4$\times$104	2.4	700	(1)
IRAS 18151–1208	18 17 58.0	$-$12 07 27	2.0$\times$104	2.9	153	(1)
AFGL 2591	20 29 24.7	+40 11 19	2.2$\times$105	3.3	363	(1)
Hot Molecular Cores
G327$-$0.6	15 53 08.8	$-$54 37 01	4.4$\times$104	3.1	1804	(6)
NGC 6334I	17 20 53.3	$-$35 47 00	1.5$\times$105	1.3	439	(5)
G29.96$-$0.02	18 46 03.8	$-$02 39 22	2.7$\times$105	5.3	599	(2)
G31.41+0.31	18 47 34.3	$-$01 12 46	8.8$\times$104	4.9	1142	(2)
Ultracompact HII Regions
G5.89$-$0.39 (W28A)	18 00 30.4	$-$24 04 02	5.1$\times$104	1.3	140	(1)
G10.47+0.03	18 08 38.2	$-$19 51 50	8.1$\times$105	8.6	2568	(3)
G34.26+0.15	18 53 18.6	+01 14 58	7.5$\times$104	1.6	421	(4)
W51N-e1	19 23 43.8	+14 30 26	1.1$\times$105	5.4	5079	(7)
NGC 7538-IRS1	23 13 45.3	+61 28 10	1.3$\times$105	2.7	433	(1)

Table 2. Table 2: Observation log for the H 18 2 superscript subscript absent 2 18 {}_{2}^{18} O pointed observations and H 2 O 987 GHz maps.

Source	Species	ObsIDa	$t_{int}$b	rmsc
			[s]	[mK]
IRAS05358	H${}_{}{}^{18}{}_2{}^{}$O	206124, 206126	3566	20
	H2O	204508	5.86	481
IRAS16272	H${}_{}{}^{18}{}_2{}^{}$O	214417, 214419	3455	20
	H2O	203166	5.86	473
NGC6334I(N)	H${}_{}{}^{18}{}_2{}^{}$O	206383	2965	22
	H2O	204523	5.86	481
W43–MM1	H${}_{}{}^{18}{}_2{}^{}$O	191670, 207372	3566	20
	H2O	215899	5.86	524
DR21(OH)	H${}_{}{}^{18}{}_2{}^{}$O	194794, 197974	3566	20
	H2O	210042	7.86	396
W3IRS5	H${}_{}{}^{18}{}_2{}^{}$O	191658, 201591	3566	20
	H2O	203160	5.86	498
IRAS18089	H${}_{}{}^{18}{}_2{}^{}$O	229882, 229883	3455	20
	H2O	218210	5.86	557
W33A	H${}_{}{}^{18}{}_2{}^{}$O	191638, 208086	3566	20
	H2O	215902	5.86	465
IRAS18151	H${}_{}{}^{18}{}_2{}^{}$O	229880, 229881	3455	20
	H2O	218212	5.86	659
AFGL2591	H${}_{}{}^{18}{}_2{}^{}$O	194795, 197973	3566	20
	H2O	210038	5.86	446
G327	H${}_{}{}^{18}{}_2{}^{}$O	214422, 214423, 214425, 214426	3428	21
	H2O	203169	5.86	490
NGC6334I	H${}_{}{}^{18}{}_2{}^{}$O	206385	2965	22
	H2O	204522	5.86	486
G29.96	H${}_{}{}^{18}{}_2{}^{}$O	191668, 191669, 229875, 229876	3700	20
	H2O	207655	5.86	475
G31.41	H${}_{}{}^{18}{}_2{}^{}$O	191671, 191672, 229873, 229874	3700	20
	H2O	207654	5.86	477
G5.89	H${}_{}{}^{18}{}_2{}^{}$O	229888, 229889, 229890, 229891	3148	21
	H2O	218201	5.86	575
G10.47	H${}_{}{}^{18}{}_2{}^{}$O	229884, 229885, 229886, 229887	3148	21
	H2O	218208	5.86	476
G34.26	H${}_{}{}^{18}{}_2{}^{}$O	191673, 191674, 229871, 229872	3700	20
	H2O	207652	5.86	493
W51	H${}_{}{}^{18}{}_2{}^{}$O	194801, 194802, 207384, 207385	3420	20
	H2O	207651	5.86	485
NGC7538IRS1	H${}_{}{}^{18}{}_2{}^{}$O	191663, 191664, 197976	3569	20
	H2O	203161	5.86	479

Table 3. Table 3: Measured velocities and derived infall rates.

Source	$V$(C18O)	$V$(H${}_{}{}^{18}{}_2{}^{}$O)	$V_{\inf }$a	${\dot{M}}_{acc}$ b
	km s-1	km s-1	km s-1	10-3 M⊙/yr
IRAS 05358	–16.0d	–18.5	+2.5	-0-
IRAS 16272	–47.0	–46.3	–0.7	1.28
NGC 6334 I(N)	–3.8	–2.7	–1.1	4.75
W43 MM1	+98.8	+99.4	–0.6	6.06
DR21(OH)	–3.5	–2.2	–1.3	1.86
W3 IRS5	–38.4	–39.9	+1.5	-0-
IRAS 18089	+32.4	+33.8	–1.4	3.28
W33A	+38.0	+33.8	+4.2	-0-
IRAS 18151	+33.6d	+33.4e	–0.2	0.07
AFGL 2591c	–5.5	-0-	-0-	-0-
G327-0.6	–45.0	–43.1	–1.9	4.19
NGC 6334I	–6.5	–7.9	+1.4	-0-
G29.96c	+98.5	-0-	-0-	-0-
G31.41	+97.3	+99.3	–2.0	9.85
G5.89c	+9.3	-0-	-0-	-0-
G10.47	+67.3	+68.7	–1.4	7.15
G34.26	+58.0	+61.2	–3.2	18.9
W51N-e1	+57.3	+58.9	–1.6	14.2
NGC 7538c	–57.6	-0-	-0-	-0-

Table 4. Table 4: Observed and deconvolved source sizes (arcsec).

Sourcea	Major axis	Minor axis	Position angle (degrees)	Deconvolvedb	Modelc	850–870 $\mu$md
IRAS 05358	12.0 (0.5)	9.3 (0.4)	-46 ( 6)	4.4	15.0	30.0
IRAS 16272	15.9 (0.7)	11.8 (0.6)	-34 ( 5)	8.8	18.0	50.0
NGC 6334 I(N)	17.2 (0.5)	12.5 (0.4)	38 ( 3)	11.3	21.6	42.0
W43 MM1	15.8 (1.6)	12.8 (1.2)	-89 (15)	15.5	11.8	27.0
DR21(OH)	12.4 (0.4)	11.1 (0.4)	-8 (11)	5.4	28.8	33.0
W3 IRS5	14.1 (0.3)	12.0 (0.3)	-9 ( 5)	7.8	40.0	57.0
IRAS 18089	11.0 (0.9)	9.5 (0.8)	-37 (21)	-0-	…	17.0
W33A	11.6 (0.7)	10.0 (0.6)	-44 (16)	5.1	…	30.0
AFGL 2591	11.3 (0.6)	9.7 (0.5)	76 (12)	3.8	25.2	25.2
G327-0.6	14.3 (0.5)	10.2 (0.4)	67 ( 4)	5.5	…	24.0
NGC 6334I	16.1 (0.5)	15.0 (0.4)	88 (13)	12.6	…	40.0
G29.96	10.7 (0.3)	9.6 (0.3)	-13 (11)	3.0	…	16.0
G31.41	9.9 (0.9)	7.5 (0.7)	-5 (13)	-0-	…	15.0
G5.89	10.9 (0.2)	10.0 (0.2)	-29 ( 7)	3.9	…	28.0
G10.47	11.9 (0.7)	10.2 (0.6)	-47 (14)	5.5	…	10.0
G34.26	22.4 (1.0)	16.9 (0.8)	-17 ( 4)	16.9	…	25.0
W51N-e1	14.5 (0.7)	10.9 (0.5)	-63 ( 6)	8.2	…	27.0

Equations

$${\dot{M}}_{\mathrm{acc}}=4\pi R^{2}m({\mathrm{H}}_2)n({\mathrm{H}}_2)|V_{\inf }|$$

Full text

*Aims. *We aim to characterize the dust and gas distribution and kinematics of the envelopes of high-mass protostars. In particular, we search for infall motions, abundance variations, and deviations from spherical symmetry, using Herschel data from the WISH program.

Methods. We use HIFI maps of the 987 GHz H2O $2_{02}$–$1_{11}$ emission to measure the sizes and shapes of 19 high-mass protostellar envelopes. To identify infall, we use HIFI spectra of the optically thin C18O 9–8 and H${}_{2}^{18}$O $1_{11}$–$0_{00}$ lines. The high-$J$ C18O traces the warm central material and redshifted H${}_{2}^{18}$O $1_{11}$–$0_{00}$ absorption indicates material falling onto the warm core. We probe small-scale chemical differentiation by comparing H2O 752 and 987 GHz spectra with those of H${}_{2}^{18}$O.

Results. Our measured radii of the central part of the H2O $2_{02}$–$1_{11}$ emission are 30-40% larger than the predictions from spherical envelope models, and axis ratios are $<$2, which we consider good agreement. For 11 of the 19 sources, we find a significant redshift of the H${}_{2}^{18}$O $1_{11}$–$0_{00}$ line relative to C18O 9–8. The inferred infall velocities are 0.6–3.2 km s*-1*, and estimated mass inflow rates range from 7$\times$10*-5* to 2$\times$10*-2* M*⊙/yr, with the highest mass inflow rates occurring toward the sources with the highest masses, and possibly the youngest ages. The other sources show either expanding motions or H${}_{2}^{18}$O lines in emission. The H${}_{2}^{18}$O $1_{11}$–$0_{00}$ line profiles are remarkably similar to the differences between the H2O $2_{02}$–$1_{11}$ and $2_{11}$–$2_{02}$ profiles, suggesting that the H${}_{2}^{18}$O line and the H2O $2_{02}$–$1_{11}$ absorption originate just inside the radius where water evaporates from grains, typically 1000–5000 au from the center. In some sources, the H${}_{2}^{18}$O line is detectable in the outflow, where no C18*O emission is seen.

Conclusions. Together, the H${}_{2}^{18}$O absorption and C18O emission profiles show that the water abundance around high-mass protostars has at least three levels: low in the cool outer envelope, high within the 100 K radius, and very high in the outflowing gas. Thus, despite the small regions, the combination of lines presented here reveals systematic inflows and chemical information about the outflows.

1 Introduction

High-mass stars ($>$8 M*⊙) play a key role in the evolution of their host galaxies, but their formation is poorly understood, especially for masses $>$20 M⊙*. The leading models of high-mass star formation involve infall from a dense protostellar core, and accretion onto the protostar via a circumstellar disk (Tan et al., 2014; Motte et al., 2018). While rotating disks have been detected around young B-type (Sánchez-Monge et al., 2013; Beltrán & de Wit, 2016) and O-type (Johnston et al., 2015; Cesaroni et al., 2017) protostars, the exact manner in (and rate at) which material is gathered from the surroundings is still a matter of debate.

In the ‘monolithic collapse’ model, a massive dense core collapses under its own gravity and forms a (cluster of) protostar(s), much like the low-mass case. This picture is supported by observations of massive collimated outflows from high-mass protostars (Beuther et al., 2002). In the alternative ‘competitive accretion’ model, the accreting protostellar core is replenished from the surroundings. Evidence supporting this model comes for example from observations of extended contracting motions in pre-protocluster regions (Pillai et al., 2011). Possibly both models are valid under different conditions, or combination models need to be developed (Peters et al., 2011). To constrain such models, observations of suitable tracers are essential.

Large-scale ($\sim$0.1 pc) infall motions have been detected toward high-mass star-forming regions in ground-based submillimeter-wave molecular emission line maps (Motte et al., 2003; Peretto et al., 2006), in redshifted NH3 line absorption at centimeter wavelengths (Sollins et al., 2005; Beltrán et al., 2006), and recently in SOFIA NH3 spectra (Wyrowski et al., 2012, 2016). Searches for infall in unbiased selections from catalogues of high-mass star-forming regions confirm the ubiquity of such motions (Fuller et al., 2005; Klaassen & Wilson, 2007; He et al., 2015; Cunningham et al., 2018).

The water molecule appears to be a promising tracer of infall motions in low-mass star-forming regions (Mottram et al., 2013), and San José-García et al. (2016) linked water observations between low- and high-mass star-forming regions. Spectra of low-$J$ line emission toward high-mass objects often exhibit inverse P Cygni profiles (Van der Tak et al., 2013), which have been modeled successfully as infall, using spherical Monte Carlo models (Herpin et al., 2016). Stronger evidence comes from maps of the luminous mini-starburst region W43 in low-energy H2O and H${}_{2}^{18}$O lines (Jacq et al., 2016): extended H${}_{2}^{18}$O absorption which is redshifted with respect to the 13CO 10-9 emission clearly indicates infall motions.

The H${}_{2}^{18}$O $1_{11}$–$0_{00}$ ground-state absorption towards W43-MM1 is remarkable because its shape closely matches that of the central absorption feature in the H2O $2_{02}$–$1_{11}$ excited-state emission line (Jacq et al., 2016). This resemblance strongly suggests that the two lines originate in the same gas, which is curious given their different excitation energies (101 vs 0 K). In order to understand the similarity of these line profiles, this paper explores whether the same effect is seen in other high-mass protostars.

Another puzzle in previous observations of H2O lines towards high-mass protostars concerns their line shape (Van der Tak et al., 2013). The profiles show narrow line cores from the protostellar envelopes, and broad line wings from the outflows, but the wings are much more pronounced at red- than at blueshifted velocities, and often the blueshifted wings are nearly or entirely missing from the profiles. This asymmetry cannot be due to continuum absorption (e.g. by a disk) which would preferentially affect background gas (i.e., receding velocities). Special geometrical configurations may explain individual cases, but not a sample of many sources. One possibility is that the H2O excitation temperature is close to the brightness temperature of the background, so that no net line emission or absorption appears in the spectra. To explore the origin of the asymmetry, this paper explores its dependence on line properties such as excitation energy and critical density.

This paper uses multi-line maps and spectra of H2O and H${}_{2}^{18}$O lines toward a sample of high-mass protostars to explore the gas distribution and dynamics. In particular, we compare H${}_{2}^{18}$O line profiles to those of C18O in order to search for velocity shifts due to infall motions. Furthermore, we use H2O maps to measure the sizes of the protostellar envelopes, and to test the assumption of spherical symmetry in previous analyses of pointed spectra. Section 2 describes our observations, and Section 3 presents the resulting maps and spectra. Section 4 compares our derived infall rates with previous observations and with models, and searches for trends with basic source parameters. Finally, Section 5 describes our conclusions.

2 Observations

2.1 Source sample

As part of the guaranteed time program WISH (Water In Star-forming regions with Herschel; Van Dishoeck et al. 2011), we have selected 19 regions of high-mass star formation for observation in lines of H2O and its isotopes with the HIFI instrument (De Graauw et al., 2010) on ESA’s Herschel Space Observatory (Pilbratt et al., 2010). The sources were selected to cover wide ranges in bolometric luminosity, mid-infrared brightness, and circumstellar mass, and to include regions with hot molecular cores and ultracompact HII regions; see Van der Tak et al. (2013) for details. Table 1 presents the source sample, where distances have been updated following König et al. (2017), and luminosities and masses have been scaled assuming a simple $d^{2}$ dependence.

Most of the updated distances are direct determinations using trigonometric maser parallax observations. The near kinematic distance for G327 seems to be broadly accepted in the recent literature. Only the case of G31.41 is more complicated. The commonly used distance for G31.41 is 7.9 kpc, based on its radial velocity from the Sun and position on the sky, coupled with a Galactic rotation model (Churchwell et al., 1990). However, such kinematically derived distances can be off by factors of ${{}_{>}\atop{{}^{\sim}}}$2 in either direction; AFGL 2591 and W33A are cases in point (Rygl et al., 2012; Immer et al., 2013). Alternatively, G31.41 may be associated with the W43-Main cloud complex, as suggested by position-velocity diagrams of the molecular gas in the surroundings (Nguyen Luong et al., 2011). For W43-Main, two distance estimates exist which are based on VLBI observations of maser parallax (see also Beltrán et al. 2018). Reid et al. (2014) report a distance of 4.9 kpc to the W43-Main core, while Zhang et al. (2014) report distances to 5 maser spots with distances ranging from 6.21 to 4.27 kpc. Given this large spread, we adopt a distance of 4.9 kpc for G31.41 in this paper, and recommend a specific maser parallax study of G31.41 itself.

2.2 Data acquisition and reduction

Maps of the H2O $2_{02}$–$1_{11}$ line at 987.927 GHz (hereafter 987 GHz) were taken with HIFI Band 4a. The maps are 1*′* on the side, and were taken in on-the-fly (OTF) observing mode. The backend was the acousto-optical Wide-Band Spectrometer (WBS) which provides a bandwidth of 4$\times$1140 MHz (1200 km s*-1*) at a resolution of 1.1 MHz (0.3 km s*-1*). Table 2 presents a detailed observation log including integration times; system temperatures were around 340 K. The FWHM beam size at this frequency is 22*′′* (Roelfsema et al., 2012), which corresponds to 0.14–0.92 pc at the distances of our sources. The maps thus cover at least part of the protostellar outflows, while the beam resolves the protostellar envelopes, but not any possible disks.

Spectra of the H${}_{2}^{18}$O $1_{11}$–$0_{00}$ line at 1101.698 GHz (hereafter 1101 GHz), the H2O $2_{11}$–$2_{02}$ line at 752.033 GHz (hereafter 752 GHz), the 13CO 10–9 line at 1101.34976 GHz, and the C18O 9–8 line at 987.560 GHz were obtained toward the same sources with HIFI, using the Double Beam Switch observing mode with a chopper throw of 3*′. The C18O and 13CO lines have been observed in the same tuning as the H2O $2_{11}$–$2_{02}$ and the H${}_{2}^{18}$O $1_{11}$–$0_{00}$ lines, respectively, and thus share the same ObsIDs. Table 2 lists the integration times of the spectra; system temperatures were around 200 and 390 K for the 752 GHz and $\sim$1 THz lines, respectively. The pointed 987 and 752 GHz spectra have been presented before by San José-García et al. (2016); the 13CO and C18O spectra were presented in San José-García et al. (2013). The DBS spectra have higher noise per second of integration than the maps at the same frequency, which represents the noise ‘penalty’ to be paid for stabilizing the system by differencing two reference positions in the DBS observing mode. For the H${}_{2}^{18}$O and C18O lines, the beam size of 20–22′′* is very similar to that of the 987 GHz maps, which permits a direct comparison of the results. The beam size of the 752 GHz observations is 28*′′*.

The data are Herschel/HIFI standard products (Shipman et al., 2017) with further processing performed in the Herschel Interactive Processing Environment (HIPE; Ott 2010) version 15; further analysis was carried out in the CLASS333http://www.iram.fr/IRAMFR/GILDAS package, version of December 2015 or later. Raw antenna temperatures were converted to $T_{\rm mb}$ scale using a main beam efficiency of 63% for both frequencies around 1 THz and 64% for the 752 GHz line444https://www.cosmos.esa.int/web/herschel/ legacy-documentation-hifi-level-2, and linear baselines were subtracted. After inspection, the data from the two polarization channels were averaged to obtain the rms noise levels reported in Table 2. The absolute calibration uncertainty of HIFI Bands 3 and 4 is estimated to be 10-15%, but the relative calibration between lines in the same spectrum should be much better, which is relevant for C18O and 13CO.

3 Results

3.1 Line profiles of H${}_{2}^{18}$O, 13CO, and C18O

Figure 1 shows the observed velocity profiles of the H${}_{2}^{18}$O $1_{11}$–$0_{00}$ and C18O 9–8 lines. For IRAS 18151, we show the H2O $1_{11}$–$0_{00}$ line as the H${}_{2}^{18}$O line is not detected. For IRAS 05358, IRAS 16272, and IRAS 18151, the C18O 9–8 line is weak, so we use the 13CO 10–9 line to measure velocities. For the other sources, the data indicate substantial optical depth in the 13CO 10–9 line, so we prefer C18O 9–8 as velocity standard.

While the C18O (or 13CO) lines appear purely in emission for all sources, the H${}_{2}^{18}$O (or H2O) profiles show absorption, in some cases mixed with emission. Despite this difference, the peak of the H${}_{2}^{18}$O absorption is seen to lie close to the peak of the C18O (or 13CO) emission, but at a measureable velocity offset. In most cases, the H${}_{2}^{18}$O absorption peak is significantly redshifted from the C18O (or 13CO) emission peak, by 0.6–3.2 km s*-1*. Table 3 reports the peak velocities of the H${}_{2}^{18}$O absorption and C18O (or 13CO) emission, as estimated directly from the HIFI spectra. We estimate the uncertainty on these velocities to be $\approx$0.3 km s*-1*. In some cases, no blueshifted absorption or no absorption at all is seen.

The C18O 9–8 line has a relatively high upper level energy (237 K) and critical density (7.7$\times$105 cm*-3*), using spectroscopy from Endres et al. (2016) and collision data from Yang et al. (2010), as provided on the LAMDA database (Schöier et al., 2005). Since in addition the C18O abundance is likely to be low ($\sim$10*-6*), this line should be an optically thin tracer of the warm dense gas close to the central protostar. For the 3 sources with weak C18O 9–8 emission, this argument seems to hold also for the 13CO 10–9 line, presumably due to a low envelope mass. These 3 sources are not the lowest-luminosity cases in our study, so the low envelope mass and weak C18O emission may be an evolutionary effect. The appearance of the H${}_{2}^{18}$O absorption at redshifted velocities thus implies infalling motions in the gas surrounding the dense warm cores seen in C18O 9–8 and/or 13CO 10–9 emission. The velocity difference between the C18O and H${}_{2}^{18}$O lines indicates approximate infall speeds between 0.6 and 3.2 km s*-1*, although these values represent line-of-sight averages.

For the sources W3 IRS5, W33A, NGC 6334I, and IRAS 05358, the H${}_{2}^{18}$O absorption peak is blueshifted from the C18O emission peak, suggesting expanding motions. The line profiles toward G5.89 and G10.47 are complex, with a mixture of infall and expansion. These two sources are not included in the analysis below.

We emphasize the importance of using a precise velocity standard, in this case the C18O 9–8 line, for the detection of infall motions. The C18O velocities in Table 3 differ from the ground-based values (Van der Tak et al. 2013, Table 1; Van Dishoeck et al. 2011) by up to 1 km s*-1*, which shows that velocity precision is often limited by source inhomogeneities, rather than by spectral resolution or other instrumental parameters.

Wyrowski et al. (2012,2016) have used SOFIA to measure the NH3 $3_{2}^{+}$-$2_{2}^{-}$ line at 1810.379 GHz toward several of our sources. They report redshifted absorption toward W43 MM1, G327, G31.41, and G34.26, implying infall, and blueshifted absorption toward W33A and G5.89, implying expansion. These results agree qualitatively with ours, and their measured velocities are similar to those reported here.

Toward G34.26, Hajigholi et al. (2016) have measured infall through multi-line NH3 line observations with HIFI, and found two infall components with velocities of 2.7 and 5.3 km s*-1*. The ground-state H${}_{2}^{18}$O and NH3 lines presented here and by Wyrowski et al. only probe the lower-velocity of these components, which may mean that the higher-velocity component mostly arises in very warm and dense gas in close proximity to the protostar. This result is consistent with a scenario where the infall velocity of the gas increases as it approaches the protostar.

3.2 Maps of H2O

Figures 2 and 6–22 show our maps of the H2O 987 GHz line emission. The greyscale and white contours denote the line core, while the blue and red contours correspond to the blue- and redshifted line wings (see the caption for details). The emission is seen to be compact (except G31.41), mildly elongated, and not to depend much on velocity interval (except NGC 6334I). For IRAS 18151, the emission is too weak to assess its morphology. The map of NGC 7538 is not shown as it suffers from mispointing, so that only limits on the emission size and shape can be obtained.

The observed morphology of the 987 GHz emission does not appear to depend much on velocity interval (Figs. 2 & 6–22). This contrasts with the low-$J$ CO emission from our sources, which shows a clear bipolar morphology, especially at velocities away from line center (see references in Van Dishoeck et al. 2011). We conclude that the bulk of the warm dense gas in the outflow as traced by the H2O 987 GHz line is confined to a small volume (${{}_{<}\atop{{}^{\sim}}}$20*′′*) from the source, unlike the outflow gas at lower temperature and density traced by low-$J$ CO lines.

We have measured the size of the 987 GHz emission by fitting a two-dimensional Gaussian plus a background offset to the images in Figs. 2 & 6–22. Table 4 reports the resulting radii, which have been deconvolved assuming that the source and beam profiles add in quadrature. The measured sizes of the H2O emission are $\sim$2$\times$ smaller than the values measured in high-$J$ CO lines with PACS (Karska et al. 2014; Kwon et al, submitted), and 2–3 times smaller than the sizes of the submm dust emission measured with JCMT/SCUBA and APEX/LABOCA (Van der Tak et al., 2013). Evidently, the H2O emission traces warm dense gas close to the protostars.

Comparing the major and minor axis values in Table 4, we see that the H2O emission is close to spherical in most cases, with axis ratios between 1.1 and 1.4. We conclude that protostellar envelopes dominate the emission, without any evidence for flattening or elongation due to rotation or bipolar outflows.

Table 4 compares the observed shape of the H2O 987 GHz emission to the predictions from radiative transfer models, assuming a constant H2O abundance, following Herpin et al. (2016). These predictions are fits to multi-line H2O (and isotopic) spectra from HIFI, using the physical structure models from Van der Tak et al. (2013). The predicted size is seen to be 30–40% larger than the observed size for most sources, which we consider good agreement given the simplifying assumption of spherical symmetry in the models. Only for the sources W3 IRS5 and W43 MM1, the predicted size is 2–4 times smaller than the observed size. As with the axis ratios, this may be due to outflows contributing to the emission. Furthermore, the models assume a single central source, whereas interferometric images of our objects often show multiple cores at the center (e.g., Hunter et al., 2014; Brogan et al., 2016; Izquierdo et al., 2018).

The line intensities in the maps are typically 70-80% of the values reported from pointed observations at the same position. This difference is as expected from the 4% larger beam size due to the OTF observing mode and the spatial regridding, assuming a small emitting area. Only for IRAS 18151 and IRAS 18089, the map intensities are substantially lower ($\approx$40% of the pointed observations) for unknown reasons. In such cases, the pointed observations are more reliable, since their calibration is more thorough, with multiple references and longer integrations. We conclude that mapping modes are useful to measure source sizes, but usually underestimate line intensities, sometimes substantially.

4 Discussion

4.1 Origin of H2O and H${}_{2}^{18}$O line emission and absorption

Figure 3 compares the observed H${}_{2}^{18}$O line profiles with those of the H2O 987 and 752 GHz lines. For the 987 GHz line, we use the pointed observations rather than convolving the map data, because of the calibration issue with the maps (§3.2) and because the map data have higher noise levels. Remarkably, the H${}_{2}^{18}$O line profile (shown in black) is very similar to the difference between the two H2O lines (shown in grey). As found before for the case of W43-MM1 by another method (Jacq et al., 2016), this close similarity implies that the H${}_{2}^{18}$O absorption originates in warm gas ($T{{}_{>}\atop{{}^{\sim}}}100$ K). Given the upper level energies of the two H2O lines (101 and 137 K), the bulk of the H${}_{2}^{18}$O absorption must arise in gas with temperatures between $\sim$100 and $\sim$140 K. These temperatures are just above the point where H2O ice sublimates from dust grains, which is expected to lead to a strong increase in the gas-phase H2O abundance (Boogert et al., 2015). The H${}_{2}^{18}$O absorption is unlikely to arise in the cold outer envelope, where the H2O abundance is too low to create detectable absorption in H${}_{2}^{18}$O (cf. Shipman et al. 2014). The success of the subtraction procedure shows that the outer envelope does not contribute to the H${}_{2}^{18}$O absorption.

For the sources W3 IRS5, NGC 7538, W33A, AFGL 2591, G29.96, G10.47, and W51N, the subtraction also reproduces H${}_{2}^{18}$O emission features. Since emission is sensitive to beam filling factors, this similarity is even stronger evidence that the H${}_{2}^{18}$O line originates between the layers where the 752 GHz line is excited and where the 987 GHz line is excited. In the models by Van der Tak et al. (2013), this zone occurs typically at radii of 1000-5000 AU, depending on the luminosity of the source. This region is small enough that it is often difficult to observe (e.g., 2–10*′′* diameter at a distance of 1 kpc).

In some cases, scaling the 752 GHz profile before subtracting it from the 987 GHz line profile improves the match of the difference to the H${}_{2}^{18}$O profile (Fig. 4), in particular for the line wings. The scaling factors that best match the observed profiles range from $\approx$1 for sources with small deconvolved sizes (Table 4) to $\approx$1.8 for the most extended sources. These values are just as expected from beam size differences between the 987 and 752 GHz spectra, assuming equal excitation temperatures. There may be other pairs of lines whose differences enable us to probe specific layers of the protostellar cores.

Toward several of the more massive sources, the H${}_{2}^{18}$O line profiles show absorption in the line wings, especially on the blue-shifted side. Clearly, the H${}_{2}^{18}$O column density is sufficient to absorb even at velocities only seen in the wings. The C18O 9–8 spectra show no such high-velocity signals, which implies that the H2O abundance is enhanced in the high-velocity gas (Herpin et al., 2016). For example, the H${}_{2}^{18}$O spectrum toward NGC 6334 I(N) shows absorption out to at least 15–20 km s*-1* from line center, which has no counterpart in C18O. For this source, the integrated H${}_{2}^{18}$O absorption from the envelope (roughly between –6 and +1 km s*-1*, which has a counterpart in C18O emission) is approximately equal to that in the high-velocity blue wing. In contrast, the C18O 9–8 line indicates ${{}_{>}\atop{{}^{\sim}}}$10$\times$ less mass at high velocities, implying an H2O abundance enhancement by more than an order of magnitude.

Similar conclusions hold for the other sources, except for G5.89 and G34.26 where weak wings are seen on the C18O 9–8 profiles. The lack of high-velocity C18O 9–8 emission for most sources is not an excitation effect, as low-$J$ C18O lines do not show wings either (Hatchell et al., 1998; Watson et al., 2003; Gibb et al., 2004; Thomas & Fuller, 2007). We conclude that H2O abundances in high-mass protostellar outflows are ${{}_{>}\atop{{}^{\sim}}}$10$\times$ higher than in the envelopes.

The H2O abundance in these sources thus appears to have at least 3 levels: low in the outer envelope, high in the inner envelope, and very high in the outflow. This is in line with the work of Van der Tak et al. (2010), who used HIFI maps of the DR21 region in 13CO 10–9 and H2O $1_{11}$–$0_{00}$ to derive H2O abundances of $\sim$10*-10* for the cool outer envelope, $\sim$10*-8* for the warm inner envelope, and $\sim$10*-6* for the shocked outflowing gas.

4.2 Infall rates and trends

The rightmost column of Table 3 gives estimates of the infall rates onto our sources. These were calculated using

[TABLE]

where $m({\rm H}_{2})$ is the mass of the H2 molecule, and the absolute value of the infall speed $V_{\rm inf}$ is taken from Table 3. Infall motion appears negative as gas is moving toward the center of the reference frame of our models. Given the similarity of the H${}_{2}^{18}$O absorption profile with the difference of the H2O 987 and 752 GHz profiles (§ 4.1), we adopt the radius of the 120 K point in the envelope models from Van der Tak et al. (2013) for $R$, and the density at that radius for $n({\rm H}_{2})$. These radii vary between 800 and 9000 au, and the densities from 7$\times$105 to 5$\times$107 cm*-3*. Our observed (deconvolved) sizes agree well (within a factor of 2) with the upper end of this range, except for W43 MM1, G10.47, and W51N, where the observed values are larger.

The resulting infall rates (Table 3, right column) are seen to range from $\sim$7$\times$10*-5* to $\sim$2$\times$10*-2* M*⊙*/yr. These values are in reasonable agreement with other observations (e.g., König et al., 2017) and with theoretical models (Tan et al., 2014; Motte et al., 2018). They should be considered order of magnitude estimates, because of our simplified treatment assuming spherical symmetry. The observational uncertainty through the measured line velocities is only a $\sim$30% effect. The derived infall rates depend only weakly on the adopted radius and density: the envelopes of our sources have density profiles which drop off approximately as $R^{-2}$, so that the effects of $R$ and $n$ on $\dot{M}$ tend to cancel each other.

For our subsamples of mid-infrared quiet and –bright HMPOs, Herpin et al. (2016) and Choi (2015) have made detailed models of the H2O distribution in the protostellar envelope, including simple step functions for the H2O abundance in the inner and outer envelope. In order to fit the line profiles of H2O, H${}_{2}^{18}$O, and H${}_{2}^{17}$O in the pointed HIFI spectra, Herpin and Choi had to include radial motions in their models. Figure 5 compares their derived infall velocities to the values found here; the dashed line indicates 1:1 correspondence. The two types of estimates of the radial velocity are seen to agree qualitatively, both in the sign of the velocity (infall or outflow) and in its magnitude. The simple estimates of the inflow / expansion velocity are on average $\approx$2$\times$ lower than those from the detailed models, although for four sources, they are actually larger. We consider this agreement as reasonable, given the differences between the two approaches.

Since our source sample covers a range of luminosities, envelope masses, and evolutionary stages, we have investigated if our derived source sizes, infall velocities, and infall rates show any trends with $L_{\rm bol}$, $M_{\rm env}$, and age. To estimate the relative ages of our sources, we have used the ratio of $L_{\rm bol}$/$M_{\rm env}$, which is straightforward to compute and appears to be a robust age tracer (Molinari et al., 2016). In addition, we use the presence of hot molecular cores and/or ultracompact HII regions as a sign of a relatively evolved stage. Third, we have looked for trends with the virial mass and the ratio $M_{\rm vir}$/$M_{\rm env}$, proposed as a stability parameter by König et al. (2017). Virial masses are calculated following Giannetti et al. (2014), using line widths from Table A.2 of Van der Tak et al. (2013).

The only significant trend that we find is between the linear sizes of our sources with their virial masses. Since virial mass depends on size, this trend probably just means that the line width is similar for all sources. In addition, the infall rates seem to increase with virial mass and with the evolutionary indicator $L_{\rm bol}$/$M_{\rm env}$, but the statistical significance of these trends is small. Even the relation with $M_{\rm vir}$ has a Pearson correlation coefficient of only $r=$0.58. For a sample size of $N$=11, this $r$-value corresponds to a probability of false correlation of $p=$6%, i.e. a $\approx$2$\sigma$ significance. We conclude that the accretion rates may increase with circumstellar mass and with evolutionary stage, but that larger source samples are required to test these claims.

5 Conclusions

Based on our measured velocity shifts between H${}_{2}^{18}$O absorption and C18O emission, infall motions appear to be common in the embedded phase of high-mass star formation, with typical accretion rates of $\sim$1$\times$10*-4* M*⊙*/yr. We find a tentative trend that the highest accretion rates occur for the most massive sources, which is globally consistent with current models of high-mass star formation (Tan et al., 2014; Motte et al., 2018). Our data do not allow to distinguish between such models, though.

In addition, the accretion rates may increase with age, unlike in the low-mass case, where accretion rates drop from the Class 0 to the Class III stage, and are highly episodic (Dunham et al., 2014). Signs of episodic accretion, which is well established in the low-mass case, have recently been reported for a high-mass star, in the form of mid-infrared variability suggesting accretion ‘bursts’ (Caratti o Garatti et al., 2017).

Our data do not allow to discern trends within specific types of sources, nor with protostellar luminosity. A study of H2O line profiles toward a large ($N\sim 100$) sample is needed to distinguish such trends and to search for episodic behaviour. Data from the Herschel open time programs by Bontemps and Wyrowski may be suitable for this purpose. In the future, such studies will be possible with ESA’s SPace Infrared telescope for Cosmology and Astrophysics (SPICA)777http://www.spica-mission.org/ (Roelfsema et al., 2018; Van der Tak et al., 2018) around 2030, and NASA’s Origins Space Telescope (OST)888https://asd.gsfc.nasa.gov/firs/ (Battersby et al., 2018) around 2040.

Acknowledgements.

This paper is dedicated to the memory of Malcolm Walmsley, who passed away on 1 May 2017 at the age of 75. We remember Malcolm as a great source of inspiration, and we will miss his sharp insight and kind manner. The authors thank the WISH team led by Ewine van Dishoeck for inspiring discussions, and the anonymous referee for useful comments on the manuscript. This research has used the following databases: ADS, CDMS, JPL, and LAMDA. HIFI was designed and built by a consortium of institutes and university departments from across Europe, Canada and the US under the leadership of SRON Netherlands Institute for Space Research, Groningen, The Netherlands with major contributions from Germany, France and the US. Consortium members are: Canada: CSA, U.Waterloo; France: CESR, LAB, LERMA, IRAM; Germany: KOSMA, MPIfR, MPS; Ireland, NUI Maynooth; Italy: ASI, IFSI-INAF, Arcetri-INAF; Netherlands: SRON, TUD; Poland: CAMK, CBK; Spain: Observatorio Astronómico Nacional (IGN), Centro de Astrobiología (CSIC-INTA); Sweden: Chalmers University of Technology - MC2, RSS & GARD, Onsala Space Observatory, Swedish National Space Board, Stockholm University - Stockholm Observatory; Switzerland: ETH Zürich, FHNW; USA: Caltech, JPL, NHSC.

Appendix: Maps of all sources