Multi-line Herschel/HIFI observations of water reveal infall motions and chemical segregation around high-mass protostars
Floris van der Tak Russ Shipman Thierry Jacq Fabrice Herpin Jonathan Braine Friedrich Wyrowski




Infall onto high-mass protostars



SRON Netherlands Institute for Space Research, Landleven 12, 9747 AD Groningen, The Netherlands; vdtak@sron.nl  
  Kapteyn Astronomical Institute, University of Groningen, The Netherlands 
  Université de Bordeaux, France  
  Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany  




The physical conditions during high-mass star formation are poorly understood. 
Outflow and infall motions have been detected around massive protostellar objects, but their dependence on mass, luminosity, and age is unclear.
In addition, physical conditions and molecular abundances are often estimated using simple assumptions such as spherical shape and chemical homogeneity, which may limit the accuracy.

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.

We use HIFI maps of the 987 GHz  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  9–8 and  1_11–0_00 lines.  
The high-J  traces the warm central material and redshifted  1_11–0_00 absorption indicates material falling onto the warm core. 
We probe small-scale chemical differentiation by comparing  752 and 987 GHz spectra with those of .  

Our measured radii of the central part of the  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  1_11–0_00 line relative to  9–8.
The inferred infall velocities are 0.6–3.2 , and estimated mass inflow rates range from 7-5 to 2-2 /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  lines in emission. 
The  1_11–0_00 line profiles are remarkably similar to the differences between the  2_02–1_11 and 2_11–2_02 profiles, suggesting that the  line and the  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  line is detectable in the outflow, where no  emission is seen.  

Together, the  absorption and  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.



Multi-line Herschel/HIFI observations of water reveal infall motions and chemical segregation around high-mass protostars 
Herschel is an ESA space observatory with science instruments provided
by European-led Principal Investigator consortia and with important participation from NASA
    F.F.S. van der Tak <ref>,<ref> R.F. Shipman <ref>,<ref> T. Jacq <ref> F. Herpin <ref> J. Braine <ref> F. Wyrowski <ref>

    Submitted 6 July 2018 / Accepted 19 March 2019
=============================================================================================================================================================================================================================================================================================






§ INTRODUCTION





High-mass stars (>8 ) play a key role in the evolution of their host galaxies, but their formation is poorly understood, especially for masses >20 .
The leading models of high-mass star formation involve infall from a dense protostellar core, and accretion onto the protostar via a circumstellar disk <cit.>.
While rotating disks have been detected around young B-type <cit.> and O-type <cit.> 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 <cit.>. 
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 <cit.>.
Possibly both models are valid under different conditions, or combination models need to be developed <cit.>.
To constrain such models, observations of suitable tracers are essential.

Large-scale (∼0.1 pc) infall motions have been detected toward high-mass star-forming regions in ground-based submillimeter-wave molecular emission line maps <cit.>, in redshifted  line absorption at centimeter wavelengths <cit.>, and recently in SOFIA  spectra <cit.>.
Searches for infall in unbiased selections from catalogues of high-mass star-forming regions confirm the ubiquity of such motions <cit.>.

The water molecule appears to be a promising tracer of infall motions in low-mass star-forming regions <cit.>, and <cit.> 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 <cit.>, which have been modeled successfully as infall, using spherical Monte Carlo models <cit.>.
Stronger evidence comes from maps of the luminous mini-starburst region W43 in low-energy  and  lines <cit.>: extended  absorption which is redshifted with respect to the  10-9 emission clearly indicates infall motions.

The  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  2_02–1_11 excited-state emission line <cit.>. 
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  lines towards high-mass protostars concerns their line shape <cit.>. 
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  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  and  lines toward a sample of high-mass protostars to explore the gas distribution and dynamics.
In particular, we compare  line profiles to those of  in order to search for velocity shifts due to infall motions.
Furthermore, we use  maps to measure the sizes of the protostellar envelopes, and to test the assumption of spherical symmetry in previous analyses of pointed spectra.
Section <ref> describes our observations, and
Section <ref> presents the resulting maps and spectra.
Section <ref> compares our derived infall rates with previous observations and with models, and searches for trends with basic source parameters. 
Finally, Section <ref> describes our conclusions. 






§ OBSERVATIONS





 §.§ Source sample



As part of the guaranteed time program WISH (Water In Star-forming regions with Herschel; ), we have selected 19 regions of high-mass star formation for observation in lines of  and its isotopes with the HIFI instrument <cit.> on ESA's Herschel Space Observatory <cit.>.
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 <cit.> for details.
Table <ref> presents the source sample, where distances have been updated following <cit.>, 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 <cit.>.
However, such kinematically derived distances can be off by factors of 2 in either direction; AFGL 2591 and W33A are cases in point <cit.>.
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 <cit.>. 
For W43-Main, two distance estimates exist which are based on VLBI observations of maser parallax (see also ). 
<cit.> report a distance of 4.9 kpc to the W43-Main core, while <cit.> 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.



 §.§ Data acquisition and reduction



Maps of the  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×1140 MHz (1200 ) at a resolution of 1.1 MHz (0.3 ).
Table <ref> presents a detailed observation log including integration times; system temperatures were around 340 K. 
The FWHM beam size at this frequency is 22” <cit.>, 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  1_11–0_00 line at 1101.698 GHz (hereafter 1101 GHz), the  2_11–2_02 line at 752.033 GHz  (hereafter 752 GHz),  the  10–9 line at 1101.34976 GHz, and the  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  and  lines have been observed in the same tuning as the  2_11–2_02 and the  1_11–0_00 lines, respectively, and thus share the same ObsIDs.
Table <ref> lists the integration times of the spectra; system temperatures were around 200 and 390 K for the 752 GHz and ∼1 THz lines, respectively. 
The pointed 987 and 752 GHz spectra have been presented before by <cit.>; the  and  spectra were presented in <cit.>.
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  and  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 <cit.> with further processing performed in the Herschel Interactive Processing Environment (HIPE; ) version 15; 
further analysis was carried out in the CLASS[http://www.iram.fr/IRAMFR/GILDAS] package, version of December 2015 or later. 
Raw antenna temperatures were converted to T_ mb scale using a main beam efficiency of 63% for both frequencies around 1 THz and 64% for the 752 GHz line[https://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 <ref>. 
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  and .




§ RESULTS





 §.§ Line profiles of , , and  



 

Figure <ref> shows the observed velocity profiles of the  1_11–0_00 and  9–8 lines. 
For IRAS 18151, we show the  1_11–0_00 line as the  line is not detected. 
For IRAS 05358, IRAS 16272, and IRAS 18151, the  9–8 line is weak, so we use the  10–9 line to measure velocities.
For the other sources, the data indicate substantial optical depth in the  10–9 line, so we prefer  9–8 as velocity standard.

While the  (or ) lines appear purely in emission for all sources, the  (or ) profiles show absorption, in some cases mixed with emission.
Despite this difference, the peak of the  absorption is seen to lie close to the peak of the  (or ) emission, but at a measureable velocity offset.
In most cases, the  absorption peak is significantly redshifted from the  (or ) emission peak, by 0.6–3.2 .
Table <ref> reports the peak velocities of the  absorption and  (or ) emission, as estimated directly from the HIFI spectra.
We estimate the uncertainty on these velocities to be ≈0.3 . In some cases, no blueshifted absorption or no absorption at all is seen.

The  9–8 line has a relatively high upper level energy (237 K) and critical density (7.75 ), using spectroscopy from <cit.> and collision data from <cit.>, as provided on the LAMDA database <cit.>.
Since in addition the  abundance is likely to be low (∼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  9–8 emission, this argument seems to hold also for the  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  emission may be an evolutionary effect.
The appearance of the  absorption at redshifted velocities thus implies infalling motions in the gas surrounding the dense warm cores seen in  9–8 and/or  10–9 emission.
The velocity difference between the  and  lines indicates approximate infall speeds between 0.6 and 3.2 , although these values represent line-of-sight averages.

For the sources W3 IRS5, W33A, NGC 6334I, and IRAS 05358, the  absorption peak is blueshifted from the  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  9–8 line, for the detection of infall motions. 
The  velocities in Table <ref> differ from the ground-based values (, Table 1; ) by up to 1 , which shows that velocity precision is often limited by source inhomogeneities, rather than by spectral resolution or other instrumental parameters.

Wyrowski et al. (,) have used SOFIA to measure the  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, <cit.> have measured infall through multi-line  line observations with HIFI, and found two infall components with velocities of 2.7 and 5.3 . 
The ground-state  and  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.







 §.§ Maps of 



 

Figures <ref> and <ref>–<ref> show our maps of the  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. <ref> & <ref>–<ref>).
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 ).
We conclude that the bulk of the warm dense gas in the outflow as traced by the  987 GHz line is confined to a small volume (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. <ref>  & <ref>–<ref>.
Table <ref> reports the resulting radii, which have been deconvolved assuming that the source and beam profiles add in quadrature.
The measured sizes of the  emission are ∼2× smaller than the values measured in high-J CO lines with PACS (; Kwon et al, submitted), and 2–3 times smaller than the sizes of the submm/THz dust emission measured with JCMT/SCUBA and APEX/LABOCA <cit.>.
Evidently, the  emission traces warm dense gas close to the protostars.

Comparing the major and minor axis values in Table <ref>, we see that the  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 <ref> compares the observed shape of the  987 GHz emission to the predictions from radiative transfer models, assuming a constant  abundance, following <cit.>.
These predictions are fits to multi-line  (and isotopic) spectra from HIFI, using the physical structure models from <cit.>.
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 <cit.>.

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 (≈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.






§ DISCUSSION






 §.§ Origin of  and  line emission and absorption



Figure <ref> compares the observed  line profiles with those of the  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 (<ref>) and because the map data have higher noise levels.
Remarkably, the  line profile (shown in black) is very similar to the difference between the two  lines (shown in grey).
As found before for the case of W43-MM1 by another method <cit.>,
this close similarity implies that the  absorption originates in warm gas (T  100 K).
Given the upper level energies of the two  lines (101 and 137 K), the bulk of the  absorption must arise in gas with temperatures between ∼100 and ∼140 K.
These temperatures are just above the point where  ice sublimates from dust grains, which is expected to lead to a strong increase in the gas-phase  abundance <cit.>.
The  absorption is unlikely to arise in the cold outer envelope, where the  abundance is too low to create detectable absorption in  (cf. ).
The success of the subtraction procedure shows that the outer envelope does not contribute to the  absorption.

For the sources W3 IRS5, NGC 7538, W33A, AFGL 2591, G29.96, G10.47, and W51N, the subtraction also reproduces  emission features.
Since emission is sensitive to beam filling factors, this similarity is even stronger evidence that the  line originates between the layers where the 752 GHz line is excited and where the 987 GHz line is excited.
In the models by <cit.>, 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  profile (Fig. <ref>), in particular for the line wings.
The scaling factors that best match the observed profiles range from ≈1 for sources with small deconvolved sizes (Table <ref>) to ≈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  line profiles show absorption in the line wings, especially on the blue-shifted side. 
Clearly, the  column density is sufficient to absorb even at velocities only seen in the wings.
The  9–8 spectra show no such high-velocity signals, which implies that the  abundance is enhanced in the high-velocity gas <cit.>.
For example, the  spectrum toward NGC 6334 I(N) shows absorption out to at least 15–20  from line center, which has no counterpart in .
For this source, the integrated  absorption from the envelope (roughly between –6 and +1 , which has a counterpart in  emission) is approximately equal to that in the high-velocity blue wing.
In contrast, the  9–8 line indicates 10× less mass at high velocities, implying an  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  9–8 profiles.
The lack of high-velocity  9–8 emission for most sources is not an excitation effect, as low-J  lines do not show wings either <cit.>.
We conclude that  abundances in high-mass protostellar outflows are 10× higher than in the envelopes.

The  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 <cit.>, who used HIFI maps of the DR21 region in  10–9 and  1_11–0_00 to derive  abundances of ∼10^-10 for the cool outer envelope, ∼10^-8 for the warm inner envelope, and ∼10^-6 for the shocked outflowing gas.

 

 



 §.§ Infall rates and trends



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

    Ṁ_ acc = 4 π R^2   m( H_2)   n( H_2)   | V_ inf |

where m( H_2) is the mass of the  molecule, 
and the absolute value of the infall speed V_ inf is taken from Table <ref>.
Infall motion appears negative as gas is moving toward the center of the reference frame of our models. 
Given the similarity of the  absorption profile with the difference of the  987 and 752 GHz profiles ( <ref>), we adopt the radius of the 120 K point in the envelope models from <cit.> for R, and the density at that radius for n( H_2). 
These radii vary between 800 and 9000 au, and the densities from 75 to 57 .
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 <ref>, right column) are seen to range from ∼7-5 to ∼2-2 /yr.
These values are in reasonable agreement with other observations <cit.> and with theoretical models <cit.>.
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 ∼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 Ṁ tend to cancel each other.

For our subsamples of mid-infrared quiet and –bright HMPOs, <cit.> and <cit.> have made detailed models of the  distribution in the protostellar envelope, including simple step functions for the  abundance in the inner and outer envelope.
In order to fit the line profiles of , , and  in the pointed HIFI spectra, Herpin and Choi had to include radial motions in their models.
Figure <ref> 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 ≈2× 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 , , and age. 
To estimate the relative ages of our sources, we have used the ratio of /, which is straightforward to compute and appears to be a robust age tracer <cit.>.
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 /, proposed as a stability parameter by <cit.>.
Virial masses are calculated following <cit.>, using line widths from Table A.2 of <cit.>.

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 /, but the statistical significance of these trends is small. 
Even the relation with  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 ≈2σ 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.





§ CONCLUSIONS



Based on our measured velocity shifts between  absorption and  emission, infall motions appear to be common in the embedded phase of high-mass star formation, with typical accretion rates of ∼1-4 /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 <cit.>.
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 <cit.>. 
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' <cit.>.

Our data do not allow to discern trends within specific types of sources, nor with protostellar luminosity.
A study of  line profiles toward a large (N∼ 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)[http://www.spica-mission.org/] <cit.> around 2030, and NASA's Origins Space Telescope (OST)[https://asd.gsfc.nasa.gov/firs/] <cit.> around 2040.



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.


aa







§ APPENDIX: MAPS OF ALL SOURCES


 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

