© J. C. Mottram et al. 2020

1 Introduction

Star formation is an inherently multi-scale phenomenon, with different processes having the potential to affect, guide, or govern its outcome from the scales of the entire Galaxy (for example Galactic gravitational potential, spiral arms, magnetic fields, metallicity) down to the scale of individual protostars (for example turbulence, magnetic fields, fragmentation, rotation, infall, outflow), both for individual stars and for stellar populations in general (for recent reviews see for example Zinnecker & Yorke 2007; Dobbs et al. 2014; André et al. 2014; Offner et al. 2014; Padoan et al. 2014; Li et al. 2014a,b; Molinari et al. 2014; Tan et al. 2014; Krumholz et al. 2014). As such, one over-arching task within star formation research is to understand what processes dominate on which scales, and what minimum range of scales and set of processes must be included and understood in order to predict the outcome of star formation.

Observational studies have noted the hierarchical density structure of star formation regions for some time, and indeed have made use of this to isolate and categorise structures on different size and density scales, ranging from giant molecular associations and giant molecular clouds on the largest scales, through clouds, clumps, filaments, and cores, to embedded discs on the smallest scales (see Table 1 for details and references). However, how the physical and chemical processes on each scale impact those on the scales above and below is still poorly quantified or understood.

Understanding the interplay between different scales and processes will give constraints to both global properties, such as the star formation timescale and efficiency, and will reveal what key factors determine the outcome of star formation. This is a challenging task given the range of temperatures, densities, scales, and processes that are relevant for star formation in general, and particularly in the complex regions where high-mass stars (≳8 M_⊙) are formed. Continuum observations may reveal powering sources or emit from regions of a certain temperature, depending on the wavelength of the observations, but they do not provide kinematic information and may include contributions from multiple distinct regions along a given line-of-sight. Velocity-resolved line observations of molecular, atomic, or ionised species can reveal gas motions, but each transition preferentially emits under certain conditions. Furthermore, abundance variations due to chemical processes and physical conditions can lead to two otherwise similar regions displaying different emission profiles. Finally, the distances to high-mass star formation regions are typically ≳2 kpc, leading to trade-offs between resolution to capture small-scale processes and both field-of-view and sensitivity to the larger scale emission needed to link sites of star formation with their neighbours and environment.

These requirements are best catered for using mm-wavelength interferometric observations of multiple molecular line species and the thermal dust continuum, ideally with complementary continuum data at other wavelengths. The time requirements to observe a statistically useful number of regions are not insignificant, and so most studies to date have focused either on the large-scale cloud and/or clump properties using single-dish observations (Motte et al. 2007; Rathborne et al. 2009; Urquhart et al. 2018) or on detailed studies of core and disc scales using interferometric observations alone (for example Beuther et al. 2012; Johnston et al. 2015; Feng et al. 2016; Ilee et al. 2016; Cesaroni et al. 2017).

The IRAM Northern Extended Millimeter Array (NOEMA) large programme named CORE: “Fragmentation and disc formation in high-mass star formation” (Beuther et al. 2018) is specifically designed to avoid such biases by observing 20 regions of high-mass star formation at 1.4 mm with three interferometric array configurations and IRAM 30 m single-dish telescope observations, making the data sensitive to emission from ~0.3−20′′ (600−40 000 AU or 0.003−0.2 pc at 2 kpc) scales. The principle scientific motivations of the CORE large programme are to study fragmentation, kinematics (for example infall, outflow, rotation), and chemistry from clump to disc scales, and to explore how these larger-scale processes and properties impact the existence and properties of discs around young embedded high-mass stars. A more extensive overview of the survey and a first analysis of the fragmentation in the targeted regions is presented in Beuther et al. (2018), while a detailed analysis of the disc properties in one case-study region, containing W3 H₂O and OH, is presented in Ahmadi et al. (2018). The following paper presents a first case-study of the kinematics and chemical diversity from clump to disc scales, for which we select the target W3 IRS4 in the W3 star formation region from the CORE survey sample.

W3 is part of the Giant Molecular Association (GMA) encosing W3, W4 and W5, which was first identified in the radio by Westerhout (1958), and is at a distance of 2.0 kpc based on maser parallax measurements (Hachisuka et al. 2006; Xu et al. 2006), placing it in the Perseus spiral arm (see Megeath et al. 2008, for a detailed review of the structure of the region). This GMA is one of the largest and most active sites of ongoing star formation in the Outer Galaxy.

Figure 1 presents a three-colour overview of W3, as well as zoom-ins on W3 IRS4 in order to provide environmental context. The three brightest far-IR sources in the W3 star forming region (see for example Rivera-Ingraham et al. 2013) are W3 OH/H₂O, W3 IRS4, and W3 IRS5 (see left panel of Fig. 1), which all host young massive stars that are still accreting material. There are also a number of H II regions of various sizes (for example Tieftrunk et al. 1997, see also the middle panel of Fig. 1) and a number of optically visible main-sequence OB stars (for example Oey et al. 2005; Bik et al. 2012). Due to the high infrared flux of W3 IRS4, originally identified in the near-IR by Wynn-Williams et al. (1972), it has been extensively studied at high spatial resolution in the infrared (see for example Bik et al. 2012) and with interferometers in the radio, mm, and sub-mm (for example Tieftrunk et al. 1997, 1998; Wang et al. 2012b). Similarly, W3 IRS5 was studied with high-resolution interferometric observations at 1.4 and 3.4 mm by Rodón et al. (2008), while W3 OH/H₂O is also part of the CORE sample, with a detailed analysis of these observations presented in Ahmadi et al. (2018).

High spatial-resolution studies of the W3 IRS4 region have shown that the region contains a number of sources and structures, as indicated on a map of the NOEMA 1.4 mm continuum emission in Fig. 2. At near and mid-infrared (IR) wavelengths, the emission is dominated by IRS4 itself (see Fig. 1), sometimes also referred to as W3 Ca, which is a massive young stellar object (MYSO) or hyper-compact (HC) H II region and the brightest continuum peak at 1.4 mm. This source exhibits a relatively rich spectrum, consistent with hot core emission (for example Wang et al. 2012b). A more evolved ultra-compact (UC) H II region, W3 C, lies less than 3′′ (6000 AU or 0.03 pc) away to the south-west, and is detectable from the mid-IR to the radio (Tieftrunk et al. 1997). The combined bolometric luminosity of W3 IRS4 and W3 C (the sources lie within one beam in the available far-IR data) is 4.5 × 10⁴ L_⊙ from SED fitting¹ using the method outlined by Mottram et al. (2011) and updated to include far-IR and sub-mm continuum fluxes from the Herschel Hi-GAL survey (Molinari et al. 2016). The 1.4 mm continuum emission also shows part of a filament that extends much further (≳40′′, 80 000 AU or0.4 pc) to the NE and N in VLA NH₃ observations(Tieftrunk et al. 1998). This filament currently ends at the location of IRS4, though its projected direction is also broadly consistent with the location of the UC HII region as well. Two cold cores are detected within the portion of the filament that is covered by the NOEMA observations.

The goals of this study are to investigate how the physical, kinematic, and chemical properties relate and vary from clump to disc scales, and how the relationships between different scales influence the evolution of forming stars. W3 IRS4 contains a number of stages in the evolutionary sequence of high-mass star formation (that is cold starless cores, a MYSO withsurrounding hot core, and an UCH II region) next to each other in one region, making it both an efficient andinteresting case study for how young massive stars affect, and are effected by, their surroundings. Achieving these goals required observations which have both high angular resolution and sensitivity to larger angular scales.

In a previous study of part of W3 including W3 IRS4, Wang et al. (2012b) were able to separate the principle structures in the W3 IRS4 region and identify some chemical differences with SMA observations at 1.4 mm. However, they only had short-spacing information from IRAM 30 m observations for ¹²CO and ¹³CO, so could notexplore the variation of temperature and abundance in the region, and their limited velocity and spatial resolution (1.2 km s⁻¹ and 2.2′′ × 1.7′′ in the lines, corresponding to ≳2000 AU scales, respectively) meant they were not really able to probe the ≲1000 AU scales wherediscs are expected around high-mass protostars (for example Krumholz et al. 2009; Peters et al. 2010a; Kuiper et al. 2011; Klassen et al. 2016; Meyer et al. 2017, 2018; Kölligan & Kuiper 2018).

The CORE NOEMA and IRAM 30 m observations presented in this paper have a factor of ~11 smaller beam and ~2.5 times higherspectral resolution, while also including short-spacing information for the whole frequency range covered by the interferometric observations. Furthermore, the point-source sensitivity in our smaller beam is more than a factor of 10 higher (see Sect. 2 for details), such that we detect more lines at the hot core position. Our new observations therefore enable the kinematic, physical, and chemical properties of this complex region to be constrained and comprehensively explored from clump to disc scales.

The structure of this paper is as follows. In Sect. 2 we discuss our NOEMA and IRAM 30 m observationsand a summary of the data reduction. Results and analyses of the physical and kinematic properties of the different structures within the W3 IRS4 region are presented in Sect. 3, while analysis of the chemical properties and variation is presented in Sect. 4. These individual results are then synthesised in trying to understand the past, present, and future of the region, as well as how it compares to other sites of high-mass star formation, in Sect. 5. We then summarise our conclusions in Sect. 6.

Fig. 1 Overview of the W3 star formation region showing the larger-scale context of W3 IRS4. Left: infrared three-colour image of W, where red is Herschel Hi-GAL (Molinari et al. 2016) 250 μm, green is Hi-GAL 70 μm, and blue is MSX 21 μm (Price et al. 2001). The three brightest far-IR sources, W3 OH/H2O, W3 IRS4, and W3 IRS5, are indicated with black circles. The region for the middle panel is indicated with a black dashed box.Middle: three-colour image of W3 IRS4 and IRS5, where red is 14.9 GHz radio free-free continuum from Tieftrunk et al. (1997) to emphasise ionised gas, green is Spitzer IRAC 4.5 μm (Ruch et al. 2007), and blue is Large Binocular Telescope (LBT, Hill et al. 2006) Ks data from Bik et al. (2012). The locations of W3 IRS4 and IRS5 are indicated with black circles, with the white dashed box indicating the region covered in the left-hand panel. Right: three-colour image focusing on W3 IRS4, where red is the NOEMA 1.4 mm continuum presented in this paper, imaged with natural weighting, while green and red are as in the middle panel. The dashed white circle shows the full-width half-maximum of the primary beam for the NOEMA continuum observations.

Fig. 2 NOEMA 1.4 mm continuum intensity maps of the W3 IRS4 region, imaged with the Clark algorithm (see text for more details) and with natural (left and middle) and uniform (right) weighting. The beam sizes are 0.94′′ × 0.82′′ and 0.45′′ × 0.32′′ respectively, shown in the lower-left of each panel. The black contours indicate intensity iso-contours at 5, 10, 15, 20, 25 and 30σrms (σrms = 0.95 and 0.27 mJy beam−1 for the left and right panels, respectively), while white-dashed contours indicate −5σrms. The left and right panels are shown in flux density units while the middle panel has been converted to brightness temperature using the Rayleigh–Jeans approximation. The principle structures within the region are indicated in the middle panel: the hot core/MYSO W3 IRS4 with a solid black circle, the UCH II region W3 C with a dashed circle, two cold continuum sources with black crosses, and the filament with white lines. The coordinates are given in Table 2. Short-spacing information is not available for the continuum emission.

2 Observations & data processing

This section describes the setup, observations and processing of the NOEMA (Sect. 2.1) and IRAM 30 m (Sect. 2.2) data used in this paper, as well as a summary of the final data products. A more detailed discussion of the data processing and results for the 30 m data is presented in Appendix A, while Appendix B discusses the details of the combination of the interferometric and single-dish data and imaging.

2.1 NOEMA

Observations were carried out with the Northern Extended Millimeter Array (NOEMA, formerly the Plateau de Bure interferometer), located in the French Alps, for all CORE sources at 1.4 mm in track-shared pairs with three array configurations so that the resulting maps are sensitive to a range of spatial scales. Details of the CORE sample, including selection criteria, are given in Beuther et al. (2018). W3 IRS4 was paired with the nearby source W3 OH/H₂O, for which 5 tracks passed quality control checks (see Table 3 for details).

The data were taken using two correlators working in parallel. The broad band correlator, WideX, covered a total bandwidth of 3.6 GHz centred around 219 GHz (see Table 4) for up to eight antennas in two linear polarisations (H and V) with a spectral resolution of ~2 MHz (corresponding to ~2.7 km s⁻¹ at 1.4 mm). The Narrow-Band correlator was used with eight spectral units providing a bandwidth of 80 MHz each with a spectral resolution of 0.312 MHz (corresponding to 0.42 kms at 1.4 mm), which were placed to cover a number of key transitions in H polarisation only (see Table 4 for details). This second correlator could only processthe signals of six antennas. During the observations, the rest-frame of W3 IRS4 was taken to be −42.8 km s⁻¹, based on theIRAM 30 m ¹²CO 2−1 observations of Wang et al. (2012b), but for most species in our observations the peak intensity at the location of IRS4 is at approximately −44.5 kms, so this latter value for thev_LSR is used in the remainder of the paper.

The data were calibrated using observations of bandpass calibrator 3C454.3 for the track in D-array and 3C84 for all other tracks. The same flux calibrator (MWC349), and the same phase and amplitude calibrator (0059+581), were used for all tracks. Calibration was carried out using the CLIC module within GILDAS². For line data, the continuum was subtracted using line-free channels in the same band, while a combined continuum dataset was formed from all WideX band channels which are completely free of line emission at all positions. The beam-size of the continuum with uniform weighting (which results in the highest resolution) is 0.45′′ × 0.32′′, increasing to 0.94′′ × 0.82′′ with naturalweighting. The identification of detected transitions towards W3 IRS4 is discussed in Sect. 4 with the properties of these lines given in Table C.2.

2.2 IRAM 30 m

In order to obtain short-spacing information for the molecular line observations, on-the-fly maps were observed for all sources using the EMIR multi-band mm-wave receiver (Carter et al. 2012) on the IRAM 30 m telescope at Pico Veleta over 5 nights: 13th–16th March and on 16th May 2015. The receiver was tuned so that the four 4 GHz windows, two in the lower side-band and two in the upper side-band, were centred at approximately 215, 219, 231, and 234 GHz respectively, with a small overlap between the two windows in the same side-band, such that the 219 GHz window covers the spectral setup of the NOEMA observations. The Fourier Transform Spectrometer (FTS) backend uses 195 kHz channels, which correspond to 0.27 km s⁻¹ at 217 GHz and 0.25 km s⁻¹ at 233 GHz.

All maps cover a total extent of 1.5′ by 1.5′ centred on the source, with full-width at half-maximum (FWHM) beam-sizes of 11.3′′ at 217 GHz and 10.6′′ at 233 GHz. In order to reduce potential scanning artefacts, a basket-weave observing strategy was employed such that one complete dataset includes a pair of maps, one with the major scanning direction in right ascension and one in declination. The spatial sampling was better than Nyquist, approximately every 4′′, with each position observed in both horizontal (H) and vertical (V) polarisation. The reference positions for each source were checked prior to starting the mapping observations to ensure that they were free from emission.

The data were reduced and processed using both standard and custom routines with the CLASS module within GILDAS. Details are given in Appendix A.1. The final maps have a typical noise of ~0.25 K in ~0.26 km s⁻¹ channels. The systematic calibration uncertainty for such data observed with the IRAM 30 m is typically of order 15% (see for example Druard et al. 2014). An overview of the 30 m data for W3 IRS4 is presented in Appendix A.2, including a sample spectrum at the central position over the full frequency range covered by the observations (Fig. A.1) and maps showing the integrated intensity, peak velocity, FWHM, and line shape for these transitions (Fig. A.2). The transitions with significant detected towards W3 IRS4 are detailed in Table C.1.

2.3 Merging and imaging

For those transitions with extended emission, the IRAM 30 m data are combined with the NOEMA visibilities to provide short-spacing information, which ensures full flux recovery and removes negative artefacts which can be induced in by spatial filtering. Short-spacing information is not available for the continuum. We compare the merged and NOEMA-only data for all lines and only use the merged data if the line is well-detected in the 30 m data and the merging results in a clear improvement in the quality of the final image. Where the merged data are used, through comparison of the original 30 m fluxes with the final merged maps, we estimate that >90% of the total flux, as measured in the single-dish observations, is recovered in the combined image.

The choice of imaging method can have a significant effect on the resulting maps. Two imaging techniques are used in this paper, Clark (Clark 1980) and Steer-Dewdney-Ito (SDI; Steer et al. 1984), depending on whether the distribution of emission is point-like or extended respectively. Transitions in the high- and low-resolution bands (see Table 4) are imaged with a velocity resolution of 0.5 and 3.0 km s⁻¹ respectively.The weighting given to visibilities on different baselines for a given dataset is chosen to balance spatial resolution against sensitivity (for example uniform or natural weighting, see Fig. 2). More details on the merging and imaging are given in Appendix B.

For W3 IRS4 the merged data are used with the SDI cleaning method for ¹³CO 2−1, C¹⁸O 2−1, SO 6₅ −5₄, H₂CO 3_0,3 −2_0,2, H₂CO 3_2,2 −2_2,1, HC₃N 24-23, and CH₃OH 4_2,2 −3_1,2, which all show extended emission, while the NOEMA only data and Clark method are used for OCS 18−17, HNCO 10_1,9−9_1,8, and CH₃CN 12−11 K = 0−6, which only show compact emission. The principle data products used in this paper will be made available on the CORE webpage³ upon publication.

3 Kinematics and physical conditions in the components of the W3IRS4 region

The emission towards W3 IRS4 is both geometrically and spectrally complex, including a number of different environments and stages in the high-mass star formation process, from cold condensations and a filamentary structure, through a mid-IR bright but radio-quiet MYSO with hot-core chemical emission, to a spatially extended UCH II region (see for example Fig. 2). The maximum brightness temperature in the NOEMA continuum image is ~4 K (middle panel of Fig. 2), so the continuum is optically thin everywhere. In the following section we start by focusing on the UCH II region in Sect. 3.1, then the cold filamentary structure in Sect. 3.2, followed by the hot core at the centre of the region in Sect. 3.3. The chemical properties of the region are discussed in Sect. 4.

3.1 The W3 C UCH II region

3.1.1 Comparison with radio continuum

The H II region which lies to the SW of the hot core IRS4, traditionally referred to as W3 C (Harris & Wynn-Williams 1976; Colley 1980; Tieftrunk et al. 1997), is resolved by our observations, showing a ring-like structure in the NOEMA continuum emission (see Fig. 2). Tieftrunk et al. (1997) mapped a large region of W3 Main with multiple array configurations of the VLA at 4.9, 14.9, and 22.5 GHz (6, 2, and 1.3 cm) with comparable resolution to our NOEMA data, identifying and analysing many of the H II regions in W3, including W3 C.

The spectral index between the three frequencies is consistent, within the uncertainties, with optically thin free-free emission for an H II region (that is α = −0.1). The free–free contribution in the NOEMA continuum data can therefore be estimated by re-scaling the VLA 22.5 GHz emission map, assuming that the intensity varies with frequency as ν^−0.1 (that is multiplying by (220/22.5)^−0.1 = 0.796), and convolving it to the same beam, grid as the NOEMA data. Figure 3 shows the NOEMA continuum emission before and after the estimate of the free–free contamination has bee subtracted, revealing that the ring of emission associated with W3 C is almost completely dominated by free-free emission from the ionised gas.

The radio data can also be used to calculate some of the properties of W3 C. For an optically thin H II region, the electron density (n_e) is given by (following Sánchez-Monge et al. 2013): (1)

where F_ν is the integrated flux density of the radio continuum emission at frequency ν, d is the distance to the source, θ_s is the angular diameter of the emission, and T_e is the electron temperature. W3 C has an integrated flux density of 0.45 Jy at 22.5 GHz, where the radio continuum emission has an angular diameter (θ_s) ~5.7′′, corresponding to 11,400 AU or 0.055 pc, and Tieftrunk et al. (1997) derive an electron temperature (T_e) of ~8000 K for W3 C from comparison between the radio continuum and H66α emission. The H II region therefore has an electron density of 1.4 × 10⁴ cm⁻³, considerably less than the ≳10⁶ cm⁻³ expected for a hyper-compact H II region (Kurtz 2005), so W3 C is considered an ultra-compact (UC)H II region despite its size placing it on the border between these two classes.

The ionising photon flux (Q₀) can also be estimated from the radio continuum flux (again following Sánchez-Monge et al. 2013) as: (2)

and is calculated to be 2 × 10⁴⁷ photons s⁻¹ for W3 C. There are no indications of an outflow from the powering source, which would suggest strong ongoing accretion, and even the most intense shocks are unlikely to produce a flux of more than a few hundred times the interstellar radiation field given the conditions in this and other star formation regions (see for example Neufeld & Dalgarno 1989; Melnick & Kaufman 2015; Karska et al. 2018; Hull et al. 2016), which given the size of the UCH II region would only generate of order 10⁴⁴ photons s⁻¹. Therefore this photon flux almost certainly comes from the photosphere of one or more high-mass stars.

Using the properties given in Table 1 of Davies et al. (2011) together with the relationship between spectral type and photospheric temperature from Martins et al. (2005) and Boehm-Vitense (1981) for zero-age main-sequence (ZAMS) O and B stars respectively,this ionising photon flux approximately corresponds to a single O9.5 zero-age main-sequence star, which has a mass of ~16 M_⊙ and a luminosity of ~2.2 × 10⁴ L_⊙. As the combined bolometric luminosity for W3 IRS4 and W3 C is ~4.5 × 10⁴ L_⊙ (see Sect. 1), this would leave ~2.3 × 10⁴ L_⊙ being contributed by IRS4 and other sources. If the powering source was in fact an equal-mass binary, then the observed ionising photons would come from two B0 ZAMS stars with individual masses of ~14 M_⊙ and luminosities of ~1.8 × 10⁴ L_⊙, leaving only ~9 × 10³ L_⊙ for other sources. Due to the steep decrease in ionising photons with decreasing stellar mass, a larger number of (individually slightly lower-mass) stars would further increase the luminosity contribution for the same ionising photon flux. However, given the relative brightness in the high-resolution near and mid-infrared of IRS4 with respect to W3 C (see Fig. 1), it seems unlikely that the UCH II region dominates the luminosity to such a large degree, so in the following we prefer the values obtained assuming a single ZAMS star powering the W3 C.

Fig. 3 Radio contribution to the NOEMA continuum emission around W3 C. Left: VLA 22.5 GHz radio continuum observations from Tieftrunk et al. (1997) rescaled to 220 GHz (see text for details). W3 C is in the centre, with weak emission from W3 IRS4, highlighted with a white dashed circle, to the NE and W3 D further away to the NW. Middle: NOEMA 220 GHz continuum emission. Right: NOEMA continuum map after the free–free contribution has been estimated from the VLA observations and removed (see text for details). All data were convolved to the same 0.8 × 0.8′′ beam. The black contours indicate intensity iso-contours at 5−30σ (σ = 0.4, 0.8 and 1.4 mJy beam−1 from left to right), increasing in 5σ increments, with the addition of a lower contour at 3σ for the VLA observations to show that the weak emission at the location of W3 IRS4 is significant. The white-dashed contours indicate −5σrms.

Fig. 4 Spatial and spectral distribution of gas and dust around the UCH II region W3 C. The colour maps show the integrated intensity in the combined NOEMA and 30 m data for 13CO 2−1 (top-left), C18O 2−1 (top-centre),c-C3H2 6−5 (middle-left), and H2CO 30,3 −20,2 (middle-centre), and the NOEMA-only continuum (top-right). The velocity integration ranges are given in the text. The grey contours indicate 5σ for the continuum map (σrms = 0.95 mJy beam−1). The beam is shown in the lower-left corner of each map. c-C3H2 and H2CO are shown using a square-root intensity scheme to emphasise the extended emission. The middle-right panel shows the fractional intensity of the continuum and gas emission as a function of radius from the centre of W3 C, averaged over θ. The centre was determined by hand and is indicated on the maps with a cyan +. The region containing emission related to the hot core W3 IRS4 and cold filament, indicated with a red polygon on the maps, is excluded from the radial analysis, while the maximum radius is indicated with a white dashed circle. The radius of peak emission for each tracer is indicated with a circle on the relevant maps. Gas spectra at three representative positions, indicated in the maps with red crosses which are numbered in the continuum panel, are shown in the lower panels.

3.1.2 Molecular line emission

Figure 4 shows how the 220 GHz continuum emission around W3 C compares to the integrated intensity and radial intensity distribution of four molecular gas tracers, ¹³CO 2−1, C¹⁸O 2−1, c-C₃H₂ 6−5, and H₂CO 3_0,3 −2_0,2, which are particularly associated with this UCH II region. The integration for the three species at low spectral-resolution (3 km s⁻¹) is from −81.8 to −0.8 km s⁻¹ for ¹³CO and C¹⁸O, and −48.8 to −28.8 km s⁻¹ for c-C₃H₂ in order to avoid contamination at the position of IRS4 from ³³SO. For H₂CO, which was observed at high spectral-resolution (0.5 km s⁻¹), −55.3 to −29.8 km s⁻¹ were used. The radial intensity distributions are calculated with offset from the approximate geometric centre of the NOEMA continuum emission (see Table 2), excluding the region encompassing the filament and W3 IRS4, indicated in Fig. 4 with a red polygon. Three example spectra at different positions are also shown, chosen to show the line profiles in the molecular gas on either side of the UCH II region, as well as in the extended, finger-like feature seen in c-C₃H₂.

These gas tracers, along with SO 6₅ −5₄ and possibly HC₃N 24−23 (see Sect. 4.1), show part or all of a ring-like structure around the UCH II region. The free–free continuum emission peaks around a radius of 1.6′′ (3200 AU), and drops off to a background level before the peak in any of the gas tracers. Of the molecular species, c-C₃H₂ peaks closest to the UCH II region at 3.6′′ (7200 AU), followed by ¹³CO, C¹⁸O, and H₂CO at 4.9, 5.3, and 6.1′′ (9800, 10 600, and 12 200 AU) respectively. Only ¹³CO shows much overlap with the free–free emission, though at lower intensity than outside it.

All lines are spectrally broad in the vicinity of W3 C (see Fig. 4), with the higher spectral-resolution H₂CO spectra showing signs of two velocity components, one around −44.5 km s⁻¹, the v_lsr of IRS4, and another around −42.5 km s⁻¹ which becomes stronger for more southern positions. Thse two velocity components could be associated with the front and the back of a shell of emission, or with different parts of the swept up gas along the line-of-sight if the geometry is more complicated than a simple spherical shell. Alternatively, the H66α radio recombination line emission from W3 C detected by Tieftrunk et al. (1997) shows a small region near the centre which is at −25 km s⁻¹, offset by ~20 km s⁻¹ from most of the ionised gas which at velocities around those we observe for the molecular gas (see their Fig. 3 and Sect. 3.2). We do not detect any gas emission at this velocity, but the redshifted ionised gas emission could be associated with the powering source or the far side of the UCH II region. Unfortunately, it is difficult to be more conclusive about the expansion velocity or depth of the shell along the line-of-sight from the NOEMA data because of the blending of the velocity components, which have a relatively small velocity offset with respect to their line-width, as well as the lack of higher spectral resolution, particularly for ¹³CO.

The presence of a ring of molecular gas outside the ionisation front, as traced by the free-free continuum, is to be expected from such an UCH II region, as a photo-dissociation region (PDR) would be formed by the ionising radiation pushing into the surrounding clump material (for example van der Wiel et al. 2009; Guzmán et al. 2011; Trevi no-Morales et al. 2014; Cuadrado et al. 2017). c-C₃H₂ peaks closerto the ionised gas than any other species, which is entirely consistent with its association with ultra-violet (UV) illuminated gas (see for example Pety et al. 2005; Fontani et al. 2012; Guzmán et al. 2015). ¹³CO and C¹⁸O then peak next, followed by H₂CO (and SO, though this is not shown). This is the opposite order to what is observed in the Orion Bar PDR (that is H₂CO slightly closer to the ionising source than C¹⁸O in that region, see for example van der Wiel et al. 2009; Cuadrado et al. 2017), but this could be due to differences in the radiation field, the temperature and density of the gas, or even the age of the region.

C¹⁸O has a lowerrelative intensity near the centre of the UCH II region than ¹³CO (see lower-right panel of Fig. 4). This could be due to the peak ¹³CO intensity being decreased at the edge of the shell due to higher optical depth with respect to the centre and C¹⁸O. However, it is also consistent with isotope-selective photodissociation by FUV photons, where ¹³CO is able to more effectively self-shield than C¹⁸O at smaller radii due to its higher abundance (see for example Glassgold et al. 1985; Röllig & Ossenkopf 2013; Shimajiri et al. 2014).

3.1.3 Powering source

One last question the NOEMA data allow us to (re-)consider is which source(s) is(are) powering W3 C. Bik et al. (2012) performed high-resolution near-IR imaging and multi-object spectroscopy of W3 Main using LUCI 1 (Ageorges et al. 2010; Seifert et al. 2010) on the Large Binocular Telescope (LBT, Hill et al. 2006). They suggested one particular source, indicated with a white circle in Fig. 5, as being the ionising source for W3 C due to its location within the radio emission and that its spectral type is consistent with the ionising flux of the H II region. However, this source is on the northern edge of the free–free emission rather than at the centre (see Fig. 5). Furthermore, in order to produce such a skewed champaign-flow like UCH II region there would need to be higher density material preventing expansion to the north and west, which is not consistent with the lack of strong molecular emission in dense gas tracers on those sides of the source (see Fig. 4). It would therefore be expected that the radio continuum emission would extend further north and west if this were truly the ionising source for W3 C.

There is another source in the LUCI K_s image near the centre of the radio emission, indicated in Fig. 5 with a red circle, which may be the powering source. This source was not part of the spectroscopic observations and is only detected in K_s (K_s = 15.86 mag., H >19 mag., Bik et al. 2014), suggesting considerable extinction consistent with being deeper within the UCH II region than the northern source. There is also a slight enhancement of emission in the continuum at 10 μm (Lumsden. priv. comm.), and a red-shifted H66α component (Tieftrunk et al. 1997) around this position. The higher extinction towards this source could come from dust inside the UCH II region which is too small and hot to be detected at mm wavelengths (see for example Everett & Churchwell 2010; Paladini et al. 2012; Anderson et al. 2012), and/or from dust in the front-part of the UCH II region which is not strong enough to be detected by NOEMA. Deep near-IR spectroscopy imaging is needed to constrain the extinction and spectral type of this alternative source, and thus ascertain whether it is consistent with the observed properties of W3 C.

Fig. 5 Comparison of near-IR Ks images from Bik et al. (2012, greyscale image, 0.9′′ resolution) with VLA 22.5 GHz radio continuum observations from Tieftrunk et al. (1997, orange contours) and free–free-subtracted NOEMA continuum at 220 GHz (magenta contours). The contour levels are the same as in Fig. 3. The cyan cross indicates the approximate geometrical centre of the UCH II region used as the centre for the radial analysis in Fig. 4, while the white circle indicates the source attributed by Bik et al. (2012) as powering W3 C and the red circle indicates the nearest source to the centre of the free–free emission, which is much redder.

Fig. 6 Maps of integrated intensity (moment 0, left, integrated between −55.3 and −29.8 km s−1), intensity-weighted velocity (moment 1, middle), and intensity-weighted FWHM (moment 2, right) for combined NOEMA and 30 m data for H2CO 30,3 −20,2, focusing on the filament. The grey contours indicate 5σ for the NOEMA 1.4 mm continuum with uniform weighting, as in the right-hand panel of Fig. 2. The cyan contours show the moment 0 of the H2CO 32,2 −22,1 transition, starting at 5σ (σ = 0.019 Jy beam−1 km s−1) and increasing to 25σ in 5σ steps. The mean velocity near the hot core is approximately −44.5 km s−1. The moment1 and 2 maps are calculated using channels above 4σ and only for pixels where the integrated intensity is also above 4σmom0, with grey used to indicate those regions that do not meet these requirements.

3.2 Filament

On the other side of W3 IRS4 from W3 C is a filament-like continuum structure that extends to the north-east of the hot core and includes two colder condensations, which show no signs of star formation and so are in an earlier evolutionary stage. Within the NOEMA continuum emission, this has a length of ~8.8′′ (17600 AU, 0.09 pc), not including IRS4 itself, and a width of ~1–4′′ (2000-8000 AU, 0.01–0.04 pc) on the plane of the sky, though as mentioned in Sect. 1 this structure is seen to extend much further, up to ≳40′′ (80 000 AU, 0.4 pc) in VLA NH₃ observations (Tieftrunk et al. 1998).

This structure is also detected in a number of the molecular gas tracers, including H₂CO, for which two transitions are included in the high spectral-resolution NOEMA bands (3_0,3−2_0,2 and 3_2,2 −2_2,1). Figure 6 shows moment maps of the combined NOEMA and 30 m data for the lower-energy transition, H₂CO 3_0,3 −2_0,2, with contours of the moment 0 of the 3_2,2 −2_2,1 transition overlaid on the first panel. The data were imaged with an intermediate weighting to balance sensitivity against spatial resolution, and primary beam correction has been performed on the moment 0 map. The two transitions are well correlated, with the lower energy transition being more readily detected over a larger region.

One of the most noticeable features of these maps is the presence of a large-scale shift in peak velocity of several km s⁻¹ which startsin the north-east and ends near the position of IRS4 near the centre of the region. The material associated with W3 C has a similar velocity to that of W3 IRS4. Indeed, this large-scale velocity shift is even seen the 30 m data for H₂CO, SO, C¹⁸O, and, to a lesser extent, ¹³CO (see Fig. A.2). The combined NOEMA and 30 m data reveal that this velocity shift is relatively smooth and associated with the filament-like structure seen in the continuum, with the gas emission being more extended than the dust by ~1′′ (2000 AU, 0.01 pc) on all sides.

3.2.1 XCLASS fitting of H₂CO

The two H₂CO transitions included in the high spectral-resolution bands can be used to trace the gas temperature in the filament. To do this, the data for the two transitions were convolved to the same spatial resolution (0.8 × 0.8′′) and fitted simultaneously using XCLASS (Möller et al. 2017). More details on the specific setup and assumptions used are given in Appendix C.1. This returns the excitation temperature, column-density of p−H₂CO, peak velocity, and FWHM for each pixel, the maps of which are shown in Fig. 7. Assuming that H₂CO is in LTE (this is discussed further in Sect. 3.2.3), the excitation temperature is equal to the kinetic temperature of the gas.

The values for the fit to the continuum peak, as shown in the spectra of Fig. 7, as well as for the two cold cores, are given in Table 5. The median values over the filament are T_ex = 50 K, N[p-H₂CO] = 5.8 × 10¹³ cm⁻², = −45.3 km s⁻¹, and FWHM = 2.9 km s⁻¹. There are hints that the lines do not consist of a single velocity component around the continuum peak, but attempts to fit multiple velocity components to the line produced worse fits than using a single velocity component. Elsewhere, the lines show only a single Gaussian component.

It should be noted that the upper-level energy of the higher energy transition, 3_2,2 −2_2,1, is 68 K, (E_up = 21 K for the 3₀₃ −2₀₂ transition) so the excitation temperature becomes increasingly uncertain above ~100 K, particularly for high column densities. In such regions the excitation temperature and column density are certainly high, but solutions with a large range in T_ex result in very similar line intensities, and thus similar fit probabilities. Fortunately, these high excitation temperatures are limited to the region around W3 IRS4 where CH₃CN is an effective tracer, which is discussed further in Sect. 3.3. In regions of the map where the upper-energy transition is not well detected (≤3σ_mom0, white dashed contours in Fig. 7), the temperature is less well constrained and probably an overestimate.

The overall temperature structure of the region shows a clear peak near the position of IRS4. Though the peak temperature is offset slightly to the SE of the continuum peak, it does coincide with the peak in intensity and column density of H₂CO. The excitation temperature decreases moving from the position of IRS4 to the NE, first steeply to the median temperature of 50 K, then becoming cooler near the position of the second cold core. The temperatureof most of the gas within the filament is therefore well above that for CO freeze-out (~20−40 K, for example Fayolle et al. 2016; Anderl et al. 2016), but below that for the sublimation of H₂ O (~100 K, for example Fraser et al. 2001). This affects the chemistry in both the gas and solid phases, by both increasing reaction rates and changing the composition of the gas and ice, compared to the typical conditions considered for dense molecular gas (10−15 K).

Fig. 7 XCLASS fitting for H2CO in the filament. The panels show: integrated intensity (between −55.3 and −29.8 km s−1) in the H2CO 30,3 −20,2 (top-left) and 32,2−22,1 (bottom-left) transitions; example spectra (black) and  XCLASS model fits (red) for H2CO 30,3 −20,2 (top-middle-left) and 32,2−22,1 (bottom-middle-left); column density in p-H2CO (top-middle-right); excitation temperature (top-right); peak velocity (bottom-middle-right); FWHM (bottom-right). The data for the two lines were resampled to a common spherical beam with a FWHM of 0.8′′, as shown in the top-right panel, and converted to the TMB scale. The grey contours in the maps show the continuum smoothed to the same resolution as the line data, starting at 5σ (σ = 0.81 mJy beam−1) and increasing to 45σ in 10σ steps. The spectra are extracted at the continuum peak (coordinates in Table 2), which is indicated in the left-hand panels with a red marker. A scale-bar is shown in the lower-left panel. The model results are only displayed for pixels where the integrated intensity for H2CO 30,3 −20,2 is above 5σmom0, with grey used to indicate those regions that do not meet this requirement. The white dashed contours in the map of excitation temperature indicate 3σmom0 for H2CO 32,2 −22,1.

3.2.2 Mass

If we assume that the dust is optically thin and in thermal equilibrium with the gas, then the gas temperature derived from H₂CO and the dust continuum emission can be used to calculate maps of N[H₂] and mass following the prescription by Hildebrand (1983) in the forms presented by Schuller et al. (2009): (3)

and (4)

Here, S_ν and F_ν are the peak and integrated flux density, R is the gas-to-dust mass ratio, B_ν(T_dust) is the Planck function at dust temperature T_dust, κ_ν is the dust absorption coefficient, θ_B is the beam solid angle, μ is the mean molecular weight, m_H is the mass of a hydrogen atom, and d is the distance to the source. Following Beuther et al. (2018) and Ahmadi et al. (2018), we use R = 150 (Draine 2011), μ = 2.8 (Kauffmann et al. 2008), and κ_ν = 0.9 cm² g⁻¹ from Ossenkopf & Henning (1994) for a wavelength of 1.4 mm assuming coagulation for 10⁵ yr at n = 10⁶ cm⁻³ with thin ice mantles, starting from the interstellar grain distribution from Mathis et al. (1977, typically referred to as the MRN distribution).

Maps of mass-per-pixel and H₂ column density were calculated using the temperature map derived from H₂CO (see Fig. 8). The associated 1σ statistical uncertainties, calculated by propagating the noise in the continuum map to the mass-per-pixel and column density are 3 × 10⁻⁴ M_⊙ pixel⁻¹ and 2.6 × 10²² cm⁻² respectively. The region within a radius of 1′′ of the hot core contains a mass of 1.1 M_⊙ of gas and dust, while the rest of the filamentary structure has a total mass of 19.5 M_⊙, including the two colder continuum peaks (see Table 6 for their individual properties). The median column density is 2.8 × 10²³ cm⁻² with a maximum value of 1.8 × 10²⁴ cm⁻². These valuesfor the masses and H₂ column densities are likely underestimated due to the large amount of missing continuum flux towards W3 IRS4 (~78%, Beuther et al. 2018).

Fig. 8 Calculated physical properties for the filament. The panels show: mass-per-pixel (left); H2 column density (middle); and volume number density (right) assuming a depth along the line-of-sight equal to the widest part of the filament (~8000 AU). The volumedensity is therefore likely an underestimate. The values are only calculated where the continuum is above 3σ (σ = 0.81 mJy beam−1). The blackcontours show the continuum smoothed to the same resolution as the line data, starting at 5σ and increasing to 45σ in 10σ steps. The beam and a scale-bar are shown in the left and middle panels respectively.

3.2.3 Volume density

The third panel in Fig. 8 shows an estimate for the volume number density, n, calculated from N[H₂] by conservatively assuming that the depth of the emission along the line-of-sight is equal to the largest width on the plane of the sky (~4′′ or 8000 AU), perpendicular to the long axis of the filament. The depth along the line-of-sight is likely a lower limit to the true value, particularly around the hot core, since the filament is narrower for much of its length by up to a factor of ~4 and N[H₂] is likely underestimated due to the missing flux in the continuum map. The median value of 2 × 10⁶ cm⁻³ is well above the effective critical density for H₂CO under the conditions in the region (n_eff ~10⁵ cm⁻³ for 10 K, ~10⁴ cm⁻³ for 100 K, Evans 1999; Gerner et al. 2014; Shirley 2015), so the assumption of LTE for the XCLASS modelling, and therefore that the gasand dust are in thermal equilibrium, is sound.

3.2.4 Stability

For regions within the filament, the thermal Jeans mass can be measured using: (5)

where γ is the adiabatic index of the gas, for which we use 5/3 since, although the gas is mostly diatomic there is unlikely to be enough energy to excite ro-vibrational states. Table 6 gives the results for W3 IRS4 and the two cold cores, which reveal that the material around the hot core is stable to further fragmentation on core scales, but the two cold cores may fragment if there are no other forces supporting them.

The stability of the filament as a whole against further fragmentation can be explored through comparison with the critical mass-per-unit-length for an infinite, isolated cylinder to be in equilibrium between thermal pressure and gravity (see for example Ostriker 1964; Kainulainen et al. 2016). This is given by: (6)

Using the median temperature of 50 K, the filament has a value of M_l,crit ~114 M_⊙ pc⁻¹. If the filament is entirely in the plane of the sky then the structure we observe has a length of ~17 600 AU and thus a mass-per-unit-length M_l of 229 M_⊙ pc⁻¹. The observed filament must be at an angle >60° with respect to the plane of the sky before the observed mass-per-unit-length is below the critical value.

However, the median FWHM of H₂CO is 2.9 km s⁻¹, while the thermal contribution to the FWHM at the median temperature of 50 K is 0.28 km s⁻¹, with similar ratios even in the cold cores, so the observed line-widths are strongly dominated by non-thermal motions. Given the large non-thermal contribution to the line-width it is prudent to consider the effect that this contribution might have on the support of the cores and filament (see for example Fiege & Pudritz 2000; Kirk et al. 2015). If the thermal sound-speed is replaced with the observed FWHM/, thus assuming that all the non-thermal contribution to the line-width provides turbulent support against gravity, then Eq. (5) becomes: (7)

and Eq. (6) becomes: (8)

Using the median FWHM for the filament from H₂CO of 2.8 km s⁻¹, M_{l,crit,nonthermal} = 658 M_⊙ pc⁻¹, so the filament could be stable against radial collapse if supported by non-thermal motions. Similarly, including the non-thermal motions shows the region around W3 IRS4 to be very stable. The second cold core, which is closer to IRS4, has a non-thermal Jeans mass very close to the observed mass (M = 2.4 M_⊙ and M_J,nonthermal = 2.2 M_⊙, see Table 6), while the mass of the first cold core is still larger than the non-thermal Jeans mass (M = 5.4 M_⊙ and M_J,nonthermal = 3.2 M_⊙).

There are no indications of outflows or infall in our data associated with either cold core, and little mid or near-infrared emission associated with them either, suggesting that star formation has not yet begun in these cold cores. Any collapse takes place on one or more free-fall timescales, given by: (9)

which is of the order of 10⁴ yr for the three regions considered (see Table 6).

3.2.5 Flow along the filament

The stability analysis discussed above does not, however, take into account the effect of the large-scale velocity gradient (see Figs. 6 and 7). Given the moderate line-widths, the lack of distinct high-velocity non-Gaussian line wings, and the relatively high densities (see Fig. 8) within the filament, the observed velocity gradient in the filament is probably not due to an outflow. The velocity difference between the top end of the filament and IRS4 is ~2.3 km s⁻¹, which if caused by a difference in kinematic distance would lead to a separation along the line-of-sight of 230 pc over less than 0.1 pc projected on the sky. This seems highly unlikely.

Were the observed velocity gradient due to solid-body rotation, this would imply a rotational period of order 10⁵ yr which is also unlikely on such size-scales, particularly given that we only detect part of a structure which extends for ≳40′′ in length (80 000 AU or 0.4 pc on the plane of the sky) to the north-east and north in the VLA NH₃ observationsof Tieftrunk et al. (1998). Those authors do not present velocity maps of the filament, but do mention in the text that they see lower line-widths and more negative peak velocities in the northern part of the filament compared to the part near IRS4, which is consistent with the velocity gradient continuing beyond the region we observe.

The filement does lie near the edge of a larger, more evolved H II region, W3 H (see Fig. 1 of Tieftrunk et al. 1998), and this may have played a role in its formation or compression (Wang et al. 2012b). However, this older H II region is unlikely to be the cause of the velocity gradient, since in that case the velocity gradient should be across the width of the filament in the direction of expansion of the bubble, rather than along the long axis around the bubble edge.

The remaining option is that the velocity gradient is due to a large-scale flow either towards or away from IRS4. In this case, the gradient of ~25 km s⁻¹ pc⁻¹ is the line-of-sight component of the velocity along the filament, while the observed distance is the plane-of-the-sky distance perpendicular to the line-of-sight. For a given angle, θ, between the axis of the filament and the plane of the sky, the mass flow rate is given by (10)

Treating the filament as a solid body results in mass flow rates along the axis of the filament of 3.3 × 10⁻⁴ M_⊙ yr⁻¹ and 1.0 × 10⁻³ M_⊙ yr⁻¹ for θ = 30° and 60° respectively, while performing the calculation pixel-by-pixel and summing over the map gives values of 2.5 and 7.6 × 10⁻⁴ M_⊙ yr⁻¹ respectively.

The direction of the flow along the filament, either towards or away from IRS4, remains uncertain. If the flow was away from IRS4, one might expect faster velocities in the lower-density material either side of the filament compared to the higher-density material traced by the continuum emission. However, no distinct changes are seen in the velocity field (see Fig. 7). Indeed, the velocity change along the filament is rather smooth, with no sign of internal jumps that might indicate shocks.

Alternatively, the flow could be moving material towards IRS4, possibly due to gravity. The emission associated with IRS4 and W3 C is at similar velocities, suggesting that there is no relative motion between them along the line of sight, and that the flow from the filament stops in the vicinity of IRS4. If the flow is towards IRS4 and the material close to the MYSO is in hydrostatic equilibrium, then one might expect heating or a sharper increase in temperature due to a shock where the flow meets the envelope of the hot core. Indeed, the peak in the H₂CO intensity and temperature is offset to the east of the continuum peak for IRS4 (see Fig. 7), although this offset could also be caused by optical depth effects and the insensitivity of the lines to higher temperatures, or by chemical and/or abundance variations. Regardless of the direction of the flow, the rate at which material is moving has a significant effect on the evolution of the region. However, for the remainder of the paper we assume that the flow along the filament is towards IRS4.

There are no sharp density or velocity changes at the edges of the cold cores which might indicate that they are de-coupled from the overall flow along the filament, so they appear to be moving with the flow. In this case, it would take only 2−6 × 10⁴ yr for the most massive cold core to travel the distance between it and IRS4 (~10 500 AU), comparable to the current free-fall timescale in the core. Thus, the conditions for any additional protostars which may be forming in the cold material are likely shaped by this motion.

3.2.6 Comparison with single-dish-only results

It is helpfulto quantify how these results from resolved observations compare to what one would obtain with only single-dish observations. A fit to the IRAM 30 m data for H₂CO (see Fig. C.1) gives T_ex = 67.5 K, N[p-H₂CO] = 9.6 × 10¹³ cm⁻², = −44.2 km s⁻¹, and FWHM = 4.2 km s⁻¹ for the central position. The values of T_ex, N[p-H₂CO]. and FWHM from the single-dish data alone are therefore higher than those for most of the filament but below that for the hot core. For the line-width, the overestimate by a factor of ~1.5 compared to the median in the filament from the map (2.9 km s⁻¹) is due to the presence of unresolved systematic motions, suggesting that single-dish surveys probably overestimate turbulent support on clump to core scales.

From the single-dish SCUBA flux of 14.99 Jy, the total mass in the region (not accounting for any free–free contamination from W3 C for simplicity) is 77.3 M_⊙ if the single-dish derived temperature (67.5 K) is used or 107.4 M_⊙ if the median temperature of the filament (50 K) is used instead. The beam-averaged H₂ column density similarly increases from 6.2 × 10²² to 8.6 × 10²² cm⁻² when using the lower temperature. Since most of the mass in the NOEMA map comes from regions at lower temperatures, the higher temperature in the single-dish results would lead to the mass of the region being underestimated and molecular abundances overestimated by a factor of ~1.4 compared to the resolved observations. However, since so much of the single-dish continuum flux is filtered out by the NOEMA observations and thus from extended clump gas which is probably warmer than the colder filament, the fraction of missing mass is probably lower than the fraction of missing continuum flux. Ultimately, though, single-dish only observations of similar high-mass star formation regions likely overestimate the temperature, meaning they underestimate the mass and H₂ column density.

3.3 The MYSO and hot core W3 IRS4

3.3.1 Radio emission and outflow direction

The MYSO and hot molecular core W3 IRS4 is the brightest structure in both continuum emission and many gas tracers in the NOEMA observations (see for example Figs. 2 and 4), as might be expected for one of the brightest infrared andsub-mm sources in W3. The combined bolometric luminosity for the W3 IRS4 region from single-dish measurements is ~4.5 × 10⁴ L_⊙ (see Sect. 1) and a contribution of ~2.2 × 10⁴ L_⊙ was estimated to this from W3 C based on its radio flux (Sect. 3.1). IRS4 dominates the infrared emission of the region (see Fig. 1), so we assign the remaining bolometric luminosity, similar to that from W3 C, to the MYSO. If the central source is a single star, then this suggests it is another 16 M_⊙, O9.5 star, while an equal-mass binary of this luminosity would be made of two 12 M_⊙, B0.5 stars.

The radio emission associated with IRS4 is unresolved and much weaker than that of W3 C (see Figs. 3and 5, as well as Tieftrunk et al. 1997), with an integrated flux at 22.5 GHz of only 1.7 mJy (~0.4% that of W3 C). Previous authors have assumed that the radio emission from IRS4 is optically thick, and thus that it is from a hyper-compact (HC)H II region (Tieftrunk et al. 1997; Wang et al. 2012b). However, HCH II and UCH II regions have similar infrared to radio luminosity ratios (cf. Fig. 6 of Hoare & Franco 2007). An alternative source for this radio emission, and one which typically results in a much fainter radio to bolometric luminosity ratio, is from a wind or radio jet (for example Martí et al. 1993; Hoare et al. 1994; Hoare & Franco 2007; Purser et al. 2016). Comparing to Fig. 6 of Hoare & Franco (2007), which compares L_bol and L_radio for known HCH II and UCH II regions, as well as wind and jet sources, W3 C lies within the group of H II regions while IRS4 clearly lies within the group of wind and/or jet sources. Thus it is much more likely that the hot core W3 IRS4 is a radio-quiet MYSO with the cm emission originating from a radio jet, rather than a source which powers an HII region. We discuss the implications this has for the central source of IRS4 in Sect. 5.

With the radio emission indicating the presence of an ionised radio jet or disc-wind, it is obvious to next try to constrain the direction and properties of the outflow from IRS4. However, neither the IRAM 30 m, NOEMA nor merged data show any noticeable sign of high-velocity line-wing emission, with similar line-shapes in the 30 m spectra between ¹²CO, ¹³CO and C¹⁸O (see for example Fig. A.2). On the other hand, infrared K_s and 4.5 μm continuum images show some emission and structures extending from IRS4 to the N and NW, and the peak in the 4.5 μm continuum band, which includes lines of shocked H₂ often associated with shocks in outflows (see for example Davis et al. 2007; Cyganowski et al. 2008; Ybarra & Lada 2009),is shifted slightly in this direction (see Fig. 9). Imaging the line-wings of ¹³CO in the merged data shows some hints of emission in this NW-SE direction in both wings, though they only extend to approximately ±13 km s⁻¹ from the v_LSR of IRS4, which is approximately −44.5 km s⁻¹ from H₂CO (see Table 5), CH₃CN and other dense-gas tracers. The only other line to show much emission at these velocities is SO, which can trace shocks and outflows (see for example Chernin et al. 1994), where the detections are compact and centred on IRS4 for both line-wings (see middle-panels of Fig. 9).

Based on this evidence, we tentatively suggest that the outflow from W3 IRS4 may be in a NW-SE direction (see purple arrows in Fig. 9), with the mid and near-IR continuum tracing the edges of thecavity (see dashed purple lines in Fig. 9, which have an opening angle of 40°, determined by eye). Given the lack of emission at larger velocity offsets from v_LSR the outflow is most likely close to the plane of the sky, in which case the line-wing emission in ¹³CO primarily traces transverse motions perpendicular to the direction of the outflow. The region of higher emission in the map of red-shifted ¹³CO to the SE of IRS4 is probably not due to the outflow, since the line-centre of the cloud shifts to slightly lower velocities in that region, but the emission to the NW could be from a bow-shock, as it lies along the suggested outflow axis, or the cavity walls. SO could be tracing the base of the jet, though it could also be due to accretion shocks onto the disc or interactions between a wind, jet or outflow and the dense envelope emission near the MYSO (cf. the observations of HD100546 by Booth et al. 2018).

3.3.2 Velocity structure of the hot core

With the existence and direction of the jet and outflow tentatively established, the structure and kinematics of material close to IRS4 can now be investigated, starting with a consideration of the dense gas and temperature tracer CH₃CN, which is often used to trace discs in high-mass star formation regions (for example Johnston et al. 2015; Cesaroni et al. 2017; Ahmadi et al. 2018). Figure 10 shows the results of XCLASS fitting to CH₃CN and CH₃¹³CN for the material around the MYSO and hot core. HNCO is also fitted simultaneously to correctly account for any overlap between the 10_1,9−9_1,8 transition of HNCO and CH₃CN 12₆ −11₆. Since CH₃¹³CN may be present but is not clearly detected, we assume a ratio of CH₃¹³CN to CH₃CN of 76, as measured by Henkel et al. (1982) for W3 and consistent with the calculations of Qin et al. (2015) for W3(H₂O). Following Ahmadi et al. (2018), we only fit the K = 4–6 transitions of CH₃CN as this givesmore consistent results.

Unlike H₂CO, the peak CH₃CN intensity coincides with the continuum peak. Temperatures of order 140 K are found at the continuum and intensity peak with little variation across the region, and a clear velocity gradient is seen in approximately the north-south direction. This is close to the axis for the outflow, but some curvature is seen in the transition region where the gas is close to v_lsr and the highest velocity gas in both the red and the blue side lies to the north-east of the proposed outflow axis. This hints at a more complicated situation that requires further investigation.

Figure 11 therefore presents the 1st and 2nd moment maps and spectra at the continuum peak of all species detected with sufficient signal-to-noise in the high spectral-resolution bands of our observations, shown in order of increasing upper-level energy (see Fig. 11 and Table C.2). W3 IRS4 does not show a single consistent velocity pattern between different molecules and transitions, in contrast to sources such as the nearby W3 H₂O (see for example Ahmadi et al. 2018). Towards IRS4, some transitions show complicated or unclear velocity structures (for example CH₃OH, HC₃N, CH₃CN 12₂ −11₂), while in others there is a shift in the axis of the velocity gradient from more NW-SE at lower energies (for example H₂CO) to N-S (for example HNCO, CH₃CN 12₃ −11₃, and CH₃CN 12₆ −11₆). The lines with N-S velocity gradients also tend to have higher line-widths compared to those with a NW-SE or unclear velocity structure.

Figure 12 shows three position-velocity (PV) diagrams, one along the outflow axis, one rotated by 45° from the outflow axis, and one perpendicular to it, for five of these transitions: H₂CO 3_2,2 −2_2,1, OCS 18-17, HNCO 10_1,9−9_1,8, CH₃CN 12₃ −11₃, and CH₃CN 12₆ −11₆. All the outflow-axis cuts (middle-left column in Fig. 12) show an additional intensity peak which is offset to the red (around −42 km s⁻¹) and the SE (around −0.5′′) of the continuum peak, though it is stronger and more extended in some lines than others. This feature is fainter but also seen in some of the cuts at 45° from the outflow axis (middle-right column in Fig. 12). Such a linear feature with velocity offset increasing with spatial offset can be caused by either solid-body rotation or a Hubble-like outflow. The latter would seem like the more likely option given the evidence presented from infrared and large-scale ¹³CO emission in Sect. 3.3.1, particularly since these species have all been detected in outflows from low-mass (for example Bachiller & Pérez Gutiérrez 1997; Codella et al. 2009; Rodríguez-Fernández et al. 2010) and high-mass YSOs (for example Palau et al. 2017).

If this is indeed the case, then the cut perpendicular to the outflow (right-hand column in Fig. 12) is therefore where one might expect to find any signature of a disc. Neither the line nor the continuum emission is strong enough to indicate that a disc is being hidden by optical depth effects (the continuum has a brightness temperature of ~4 K). However, these PV-diagrams do not show a clear signature of Keplerian rotation, with emission in all four quadrants of the map possibly suggesting a combination of infall and rotation (see for example Ohashi et al. 1997; Tobin et al. 2012; Oya et al. 2016). The distribution of the regions of red and blue-shifted material in the moment 1 maps of some species (most notably OCS, HNCO, and CH₃CN 12₃ −11₃, see Fig. 11), as well as the different shapes of the transition regions around v_lsr and the differences between different species, also suggest that there are multiple contributions to the velocity distribution across some or all of the region around IRS4.

Rotation or warping of the axis of the red-blue velocity shift in moment 1 maps has been seen in other sources. For example, Wang et al. (2012a) saw a warping of the red-blue axis in their SMA data of the MYSO AFGL2591-VLA3 which is similar to what we observe in the OCS map. They showed through radiative transfer models that such behaviour is consistent with a combinationof Keplerian rotation and Hubble-law type expansion from an outflow. However, rotation as a function of radius can also be seen where a rotating disc is embedded in an infalling envelope, as shown by van ’t Hoff et al. (2018), for example, in their modelling of ALMA data for the Class 0 low-mass YSO L1527 (see particularly their Figs. 5 and D.1). This shift is quite abrupt in their models, close to the outer radius of the disc, which is again reminiscent of the OCS and CH₃CN data for W3 IRS4.

The rotation in the direction of the velocity gradient in the moment 1 maps of the denser and higher-temperature tracing lines suggests that what we are seeing is a combination of outflow and infall on larger scales (see Sects. 3.3.1 and 3.2), with rotation on small scales. Since the radius at which the shift is most prominent in OCS is on the size of the ~0.3′′ beam, any possible disc around W3 IRS4 is most likely unresolved in our data and thus has as size ≲600 AU. The characterisation of the disc around W3 IRS4 therefore requires higher angular resolution observations in molecules which are not enhanced in outflows and trace denser material, possibly with the aid of radiative transfer modelling.

Fig. 9 LBT LUCI Ks (top-left, data from Bik et al. 2012) and Spitzer IRAC 4.5 μm (bottom-left) continuum maps, as well as integrated emission in the blue (top) and red (bottom) line-wings of combined NOEMA and 30 m maps of 13CO 2−1 (middle) and SO 6−5 (right). The NOEMA continuum emission with natural weighting is shown as green (left) or grey (middle and right) contours from 5 to 30σ in 5σ intervals (σ = 0.95 mJy beam−1). The velocity integration windows for the 13CO and SO maps are indicated in the right-hand panels. The purple arrows indicate the tentatively identified direction of the outflow, with dashed lines used to highlight features which may trace the outflow cavity walls (see text for details).

Fig. 10 XCLASS fitting for CH3CN in the hot core. The panels show: integrated intensity in the CH3CN 124 −114 transition (top-left); column density of CH3CN (top-middle); excitation temperature (top-right); peak velocity (bottom-middle); FWHM (bottom-right); example spectra (black) and  XCLASS model fits (bottom-left) for CH3 CN (red), CHCN (cyan), and HNCO (blue). As k = 0−3 of CH3 CN were notincluded in the fit, the resulting fit solution for those transitions is shown with a dashed line. The data were imaged with an intermediate weighting between uniform and natural in order to provide a good balance between sensitivity and resolution, resulting in a beam size of 0.67′′ × 0.56′′ (PA = −83.2°), as shown in the top-left panel with a scale-bar, and converted to the Tmb scale. The black contours in the maps show the continuum at the same resolution as the line data, starting at 5σ and increasing to 45σ in 10σ steps (σ = 0.48 mJy beam−1). The spectra are extracted at the continuum peak, which is indicated in the top-left panel with a red square. All maps are clipped to show only results where the integrated intensity of the CH3 CN 124 −114 is above 5σ (σ = 1.27 K kms). The purple arrow indicates the suggested direction of the outflow as in Fig. 9, while the green dashed line is perpendicular to this axis.

3.3.3 Temperature structure

Returning to the XCLASS results, the median temperature obtained from fitting CH₃CN is 132 K, with the lowest temperature near the expected sublimation temperature of ~91 K (Yamamoto 1985). This therefore gives us a resonable idea of the extent of the hot core region around W3 IRS4 (a radius of ~0.8′′, corresponding to ~1600 AU), where H₂O and other solid-phase species are thermally sublimated into the gas phase at temperatures of ~100 K.

Figure 13 shows the temperature along the filament as a function of offset from the MYSO W3 IRS4, including a comparison between those derived from CH₃CN (see Fig. 10) and H₂CO over the small region around IRS4 where both are well detected (middle-left panel). The uncertainties increase in the H₂CO temperaturemap for values well above the upper energy level of the higher excited transition (68 K). In addition, CH₃CN does not show the high-temperature spot to the SE of the hot core that is seen in the H₂CO map, though this is where the CH₃CN intensity becomes lower so the temperature from CH₃CN becomes more uncertain due to lower sensitivity. The hot spot in H₂CO could also be due to chemical effects or enhancement in the outflow, and CH₃CN likely traces denser conditions, on average, along the line of sight than H₂CO. Overall, however, the results are broadly consistent, with differences of ≲30 K around the continuum peak, giving some credence to the higher-temperature regions of the H₂CO map.

From these data it is possible to constrain, in a simplified way, the shape of the temperature profile in the envelope of IRS4. Given the spatial variation and complexity, a spatial average is not possible, but it is possible instead to take a slice in the temperature map along the direction of the filament, averaged over a beam, as shown in Fig. 13). This direction is also roughly perpendicular to the outflow. The temperature profile shows the decrease in temperature from IRS4 towards the filament, with a dip around the further cold core.

Such profiles are typically fitted with a power-law profile, as shown in the right-hand panel of Fig. 13. The inner region is relatively flat in temperature, probably due to it being close to the beam size, so this is excluded from the fitting, instead starting at an inner radius of 1200 AU. There is a change in the slope around 2000 AU, probably because the temperature averaged in a beam from 1200−2000 AU is increased by the presence of the hot spot. Beyond ~3000 AU the temperature becomes more constant, so this is excluded as well, though there’s not much indication in the map or profile that the first continuum peak has a distinct effect on the temperature (see Fig. 13).

The results for the fits for the inner, outer and whole region are given in Table 7, with the power-law exponent q for the whole region being 0.91±0.05. Values of 0.34−0.40 have been found for massive dense cores (for example Palau et al. 2014), 0.4–0.5 is typically found from hydrostatic radiative transfer modelling of the envelopes and discs of embedded YSOs (Chiang & Goldreich 1997; Beuther et al. 2007; Persson et al. 2016), observations of more evolved T-Tauri discs result in q = 0.4−0.75 (Andrews & Williams 2005), while recent ALMA results for the hot core source G31.41 found a value of 0.77 (Beltrán et al. 2018).

That a comparatively steep slope is found for W3 IRS4 compared to less dynamically complicated and lower-mass regions may be an indication that the central region is not in hydrostatic equilibrium, possibly due to the flow of material between the filament and the hot core or vice-versa. Alternatively, the temperature profile may be increased in the central region due to part of the H₂CO along the lines-of-sight close to the hot core tracing outflow or shock-heated gas, rather than the true envelope temperature, leading to a steeper profile. The profile could also be steeper compared to lower-mass sources because the higher densities lead to more efficient cooling and confinement of the region heated by the central source. Regardless of the cause, the fact that the profile is also higher for G31.41 may indicate that the requirements of forming a high-mass star lead to steeper temperature gradients than in their lower-mass cousins.

Fig. 11 Intensity-weighted velocity (moment 1, top, third, and fifth rows) and intensity-weighted FWHM (moment 2, second, fourth, and sixth rows) maps of key species in the high velocity-resolution data, along with their spectra at the continuum peak, as indicated with a red box on the moment 1 maps and a cyan box on the moment 2 maps. The lines are shown from left to right and top to bottom in order of increasing transition upper energy level (Eup). The NOEMA continuum emission with similar beam-size to the data are shown as grey contours from 5 to 45σ in 10σ intervals (σ = 0.43 mJy beam−1) in the maps, which are clipped at 5σ of the corresponding integrated intensity map. The purple arrows indicate the suggested direction of the outflow as in Fig. 9.

Fig. 12 Moment 1 (left) and position-velocity (PV) maps for molecular gas around W3 IRS4. From top to bottom, each row shows the moment 1 map and PV cuts along three different axes (middle-left, middle-right and right), for H2CO 32,2 −22,1, OCS 18-17, HNCO 101,9−91,8, CH3 CN 123 −113, and CH3CN 126 −116. The three PV cuts, taken along the outflow axis (1), rotated by 45° (2), and perpendicular to it (3), are shown on the moment 1 maps with grey lines numbered in red, where each cut is averaged over a width similar to the beam perpendicular to the cut. The 5σ contour of the data in each PV diagram is indicated with a white solid contour, while white dashed horizontal and vertical lines indicate the vLSR and zero offset from the continuum peak. The spatial and spectral resolution of the PV diagrams is indicated with a white cross in the lower-left corner of PV diagram 1 in each row. The black dashed lines indicate the results of the same cut through the moment 1 map, while the red dashed lines indicate what would be expected for Keplerian motion with a central mass of 10 and 20 M⊙. The black contours in the moment 1 maps show the continuum, as in Fig. 11.

Fig. 13 Temperature profile of H2CO as a function of offset from the MYSO and hot core W3 IRS4 along a slice which passes through the filament. Left: the position of the slice is indicated as a grey line on the temperature map from the  XCLASS fitting of H2CO (see Sect. 3.2.1 and Fig. 7). The temperature at each position is calculated as an average over the beam FWHM of the map (shown in the lower-left corner, 0.8′′), indicated with the grey dashed lines. This profile for H2CO (black), and that for CH3CN (red), are shown in the middle-left panel. The middle-right panel shows a zoom-in focusing on IRS4, with the profile shown in logarithmic units in the right-hand panel. The green dashed lines in both panels indicate three boundaries, between which we perform power-law fits to the temperature profile. The solid magenta line shows the fit over the whole region, while the cyan and blue lines show the fits to only the inner and outer regions. These latter two fits are shown as solid lines in the regions where they are fit and dashed lines where they are extrapolated to the same region as the magenta line. The results of these fits are given in Table 7. The vertical dashed black lines in the middle-left and right-hand panels indicate the beam size. The black contours in the maps show the continuum at the same resolution as the H2CO data, starting at 5σ and increasing to 45σ in 10σ steps (σ = 0.81 mJy beam−1), while the purple arrow in the middle-right panel indicates the suggested direction of the outflow as in Fig. 9.

4 Chemical conditions in the W3IRS4 region

While not as chemically rich as some other sources in the CORE sample (for example AFGL2591 or W3 OH/H₂O, see Beuther et al. 2018; Gieser et al. 2019; Ahmadi et al. 2018), the W3 IRS4 region does show a reasonable degree of chemical complexity, particularly at the position of IRS4 itself (see Fig. 14). This section begins by discussing the spatial distribution of intensity in the different species (Sect. 4.1), then discusses the results of fitting spectra at four representative positions (Sect. 4.2), before discussing the spatial abundance distributions of a few species (Sect. 4.3) and how the abundances vary with gas temperature (Sect. 4.4).

Fig. 14 WideX spectrum centred on the position of the hot core for both the combined NOEMA and 30 m data imaged with the SDI algorithm (cyan) and the NOEMA only data imaged with the Clark algorithm (black).

4.1 Spatial intensity distributions

Figure 15 shows integrated intensity maps for 22 of the brightest unblended transitions detected in the low spectral-resolution WideX data towards W3 IRS4, arranged in approximate order of chemical complexity. Eight species and the continuum show extended emission, with a range of morphologies associated with the molecular clump, outflow, hot core, and/or H II region. In contrast, the other 14 species only show compact emission associated with the hot core.

The identifications were performed by considering not only the species with transitions around the frequency ofa given line, but also by then attempting fits with XCLASS to those species to see whether other transitions within the WideX coverage were detected when they should be. For a few of the rarer species (H₂ CCO, t−HCOOH, NH₂CHO) only one transition is covered, so these identifications are considered tentative, while noting that other options at these frequencies have been checked and either ruled out or less plausible due to needing very high abundances. Most transitions have FWHM line-widths in the range 4−8 kms, and so are not well resolved by the 3 km s⁻¹ channels of the WideX data.

The filament and ring of gas emission in simpler species around the H II region W3 C have already been discussed in Sects. 3.2 and 3.1.2 respectively. The morphology of three species shown in Fig. 15, DCN, C¹⁸O, and c−C₃H₂, deserve further comment.

DCN is barely detected at the position of IRS4, instead showing peaks in cold core 1 and a point to the north of W3 C. The possible reasons for this are discussed further in Sect. 4.4

C¹⁸O is present along the filament and around the H II region, and peaks strongly to the north of W3 C, but shows no enhancement or peak at the position of IRS4. This could be due to optical depth and/or excitation due to the low upper energy of the transition (15.8 K) meaning that the J = 2−1 line is a poor tracer of the dense, warm conditions near the hot core. The enhancement for ¹³CO near IRS4 is offset slightly towards the SE, similar to H₂CO, so the ¹³CO emission near IRS4 may be outflow related, which would explain why it doesn’t show the same distribution as C¹⁸O.

c−C₃H₂, the simplest possible cyclic hydrocarbon, generally shows a very different morphology compared to all other detected species. It is distributed around the edge of the H II region (see also Sect. 3.1.2), in a linear feature across the hot core position which aligns with our proposed outflow direction, and extends along the outflow cavity to the NW. Similar to C¹⁸O, c−C₃H₂ does not show much enhancement at the position of hot core. The absence of an enhancement of c−C₃H₂ at the hot core position is complicated by the fact that the brighter of the two transitions in the covered frequency range is next to a stronger transition of ³³SO. Care was taken in calculating the integrated intensity maps for Fig. 15 to separate these two transitions, so this lack of an enhancement seems robust.

c−C₃H₂ is often detected in UV-illuminated gas in photo-dissociation regions (PDRs) and H II regions, together with C₂ H, (see for example Pety et al. 2005; Fontani et al. 2012; Pilleri et al. 2013; Guzmán et al. 2015). More recently, c−C₃H₂ has been found to trace the outflow cavity walls in two prototypical low-mass Class 0 YSOs (Murillo et al. 2018), presumably tracing gas which is illuminated by UV radiation generated either in the jet and/or outflow or from accretion. Overall, the distribution of c−C₃H₂ is consistent with with an origin in UV-illuminated gas, and provides additional support for the proposed direction of the outflow, as also seen in lower-mass Class 0 YSOs. The lack of enhancement at the hot core position could be due to shielding of the dense gas, such that only the cavity wall shows an enhancement, or if W3 C is in front of IRS4 then self-absorption in the ring around the H II region could lower the observed intensity.

Fig. 15 Integrated intensity maps of the strongest unblended line of the principle detected molecular species, as well as the continuum (bottom-right), imaged with natural weighting and resampled to a circular 1.2′′ beam. The combined NOEMA and 30 m data, imaged using this SDI algorithm, were used for those species indicated with a • before the species label, while the NOEMA only data, imaged with the Clark algorithm, were used for all others. The white contour indicates 5σ on the continuum (σ = 1.56 mJy beam−1). In the continuum panel (bottom-right), the black dashed circle indicates the ring of continuum emission associated with W3 C, the black circle indicates the position of W3 IRS4 and the black crosses indicate the two continuum peaks associated with cold cores, as shown in Fig. 2. WideX spectra are fitted in (see Sect. 4.2) at the locations of the cold cores, IRS4 and a position in the outflow with a knot of continuum emission, indicated by the green circle.

Fig. 16 XCLASS fitting results for the Clark-imaged WideX data at the hot core position. Results for the outflow and two cold core positions, as well as SDI-imaging for the hot core, are shown in Figs. C.2–C.8. The data are shown in black and the combined fit in red, while the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen). The data is cut off below the peak intensities for 13CO and SO in order to improve the visibility of the fainter lines. Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey band on the data spectrum, which extends beyond the range of the data on the right and left, indicates 3σ (σ = 0.044 K).

4.2 Chemical abundances

In order to proceed in a more quantitative fashion, spectra were extracted at four positions from the WideX data in both the NOEMA-only data which were imaged with the Clark algorithm and the merged NOEMA and 30 m data imaged with the SDI algorithm (see Sect. 2.3 and Appendix B.3). These four positions are indicated in Fig. 15, and were chosen to probe different conditions: the two cold continuum emission peaks emission in the filament, the bright knot which is also detected in the continuum in the direction of the outflow, and the MYSO IRS4.

XCLASS fitting was then performed on each spectrum to identify the emission lines and quantify column densities (see Appendix C.1). Since the continuum is detected towards each position, we can also assume that this gives a good measure of N[H₂], using the temperature derived from H₂CO for that spectrum, and thus obtain absolute molecular abundances. The resulting fits to the Clark spectrum at the hot core position are shown in Fig. 16 as an example, while the remaining fits are shown in Figs. C.2–C.8.

A total of 21 molecular species are detected towards the hot core, not including vibrationally excited transitions (see Table C.2 for the details of all detected transitions). For many species, there are not enough transitions within the frequency range of the data for accurate temperatures to be determined, so the temperature was assumed to be the same as derived for H₂CO or a more common isotopologue in the same WideX spectrum. This assumption may not hold in all cases, particularly for common species (for example ¹³CO, C¹⁸O, SO), as the intensity in each molecule and transition is affected differently by temperature, density, and abundance variations along the line of sight. However, it is probably not too bad for the majority of species compared to using a temperature which is very poorly constrained.

For SO₂ and HC₃N at the hot core position, the fits were performed including vibrationally excited lines under the assumption that these have the same properties (column density, excitation temperature, peak velocity, and FWHM) as the main species. The same assumption is also used for A-CH₃OH, where only vibrationally excited transitions lie within the frequency coverage of the data. The resulting temperatures, column densities and abundances with respectto H₂ from this fitting are presented in Table 8.

A detailed calculation of the uncertainties on the abundances is beyond the scope of this work, but we can use the work of Gieser et al. (2019, see Table 2 and Appendix B) on CORE data of another source, AFGL 2591, to gain some insight. The abundances are calculated using the column density of H₂, which is calculated using the temperature from XCLASS, the continuum intensity and an assumed opacity, and the column density of the individual species, also from XCLASS. We have been more restrictive on temperature so we assume 20% for the temperatureand 20% for the continuum intensity. The range in uncertainty in N for the species detected in AFGL 2591 vary up to ~60% for molecules with only one transition within the frequency coverage of the CORE data, such as C¹⁸O. Adding these in quadrature would give a conservative upper uncertainty on the abundances of ~70% before opacity is considered. This is difficult to constrain, but might lead to a conservative overall uncertainty of about a factor of two (that is ~100%).

Table 9 presents a number of abundance ratios for species where we detect multiple isotopologues. The spectral setup of the NOEMA observations for the CORE survey does not include ¹²CO, so in order to be able to explore the ¹²C/¹³C ratio using the observed ¹³CO abundances these are ratioed to the ¹²C abundance for the solar system from Asplund et al. (2009). Since ¹²C is less abundant than ¹⁶O in the solar abundance distribution, we make the simple assumptions that the availability of ¹²C sets the abundance for ¹²CO, that ¹²CO is the main source of ¹²C in the W3 IRS4 region, and that the cause of the change in ¹²C/¹³C with Galactocentric radius is due to variation of ¹³C rather than ¹²C.

The global ¹²C/¹³C ratio is expected to be ~76 for the W3 IRS4 region (Henkel et al. 1982; Wilson & Rood 1994; Milam et al. 2005; Qin et al. 2015). However, both the ratio of solar ¹²C to ¹³CO and H₂CO to HCO for the hot core give a ratio of <40. The determination of N[H₂] for the hot core cancels out for the ratio of H₂CO to HCO, since the same value comes into the calculation of both abundances, but does not for the ratio of solar ¹²C to ¹³CO since it only enters the calculation for X[¹³CO]. However, the resonable agreement between the two measured ratios suggests that the value of N[H₂] used for the hot core position is relatively good and therefore that the derived abundances are reasonable.

The lower ¹²C/¹³C ratio in these species towards the hot core would therefore seem to be real, though it may not be true for the whole region. Using the lower ratio for CH₃¹³CN results in a slightly worse fit, but the difference is not large enough to be conclusive as to whether this is true for all carbon-bearing species in the hot core. Low values of ¹²C/¹³C compared to those expected for the host cloud have recently been observed in complex species in two other hot core regions observed at high resolutions, IRAS 20126+4104 (~30 in CH₃OH and CH₃CN compared to an expected ratio of ~60–70, Palau et al. 2017) and G31.41+0.31 (≲10 in CH₃CN and CH₃OCHO compared to an expected ratio of ~47, Beltrán et al. 2018). Such low values may therefore be a feature of hot cores.

A possible explanation is that, similar to the case of deuterium enhancement, ¹³ C can be enhanced in CO at low temperatures due to the reaction (11)

and depleted in species which form through reactions with C⁺ (for example Langer et al. 1984). Low-temperature gas-phase enhancement of ¹³CO relative to ¹²CO is preserved in the ices in cold regions where CO is frozen out, which is then passed on to species formed from CO and maintained for some time once the ices are sublimated. H₂CO can be formed at warmer temperatures in the gas-phase via (12)

or at cold temperatures on grain surfaces via sequential hydrogenation of CO a process that also leads to the formation of CH₃OH (for example Walsh et al. 2014): (13)

The similar enhancement of ¹³C in CO and H₂CO, with the ratio in CO slightly lower than in H₂CO as predicted by Langer et al. (1984), is consistent with a cold temperature origin for H₂CO, formed in the solid phase on dust grains.

Indeed, the ratio of E- to A-type methanol at the hot core position tells a similar story. These two versions of methanol have different nuclear spin alignments in the hydrogen atoms of the methyl group, effectively acting as two distinct species because there are no radiative or collisional transitions between the different types. The ratio of E- to A-type CH₃OH is therefore set during formation, with a ratio of 1 indicating formation at temperatures above ~37 K and colder temperatures resulting in an overabundance of A-type methanol (Wirström et al. 2011). The value of E/A ~0.3 obtained for the hot core position, is equivalent to a nuclear spin temperature of ~6 K, again indicating a cold-temperature origin for the observed CH₃OH in W3 IRS4. The ratio of A-CH₃OH ν_t = 1 to A-CH₃OH ν = 0 is unlikely to be greater than 1. If it is less than 1 as assumed here then this value of E/A is an overestimate.

The abundance ratio of ¹³CO/C¹⁸O increases slightly from cold core 1 (4.5), to cold core 2 (6.1), to the hot core (7.7), that is with increasing temperature of the region, though the ratio for the outflow position (5.6) is slightly lower than for cold core 2 despite having a warmer temperature. The values all lie between the solar value (3.0) and the value expected for W3 of ~7.3 (Wilson & Rood 1994). They are all below those observed for PDR regions of Orion (≥9.6; Shimajiri et al. 2014) where preferential destruction of C¹⁸O from isotope-selective photodissociation is thought to play a role. Ultimately, the differences are small enough that it is not clear whether they are truly significant or not given the likely uncertainties in the data and fitting.

The large solar ¹²C to ¹³CO ratio at cold core 1 (~400), and to a lesser extend for cold core 2 (~180), as well as the fact that the abundances for all species are approximately one to two orders of magnitude lower than those at the warmer positions, is consistent with freeze-out of CO onto the dust grains, particularly since the temperature there is close to the CO freeze-out temperature, depending on the ice composition and structure (Anderl et al. 2016; Fayolle et al. 2016).

The lower solar ¹²C to ¹³CO abundance ratio at the outflow position (~11) could be due to N[H₂] being underestimated (only ~13% of the single-dish continuum flux is recovered in the 1.4 mm NOEMA continuum data, Beuther et al. 2018). Alternatively, if the result is a true reflection of the ratio, this might indicate that the 1.4 mm continuum is not tracing the full column of molecular hydrogen, perhaps due to destruction of dust grains in the outflow cavity. Isotope-selective photodissociation is another possibility, but in that case one might expect a larger ¹³CO/C¹⁸O ratio than is observed. If the ¹²C/¹³C ratio for this position is actually closer to the ~35 for the hot core position or the ~76 suggested by cloud-scale observations, then this would imply that our measurement of N[H₂] is underestimated by a factor of ~3.5 or 7.6, respectively. This would lead to the abundances at this position decreasing by the same factor.

The ratios for the various sulphur-bearing species all fall well below the solar ratios (see Table 9), but the values for SO₂ are similar to those obtained by Esplugues et al. (2013) for various regions around Orion KL. As in that study, our value for ³⁴SO₂/³³SO₂ = 5.8 is close to the solar abundance ratio of 5.6. The extremely low value for SO/³³SO is likely due to optical depth effects in the main isotopologue, but it seems less likely that SO₂ is suffering from such issues. It is possible that the isotopologues of SO₂ experience selective photodissociation due to self-shielding effects (Lyons 2007), which may account for the observed ratios. However, since sulphur chemistry is still relatively poorly understood and shocks may contribute to some or all of the emission, it is difficult to be more conclusive.

The abundances of a number of species are comparable to values found for other hot core regions. For example, Palau et al. (2017) found similar abundances for SO, HNCO, HCO, HCOOH, and CH₃CN towards IRAS 20126+4104. However, the values for H₂CCO, CH₃OH, CH₃OCHO, and CH₃COCH₃ are approximately an order of magnitude higher in IRAS 20126 (Palau et al. 2017), and both CH₃CN and CH₃OCHO are more than one order of magnitude higher in G31.41+0.31 (Beltrán et al. 2018) than in IRS4, suggesting some difference in conditions or age between these sources.

Fig. 17 Abundance distribution of H2CO (left), DCN (middle), and CH3CN (right). The black contours in the maps show the continuum smoothed to the same resolution as the line data, starting at 5σ and increasing to 45σ in 10σ steps (σ = 0.81 mJy beam−1). The beam is shown in each panel and a scale-bar is shown in the middle panel. Abundances are only calculated where both the continuum emission and integrated intensity of the H2CO 303 −202 and DCN 3−2 lines are above 3σ and where CH3CN 124 −114 is above 5σ. The uncertainty in the abundances is conservatively estimated to be around a factor of two (see text for details).

4.3 Spatial abundance distributions

The maps of N[H₂] and N[p-H₂CO] can be used to produce a map of the abundance of H₂CO by assuming a value for the ortho-to-para ratio for H₂CO. If H₂CO is produced in the gas-phase, this should be equal to the thermal limit of 3 while a lower value, typically in the range 1.5−2, is expected if the formation takes place in ice on the grain surface, followed by sublimation to the gas-phase (see for example Kahane et al. 1984; Dickens & Irvine 1999; Jørgensen et al. 2005). Given the suggestion from the H₂CO/HCO ratio that the formaldehyde in the hot core was formed in cold conditions on the ice, we assume an ortho-to-para ratio of 1.6, as obtained by Jørgensen et al. (2005). Our abundances would increase by a factor of 1.5 if the high-temperature value of 3 was used instead.

The resulting map is shown in the left panel of Fig. 17. X[H₂CO] varies between 2.4 × 10⁻⁹ and 1.1 × 10⁻⁷ around the hot core IRS4 with a median of 1.0 × 10⁻⁸. The median abundance over the rest of the filament is 7.4 × 10⁻¹⁰, though these are upper limits since the H₂CO data contain the large-scale information but the continuum (used for N[H₂]) do not. There is a clear increase in abundance towards the hot core and lower abundances towards the colder continuum peaks, supporting the idea that H₂CO is indeed formed on the dust grains and thermally desorbed near IRS4.

These can be compared with the results of Gerner et al. (2014), who observed a large sample of high-mass star forming regions with the IRAM 30 m and obtained abundances for a number of species, including H₂CO. They divided these sources into four classifications, infrared dark clouds (IRDCs), high-mass protostellar objects, hot cores, and UCH II regions, and found median abundances for H₂CO of 5−7 × 10⁻¹⁰, 5−7 × 10⁻¹⁰, 2−9 × 10⁻⁹, and 0.7−5 × 10⁻⁹ respectively for these stages.

Recreating this analysis for W3 IRS4, we take the peak 850 μm continuum flux for W3 IRS4 from the SCUBA legacy archive catalogue (14.99 Jy, Di Francesco et al. 2008) convert to the expected flux at 220 GHz using a ν⁴ scaling, and use this to calculate N[H₂] = 6.2 × 10²² cm⁻² using the temperature from the fit to the IRAM 30 m H₂CO data (67.5 K, see Sect. 3.2.6).

The single-dish observations for W3 IRS4 then result in an abundance of X[H₂CO] = 4.1 × 10⁻⁹, in line with the results for hot cores from Gerner et al. (2014). The values obtained for the filament from the combined NOEMA and 30 m data are roughly in line with the IRDC results from Gerner et al. (2014). However, the single-dish abundance is a significant underestimate for the hot core and a factor of five higher than the abundance in the filament (cf. similar size-scale issues with temperature in Sect. 3.2.6). Gerner et al. (2014) found that their hot core models produced considerably higher abundances than they observed, but the abundances we obtain for the hot core from the combined NOEMA and 30 m data agree with their models for hot cores at some timescales. Thus the discrepancy found by Gerner et al. (2014) between their observations and models for H₂CO is due to a mixing of conditions within the single-dish beam, which results in an intermediate value between the colder and warmer gas.

The abundance map of DCN (middle panel of Fig. 17) was calculated using the temperature map from H₂CO and the integrated intensity in the DCN J = 3−2 transition. It shows a decrease near the cold density peaks and a similar spread of ~1.5 dex over the map compared to H₂CO. The values are broadly within the range predicted by the models of Gerner et al. (2015). While the intensity of DCN is too low to be well detected at the hot core position in this higher-resolution map, the transition is detected in the WideX spectrum in a 1.2′′ beam (see Fig. 16 and Table 8).

The abundance of CH₃CN varies in the map between 9.5 × 10⁻¹⁰ and 9.6 × 10⁻⁹ with a median of 2.7 × 10⁻⁹, slightly below the value derived using the WideX spectrum at the hot core position (see Table 8). Both the higher and lower extreme values lie near the edges of the abundance map, where the intensity of the emission is lower and so there is likely more uncertainty in both the column density and the temperature. The abundances relative to CH₃OH (~10⁻²) are similar to those found for other hot cores and MYSOs (see for example Fayolle et al. 2015) and in the chemical models of hot cores by Gerner et al. (2014).

Fig. 18 Abundance vs. temperature from the maps shown in Fig. 17 (dots) and the fits to the WideX spectra discussed in Sect. 4.2 (connected crosses). The different colours denote different species (see legend), using the same colours as in Fig. 16. The dashed lines indicate a rolling median of map points within ± 10 K, measured every 10 K as long as at least 10 points are included. The downward vertical arrows show the effect of increasing N[H2] at the outflow position (which is likely underestimated, leading to the molecular abundances being overestimated) such that the 12C/13C ratio would be 76. The uncertainty in the abundances is conservatively estimated to be around a factor of two (see text for details).

4.4 Variation of abundance with temperature

Figure 18 presents the abundance distribution from the maps in Fig. 17 as a function of temperature. Also shown are the results from fitting the WideX spectra from Table 8 for those species which are detected in at least two positions, that is cold core 1, cold core 2, the outflow position, and the hot core IRS4. Taking into account the possible overestimation of the abundances for the WideX spectrum at the outflow position due to questions over the determination of N[H₂] (see Sect. 4.2 and the downward arrows in Fig. 17), most species show an increasing abundance with increasing temperature. The abundance of H₂CO rises more quickly below ~90−100 K than above it, and CH₃CN, CH₃OH, HNCO, c−C₃H₂, and SO allshow increases of ~1−2 orders of magnitude between positions below and above this temperature. Indeed, the enhancement in SO is a lower-limit due to the likelihood of optical-depth issues at the hot core position (see Sect. 4.2).

In the case of H₂CO, CH₃CN, CH₃OH, and HNCO, this is entirely consistent with their being locked into the water ice mantles on dust grains and being released as water sublimates at ~100 K (for example Fraser et al. 2001; Burke & Brown 2010). The gas-phase abundances in the cold gas are likely maintained by non-thermaldesorption mechanisms, possibly through cosmic-ray induced photodesorption similar to or together with water (see for example Caselli et al. 2012; Mottram et al. 2013; Schmalzl et al. 2014), or via chemi-desorption. That the H₂CO abundance rises relatively smoothly with temperature, rather than sharply increasing at a certain temperature, suggests that the thermal sublimation process is relatively gradual. The enhancements in c−C₃H₂ and SO maybe due to the stronger UV radiation and shocks close to the base of the outflow and the MYSO powering the hot core, since they are not typically thought of as being formed in ice mantles.

The details for ¹³CO and C¹⁸O depend on the correction needed for N[H₂] at the outflow position, but their trends are broadly consistent with the lower sublimation temperature of CO (~20−40 K, for example Fayolle et al. 2016; Anderl et al. 2016). DCN and HC₃N show modest variation with temperature at best, suggesting that their abundances are maintained through gas-phase reactions which are efficient even at lower temperatures down to ~40−50 K. In particular, the abundance of DCN at the hot core position is entirely consistent with that for the colder gas, indicating that the lower intensity in the J = 3−2 transition around the hot core (detected in the WideX spectrum at 1.2′′ resolution in Fig. 16 but below the sensitivity limit imposed on the 0.8′′ resolution map in Fig. 17) is primarily an excitation effect due to the low E_up (~21 K) of this transition, rather than destruction at higher temperatures.

5 Discussion

5.1 Links between spatial scales

The CORE data and the analysis presented in Sects. 3 and 4 have shown that the W3 IRS4 region contains a number of structures which span a range of scales and environments. Table 10 gives a qualitative breakdown of the species which trace each physical component of the system, as well as approximate physical properties as determined in the previous sections. While each region can be identified through a combinationof physical, chemical, and spatial differences, there is also considerable overlap in properties between them, with no indication of sharp boundaries between clump and filament or filament and core. This suggests that any flow from clump to core scales does not cause a shock, and may instead serve to blur or wash out any transitions.

The kinematics do suggest some interplay between the different components along the filament, such that they cannot be treated entirely in isolation from the neighbouring and larger-scale structures. The likely evolutionary consequences are discussed in Sect. 5.2. This does, however, suggest an important note of caution for those interested in disc or core-scales in active star formation regions. If such studies were only to include information on smaller spatial scales then they might incorrectly asume a more static and isolated picture than is revealed when larger-scales are included. Furthermore, the properties obtained from such spatially filtered data may be biased in an analogous way to what happens when only large scales are observed with single-dish observations. If we truly seek to obtain a comprehensive picture of star formation then this requires observations with information on a range of spatial scales.

5.2 The past, present and future of the W3 IRS4 region

When considering the different evolutionary stages of the different sources within the region, the O-star powering the UCH II region W3 Clikely formed first by virtue of it’s more advanced evolutionary state, potentially fed by the filament which could originally have extended further to the SW. The formation and expansion of the UCH II region could have partially disrupted the filament and led to the formation of W3 IRS4 at the point where these to flows balance between the flow along the filament and the expanding H II region.

W3 IRS4 does not power an H II region despite being luminous enough (L_bol ~2.2 × 10⁴ L_⊙, corresponding to a ~20 M_⊙ O9.5 star). The central star may be a bloated beyond what would be expected for a main-sequence star of the same luminosity and/or mass (see for example Hosokawa et al. 2010; Kuiper & Yorke 2013; Palau et al. 2013; Haemmerlé & Peters 2016), a process which can affect stars between ~5 and ~30 M_⊙ experiencing rapid accretion, consistent with the observed large mass infall rate of 10⁻⁴−10⁻³ M_⊙ yr⁻¹ towards IRS4. Alternatively, the H II region may be quenched by the large infall rates (see for example Peters et al. 2010b,c).

It is clear from the chemistry that the gas around W3 IRS4 spent a significant amount of time at cold enough temperatures that not only H₂ O but also CO was frozen out onto the grain mantles (that is T ≲ 35 K). However, the bulk temperature in the dense gas within the region is currently warm enough that CO should be in the gas phase. Given that star formation has been ongoing in W3 in general and near W3 IRS4 in particular for some time, as evidenced by the numerous older H II regions and OB stars nearby (for example Tieftrunk et al. 1997; Oey et al. 2005; Bik et al. 2012), the radiation field experienced by the filament has probably not changed much in the past 10⁵−10⁶ yr. Since the timescales for high-mass prestellar cores has been observationally constrained to be very short (10⁴−10⁵ yr, see for example Tackenberg et al. 2012), the whole clump or cloud may have spent a significant amount of time in a cold dark phase before star formation began.

Alternatively, material with a cold origin is being deposited in and around W3 IRS4 through flow along the filament from beyond our field of view. As discussed in Sect. 3.2.5, the flow is certainly fast enough that material from cold cores 1 and 2 will probably reach the envelope of IRS4 within ~2−6 × 10⁴ yr. Either way, this provides reassurance that the standard practise of having a long cold, dark prestellar phase when running chemical models is reasonable, even in complex environments.

5.3 The W3 IRS4 region in context

It is now important to place W3 IRS4 within the context of observations of other regions, to consider why any differences may occur, and to explore what implications these might have for our understanding of star formation as a whole. W3 IRS4 shows a steeper temperature profile in the north-west part of the surrounding envelope than those typically found for embedded sources (cf. Chiang & Goldreich 1997; Beuther et al. 2007; Palau et al. 2014; Persson et al. 2016). Furthermore, the material close to IRS4 is much more disrupted and confused down to the resolution of the NOEMA observations (~600 AU) when compared to other sources of similar luminosity and observed at similar resolution, where stable structures or at least consistent kinematics between different species can be seen on scales up to ~1000 AU (for example Beuther et al. 2012; Johnston et al. 2015; Ilee et al. 2016; Cesaroni et al. 2017; Ahmadi et al. 2018).

The high mass-infall rate from the large-scale, one-sided flow may explain both of the steeper temperature profile and the more disrupted nature of the material around W3 IRS4, to the point where there is not indication of any disc-like structure down to the spatial resolution of the observations (that is 600 AU). Flow along filamentary structures has been seen in other star formation regions (for example the hub-filament system SDC13, Peretto et al. 2014; Williams et al. 2018) or the spiral flows towards W33A (Maud et al. 2017), suggesting that W3 IRS4 is unlikely to be unique, though the one-sided nature of the flow may make it more extreme.

It has been suggested in both low and high-mass star formation contexts that the rotationally stable part of discs around YSOs may start out very small and grow over time (for example Stahler et al. 1994; Williams & Cieza 2011; Kuiper et al. 2011). An alternative, though not mutually exclusive, reason for the non-detection of a disc in IRS4 is therefore if it is younger than the sources where disc-like structures have been detected. Indeed, Csengeri et al. (2018) observed a young source (G328.2551-0.5321) with ALMA which is in the process of forming a massive stars, where the disc is only seen on ~250 AU scales, much smaller than the 600 AU resolution of our NOEMA observations. Testing such a hypothesis for IRS4 will require observations with higher resolution by a factor of at least >3 in order to be conclusive, though comparative chemical modelling of the different sources within the CORE survey would also help to establish whether there are any grounds to believe IRS4 is truly younger than W3 H₂O, for example.

Whatever the cause, the dynamics of the region imply that it will produce fewer, more massive stars than it would if there were no significant flow along the filament, which might allow further fragmentation and star formation to proceed in the cold cores. Thus, the dynamics of the environment on intermediate scales have the potential to dramatically affect the resultant stellar distribution on clump to core scales given the same starting mass. Indeed, suppressing fragmentation and increasing the mass reservoir via rapid mass infall along filaments may well be needed to form the most massive stars.

6 Summary and conclusions

This paper has presented the data, results and analysis of 1.4 mm continuum and molecular line NOEMA and IRAM 30 m observations of the high-mass star forming region W3 IRS4 as a case-study for the potential of such data from the CORE survey (Beuther et al. 2018). Our conclusions are summarised below.

Firstly, inclusion of short-spacing observations are crucial for correctly understanding the environment and context within which high-mass star formation takes place. This is true not only for correctly constraining the intensity, but also the structure and kinematics of the available gas reservoir. Attention should also be paid to the assumptions implicit in using certain imaging algorithms (for example point vs. extended emission).

Secondly, the W3 IRS4 region contains a range of evolutionary stages of high-mass star formation (cold filament and cores, MYSO with surrounding hot core, UCH II region) on a range of spatial scales from 0.01−0.1 pc which are interacting with each other and their clump-scale (≥0.1 pc) environment, affecting each other’s evolution.

W3 C, ~3′′ to the SW of W3 IRS4, is an UCH II region powered by a ~16 M_⊙ O9.5 star, though the location of the powering source is now in doubt. It exhibits a ring of ionised free-free continuum emission at 1.4 mm which is surrounded by a gas ring, detected in ¹³CO, C¹⁸O, H₂CO, c-C₃H₂, and SO. The radial intensity distribution of these species suggest that isotope-selective photodissociation by FUV photons plays a role in the structure of the surrounding PDR. This may be compressing some of the gas on the SW side of W3 IRS4.

There is a large velocity gradient along the cold filament to the north-east of W3 IRS4, seen in a number of species. This filament has a mass-per-unit-length of 229 M_⊙ pc⁻¹, median temperature of ~50 K, median volume density of 2 × 10⁶ cm⁻³, and is thermally unstable to gravitational collapse but may be supported by turbulence. Gas within the filament is flowing towards W3 IRS4 with a mass infall ratio of 10⁻⁴−10⁻³ M_⊙ yr⁻¹ (depending on inclination), fast enough that the few cold cores within it are unlikely to form stars before being accreted.

W3 IRS4 is a MYSO embedded within a hot core, possibly in a bloated phase, since the radio emission associated with it is too weak and more consistent with being formed in a jet. The axis of the jet and/or outflow is suggested to be in the plane of the sky in the NE-SW direction. The velocity gradient seen in a number of species around the W3 IRS4 is clear, but the axis changes depending on the species and so likely traces different combinations of infall, rotation and outflow, depending on the tracer. Any disc around W3 IRS4 must be smaller than 600 AU.

The temperature in the centre of W3 IRS4 is warm, with peak temperatures of 160−170 K in most tracers. We derive a radial temperature profile which has a power-law exponent of between q = 0.6 and 1.1 depending on the radii used. This is steeper than typically found in hydrostatic simulations, possibly due to the high asymmetric infall of material from the filament.

Finally, the chemical composition of the gas around W3 IRS4 suggests that it originated in regions where CO was frozen out, despite the median temperature in the surroundings being above the CO sublimation temperature. Most species show significant enhancements above ~90−100 K, the exceptions being DCN and possibly HC₃N, which therefore likely have a gas-phase origin. The most chemically complex species detected is acetone (CH₃COCH₃), while the spatial distribution of c−C₃H₂ stands out from the other species, likely tracing UV-illuminated gas, as suggested for PDRs and outflows from low-mass Class 0 protostars.

Taken together, these observations and analysis have revealed the ongoing star formation in the W3 IRS4 region to be even more inter-connected than previously thought, revealing the wealth of information that can be obtained when considering the links between clump and disc scales and significantly improving the interpretation of observational results for this region. We have also been able to map the physical conditions, particularly temperature, over large regions for the first time, providing more insight into their variation and more accurate mass and abundance determinations. More generally, this study has revealed the potential that dynamics have to change the distribution and properties of star formation within a region, providing a warning to observations which only focus on the continuum or highest resolution. A number of the lines of enquiry explored in this paper (for example large-scale temperature maps from H₂CO) will provide fruitful grounds for further studies when expanded to the to full CORE sample.

Acknowledgements

The authors would like to thank the anonymous referee for their helpful comments, which improved the paper, the IRAM and GILDAS staff for their help during observations and data reductions, Thomas Möller for his technical support with XCLASS, Arjan Bik for helpful correspondence about the ionising source for W3 C, and Ewine van Dishoeck, Lars Kristensen, Magnus Persson and Yuan Wang for fruitful discussions. J.C.M., H.B. and A.A. acknowledge support from the European Research Council under the European Community’s Horizon 2020 framework programme (2014-2020) via the ERC Consolidator grant “From Cloud to Star Formation (CSF)” (project number 648505). T.Cs. acknowledges support from the Deutsche Forschungsgemeinschaft, DFG via the SPP (priority programme) 1573 “Physics of the ISM”. R.K. acknowledges financial support via the Emmy Noether Research Programme funded by the German Research Foundation (DFG) under grant no. KU 2849/3-1. A.P. acknowledges financial support from CONACyT and UNAM-PAPIIT IN113119 grant, México. REP is supported by an NSERC Discovery grant. S.E.R. acknowledges support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement # 706390. This work is based on observations carried out under large project number L14AB with the IRAM NOEMA Interferometer and the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This research made use of MATPLOTLIB (Hunter 2007), APLPY, an open-source plotting package for PYTHON (Robitaille & Bressert 2012), and Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2018).

Appendix A Details for the IRAM 30 m data

This section provides details of the IRAM 30 m data reduction (Appendix A.1), as well as an overview of the data for W3 IRS4 (Appendix A.2).

Fig. A.1 IRAM 30 m spectra at 11′′ resolution towards the hot core position of W3 IRS4, resampled to 3 km s−1 channels. The location of the main detected transitions within these frequency ranges are labelled and indicated with red dashed vertical lines (see Table C.1 for details). The location of SiO 5−4, which is not detected, is indicated in blue. The frequency axis has been shifted to the rest velocity using a source velocity of vLSR = −44.5 km s−1.

A.1 Reduction

The IRAM 30 m spectral line maps were reduced and processed using the following steps. First, each spectrum was converted from antenna temperature () to main-beam temperature(T_mb) using values of η_mb extrapolated from the measurements of 0.670 at 210 GHz and 0.641 at 230 GHz provided on the IRAM website⁴ to the central frequency of each window. The specific values used were 0.663, 0.657, 0.640, and 0.635 at 215, 219, 231, and 234 GHz, respectively.

Next, initial line-free regions in each spectrum were identified by excluding all channels in that spectrum above 5σ_rms. Any group of blanked channels that was narrower than 50 channels (~13 km s⁻¹) was increased to 50 in an attempt to account for faint line-wings. The detector covers each window using three frequency units, so a rough baselining was performed separately for each unit in order to remove any platforming or discontinuities in the baseline of the spectrum between the different units. Polynomials of order 0 to 9 were attempted, with a higher order polynomial fit only accepted as superior if σ_rms decreased by more than 10% compared to a lower-order fit. For most spectra, this resulted in a 0th or 1st order polynomial being used.

Following rough baselining, new line-free regions were identified in the same way as before, but excluding all channels above 3σ_rms and with a size minimum of 30 channels (~8 km s⁻¹). A deeper search was used for this second pass because the worst of any baseline artefacts, which might otherwise be flagged as emission, should have been removed in the previous steps. The baselining was then repeated on the original data using the updated regions.

Spikes were then identified as single channels which were above 4σ_rms with their neighbours being below 4σ_rms. These were replaced with noise. The strong ¹²CO and ¹³CO 2−1 lines caused artefacts due to photon leaks into the opposite sideband, mirrored in frequency around ~225 GHz. Small regions around these artefacts were therefore also replaced with noise.

Additional spectral artefacts were identified where the difference between the H and V polarisation spectra for the same position in the same scan differed by more than 10σ_rms. These channels were also replaced by noise. The spectra for a given source and frequency window were then collected together and gridded to produce final output maps.

The final product of this process is a position-position-velocity cube for each source in each of the four 4 GHz spectral units, centred at approximately 215, 219, 231, and 234 GHz respectively.

A.2 Overview of the data

This section provides an overview of the IRAM 30 m observations towards W3 IRS4. Figure A.1shows the spectrum for the central position for each of the four frequency bands. The position of prominent or expected lines are indicated, with details given in Table C.1. The broad absorption features which are due to atmospheric absorption. Figure A.2 presents moment-maps for the brightest transitions.

Fig. A.2 Overview of the IRAM 30 m data for selected lines, one per row. For each transition the panels show a spectrum at the source centre (left), as well as maps of the integrated intensity (moment 0, middle-left, between −87.5 and 2.2 km s−1), intensity-weighted velocity (moment 1, middle-right), and intensity-weighted line FWHM (moment 2, right). The position of the spectrum is indicated on the left-hand panel by a black square. The contours in the three maps indicate the NOEMA 1.4 mm continuum with uniform weighting, as in the right-hand panel of Fig. 2. The beam-size and a scale bar are given in the middle-left panel. For the middle-left and middle-right panels, values are only calculated where the spectrum has a peak intensity above 3σ, with grey used to indicate those regions that do not meet this requirement.

Fig. A.2 continued.

Appendix B Data combination and imaging

This section provides more details on the combination of the NOEMA and IRAM 30 m data (Appendix B.1) and on the choices and underlying assumptions made during the imaging process (Appendix B.2). A further subsection (Appendix B.3) provides a detailed discussion of an example where the combination of both strong point-like and significant extended emission make the choice of approach more difficult.

B.1 The combination of interferometric and single-dish data

For those transitions with extended emission, the 30 m data are used to provide short-spacing information for the NOEMA observations. In general there are two approaches to this task: to convert the single-dish flux to the uv-plane and combine with the interferometric visibilities before imaging the combined uv-plane data; or image the interferometric data alone and then feather them with the single-dish map through a linear combination in Fourier space. We use the former approach, as implemented using the GILDAS task UVSHORT within the MAPPING module.

In this method, the 30 m data are fits manually resampled to the same frequency channels as the NOEMA data, after which the flux units are converted from main-beam antenna temperature to Jy beam⁻¹. The re-gridded and re-scaled data are then Fourier-transformed into the uv-plane, deconvolved of the 30 m single-dish beam and Fourier-transformed back to the image plane where they are then convolved with the primary beam of the 15 m NOEMA dishes. The result is then Fourier-transformed back into the uv-plane again and sampled onto the same regular physical grid as the interferometric visibilities. Finally, the interferometric and short-spacing visibilities are merged, ready for imaging.

Extended emission is filtered out by the interferometer, which can lead to artificial negative features or artificially low flux values in regions where there is genuine self-absorption on small spatial scales. The inclusion of short-spacing information from the 30 m data can therefore have a marked effect on the resulting map, increasing flux levels and removing negative features where the interferometer fails to recover the larger-scale emission correctly. Indeed, as can be seen in the example in Fig. B.1, in some cases including the single-dish data not only changes the flux levels but also the spatial distribution of the emission.

If the emission is compact, and so correctly recovered by the interferometer, then adding the short-spacing information does not improve the data. Indeed,the line may not be detected in the single-dish data due to beam-dilution, and so one is essentially combining the interferometric visibilities with noise. The resulting merged data may then have slightly higher noise than the non-merged version and the peak flux may be decreased as emission is redistributed during imaging into a very weak but extended component.

For this reason, we compare the merged and NOEMA-only data for all lines and only use the merged data if the line is well-detected in the 30 m data and the merging results in a clear improvement in the quality of the final image. Where the merged data are used, through comparison of the original 30 m fluxes with the final merged maps, we estimate that >90% of the total flux, as measured in the single-dish observations, is recovered in the combined image.

Fig. B.1 Integrated intensity (0th moment, integrating between −55.3 and −29.8 km s−1) maps of H2CO 30,3 −20,2 imaged with the SDI algorithm (see text for more details) and uniform weighting for NOEMA only (left) andmerged NOEMA and 30 m data (right). The black-solid and white-dashed contours indicate 5σ and −5σ, where σ = 12.5 and 12.6 mJy beam−1 km s−1 for the NOEMA only and the merged data respectively. The black circle, dashed circle and crosses indicate the positions of W3 IRS4,W3 C and the two cold cores, as indicated in Fig. 2.

B.2 Imaging interferometric and merged visibilities

Two further choices during the imaging process have a bearing on the resulting data. Firstly, the weighting given to visibilities on different baselines effects the beam-size of the final image and the noise level of those images, essentially trading spatial resolution against sensitivity (for example uniform or natural weighting, see Fig. 2). Secondly, the cleaning and imaging of interferometric visibilities requires an assumption of how the emission is distributed (for example whether most of the emission is point-like or extended). This latter choice can have a non-trivial effect on the resulting maps.

For example, Fig. B.2 shows a comparison of two methods implemented in GILDAS: Clark (Clark 1980); and Steer-Dewdney-Ito (SDI; Steer et al. 1984). Both clean algorithms start from the dirty image and work by iterating over the following steps until some criteria is met: first, source emission is identified in the residual image; next, the identified emission is removed from the residual image in Fourier space after convolution with the beam to produce a new iteration of the residual image; the subtracted source emission is then added to the “clean” image. The principle difference between these two methods is that the Clark algorithm assumes a point-source centred on the brightest pixel in the current residual image while the SDI algorithm includes all emission above a threshold based on the value of the brightest pixel in the residual image. SDI therefore assumes that the emission is extended.

In the case of the compact OCS emission, the SDI cleaning method distributes the flux over a larger area than the Clark, as is to be expected, but given the compact nature of the emission the Clark imaging is the best choice. In contrast, the Clark imaging of C¹⁸O emission results in a much more broken-up appearance of the emission compared to the SDI imaging, which shows a smoother large-scale distribution. The choice of imaging method therefore has a strong impact on the resulting maps in cases where there is strong extended emission. In this case it is clear that the assumption of point-like emission inherent in the Clark algorithm does not hold and so the SDI method is likely to return a result which is closer to the “true” distribution of emission in this line. An example of a more complex case, where both extended and strong point-like emission is present (SO 6₅ −5₄ in W3 IRS4), is discussed in more detail in Appendix B.3.

The middle and right-hand panels of Fig. B.2 also show the impact of correcting the data for the decreasing sensitivity of the primary beam of the interferometer. Correcting for the Gaussian decrease in intensity caused by the primary beam leads to an increase in the flux density of emission around the edge, but also leads to increasing noise towards the edge of the map. The corrected data are therefore typically truncated at some radius, in this case at the FWHM of the primary beam (~23′′).

B.3 Imaging lines with both extended and bright point-like emission

As discussed in Appendix B.2, the two imaging algorithms employed in this paper, Clark and SDI, each invert and image the interferometric visibilities assuming either that the emission is point-like or extended, respectively. For the majority of lines, this provides good results. However, for some transitions there is both extended and strong point-like emission, and so care must be taken to check the imaging results for artefacts induced by this more complicated situation. The most extreme example within the CORE NOEMA data for W3 IRS4 is the SO 6₅ −5₄ transitionataround 219.95 GHz, for which Fig. B.3 presents integrated intensity maps of the data after imaging with either the Clark or SDI methods, with and without including the 30 m short-spacing information, and for three different weightings (uniform, natural and an intermediate option).

Fig. B.2 Integrated intensity (0th moment) maps of OCS 18−17 (left, between −55.3 and −29.8 km s−1) and C18O 2−1 (middle and right, between −81.8 and −0.8 km s−1) for data images with the Clark (top) and SDI (bottom) clean algorithms. All are imaged with uniform weighting. The OCS map is from only NOEMA data while the C18O is merged NOEMA and 30 m data. The C18O maps are without (middle) and with (right) primary beam correction, where in the latter case the data have been truncated at the FWHM of the primary beam. The black contours indicate 5σ, where σ = 13.8 (OCS, Clark), 12.6 (OCS, SDI), 19.0 (C18O before primary-beam correction, Clark), and 20.0 mJy beam−1 km s−1 (C18O before primary-beam correction, SDI). The black circle, dashed circle and crosses indicate the positions of W3 IRS4, W3 C and the two cold cores, as indicated in Fig. 2.

Fig. B.3 Comparison of the integrated intensity of SO 65−54 between −81.8 and −0.8 km s−1 with different imaging choices. The left and middle-left columns show results with only the NOEMA visibilities while the middle-right and right columns show the results with the 30 m short-spacing information included. The maps in the left and middle-right columns were imaged with the Clark algorithm while those in the middle-left and right columns were imaged with the SDI algorithm. The top row was imaged with uniform weighting, the middle row using a GILDAS-robust parameter of 3 to give an intermediate beam-size, while the bottom row using natural weighting. The intensity scale is chosen to show the extended emission, with the peak intensity for the hot core being between ~240 and ~790 K km s−1, depending on the method.

This transition is strongly detected and extended in the IRAM 30 m data (see Fig. A.2) and has extended emission in the NOEMA visibilities even if the 30 m data are not included, though not including the short-spacing information leads to less extended emission and considerably more negative features, particularly near the hot core. The short-spacing information is therefore required in the imaging. The SDI maps show some ring-like features centred on the very bright, compact emission from the hot core which are not seen in the Clark maps, which are likely a result of side-lobes from the point-source not being properly accounted for in the SDI cleaning. The Clark imaging still introduces a mottling or break-up of the extended emission due to the point-source assumption (see Sect. 2.3), but comparison of the two makes clear that the majority of the emission structures are real.

Ideally, one would like to use a hybrid scheme for such data which started with a small number of iterations of the Clark algorithm to minimise the side-lobes of any point sources, then continue with the SDI method to better recover the extended emission structure, but such a solution is not currently implemented in GILDAS and beyond the scope of this paper to develop. We therefore choose to use the SDI imaging in our analysis of the extended emission (see Sect. 3.1.2). However, the regular peaks seen in the radial intensity profile for SO in Fig. 4, and possibly H₂CO to a lesser extent, are likely due to this imaging artefact. The “true” radial intensity profile is likely smoother, since flux is conserved during the imaging.

Appendix C Detected lines and XCLASS fitting

C.1 XCLASS setup and assumptions

Spectral line data are fitted by the XCLASS software under the assumption that: the emission comes from an isothermal slab, so there are no density gradients and one temperature can reproduce all transitions of a given species; the kinematics (that is peak velocity and FWHM) of all transitions of a given molecule are the same; the molecules are in local thermodynamical equilibrium (LTE) with no contribution from radiative excitation, such that the excitation temperature is equal to the kinetic temperature of the gas; and that the emission is resolved such that the source size is fixed to a value which is larger than the beam (typically 20′′).

The remaining four free parameters (peak velocity, FWHM, column density, and excitation temperature) are allowed to be free within reasonable and generous ranges unless otherwise stated. In most cases, we fit isotopologues separately rather so that their ratio can be constrained. After extensive testing, the best-fit solution for most species was found using 300 iterations of the Genetic algorithm to ensure that the global maximum is found, followed by 100 iterations of the Levenberg-Marquardt algorithm to ensure the best solution within that minimum (see Möller et al. 2017, for further details and references). For fits to high spectral-resolution H₂CO data, this was found to result in unreasonably large pixel-to-pixel variations, while using a 300 iterations of the Levenberg-Marquardt algorithm provided good results, so this latter approach was adopted.

For the fitting of the WideX spectra, each species was fitted separately, starting with the simpler and brighter species, except in cases where some lines are close enough to affect each other, in which case these lines were fitted together simultaneously (for example CH₃CN and HNCO). Once a good fit was achieved for a given molecule, the solution was held constant dring the fitting of subsequent molecules.

C.2 Detected transitions

Figure C.1 presents the results of XCLASS fitting to the IRAM 30 m H₂CO data for comparison with the results obtained with the merged NOEMA and IRAM 30 m data in Sect. 3.2. Figures C.2–C.8 presents the results of XCLASS fitting to the merged NOEMA and 30 m WideX spectra imaged with the SDI algorithm at the hot core, cold core 1, cold core 2 and outflow positions, as well as the NOEMA only WideX spectra imaged with the Clark algorithm at the cold core 1, cold core 2 and outflow positions.

Fig. C.1 XCLASS fitting for the IRAM 30 m H2CO data. The panels show: integrated intensity (between −87.5 and 2.2 km s−1) in the H2CO 30,3 −20,2 (top-left) and 32,2−22,1 (bottom-left) transitions; column density in p-H2CO (top-middle-left); excitation temperature (top-middle-right); peak velocity (bottom-middle-left); FWHM (bottom-middle-right); example spectra (black) and  XCLASS model fits (red) for H2CO 30,3 −20,2 (top-right) and 32,2−22,1 (bottom-right). The beam is shown in the top-left panel. The black contours in the maps show the NOEMA 1.4 mm continuum smoothed to a FWHM of 0.8′′, as used for the analysis of the combined H2CO data in Sect. 3.2, starting at 5σ and increasing to 45σ in 10σ steps. The spectra are extracted at the position covering the hot-core, which is indicated in the right-hand panels with a red square. A scale-bar is shown in the lower-left panel.

Fig. C.2 XCLASS fitting results for the SDI-imaged WideX data at the hot core position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey band on the data spectrum indicates 3σ (σ = 0.092 K).

Fig. C.3 XCLASS fitting results for the Clark-imaged WideX data at the outflow position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey bands on the data spectrum indicate 3σ (σ=0.045 K).

Fig. C.4 XCLASS fitting results for the SDI-imaged WideX data at the outflow position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey bands on the data spectrum indicate 3σ (σ=0.084 K).

Fig. C.5 XCLASS fitting results for the Clark-imaged WideX data at the cold core 1 position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey bands on the data spectrum indicate 3σ (σ=0.052 K).

Fig. C.6 XCLASS fitting results for the SDI-imaged WideX data at the cold core 1 position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey bands on the data spectrum indicate 3σ (σ=0.085 K).

Fig. C.7 XCLASS fitting results for the Clark-imaged WideX data at the cold core 2 position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey bands on the data spectrum indicate 3σ (σ=0.051 K).

Fig. C.8 XCLASS fitting results for the SDI-imaged WideX data at the cold core 2 position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey bands on the data spectrum indicate 3σ (σ=0.086 K).

References

All Tables

All Figures

Fig. 1 Overview of the W3 star formation region showing the larger-scale context of W3 IRS4. Left: infrared three-colour image of W, where red is Herschel Hi-GAL (Molinari et al. 2016) 250 μm, green is Hi-GAL 70 μm, and blue is MSX 21 μm (Price et al. 2001). The three brightest far-IR sources, W3 OH/H2O, W3 IRS4, and W3 IRS5, are indicated with black circles. The region for the middle panel is indicated with a black dashed box.Middle: three-colour image of W3 IRS4 and IRS5, where red is 14.9 GHz radio free-free continuum from Tieftrunk et al. (1997) to emphasise ionised gas, green is Spitzer IRAC 4.5 μm (Ruch et al. 2007), and blue is Large Binocular Telescope (LBT, Hill et al. 2006) Ks data from Bik et al. (2012). The locations of W3 IRS4 and IRS5 are indicated with black circles, with the white dashed box indicating the region covered in the left-hand panel. Right: three-colour image focusing on W3 IRS4, where red is the NOEMA 1.4 mm continuum presented in this paper, imaged with natural weighting, while green and red are as in the middle panel. The dashed white circle shows the full-width half-maximum of the primary beam for the NOEMA continuum observations.

In the text

Fig. 2 NOEMA 1.4 mm continuum intensity maps of the W3 IRS4 region, imaged with the Clark algorithm (see text for more details) and with natural (left and middle) and uniform (right) weighting. The beam sizes are 0.94′′ × 0.82′′ and 0.45′′ × 0.32′′ respectively, shown in the lower-left of each panel. The black contours indicate intensity iso-contours at 5, 10, 15, 20, 25 and 30σrms (σrms = 0.95 and 0.27 mJy beam−1 for the left and right panels, respectively), while white-dashed contours indicate −5σrms. The left and right panels are shown in flux density units while the middle panel has been converted to brightness temperature using the Rayleigh–Jeans approximation. The principle structures within the region are indicated in the middle panel: the hot core/MYSO W3 IRS4 with a solid black circle, the UCH II region W3 C with a dashed circle, two cold continuum sources with black crosses, and the filament with white lines. The coordinates are given in Table 2. Short-spacing information is not available for the continuum emission.

In the text

Fig. 3 Radio contribution to the NOEMA continuum emission around W3 C. Left: VLA 22.5 GHz radio continuum observations from Tieftrunk et al. (1997) rescaled to 220 GHz (see text for details). W3 C is in the centre, with weak emission from W3 IRS4, highlighted with a white dashed circle, to the NE and W3 D further away to the NW. Middle: NOEMA 220 GHz continuum emission. Right: NOEMA continuum map after the free–free contribution has been estimated from the VLA observations and removed (see text for details). All data were convolved to the same 0.8 × 0.8′′ beam. The black contours indicate intensity iso-contours at 5−30σ (σ = 0.4, 0.8 and 1.4 mJy beam−1 from left to right), increasing in 5σ increments, with the addition of a lower contour at 3σ for the VLA observations to show that the weak emission at the location of W3 IRS4 is significant. The white-dashed contours indicate −5σrms.

In the text

Fig. 4 Spatial and spectral distribution of gas and dust around the UCH II region W3 C. The colour maps show the integrated intensity in the combined NOEMA and 30 m data for 13CO 2−1 (top-left), C18O 2−1 (top-centre),c-C3H2 6−5 (middle-left), and H2CO 30,3 −20,2 (middle-centre), and the NOEMA-only continuum (top-right). The velocity integration ranges are given in the text. The grey contours indicate 5σ for the continuum map (σrms = 0.95 mJy beam−1). The beam is shown in the lower-left corner of each map. c-C3H2 and H2CO are shown using a square-root intensity scheme to emphasise the extended emission. The middle-right panel shows the fractional intensity of the continuum and gas emission as a function of radius from the centre of W3 C, averaged over θ. The centre was determined by hand and is indicated on the maps with a cyan +. The region containing emission related to the hot core W3 IRS4 and cold filament, indicated with a red polygon on the maps, is excluded from the radial analysis, while the maximum radius is indicated with a white dashed circle. The radius of peak emission for each tracer is indicated with a circle on the relevant maps. Gas spectra at three representative positions, indicated in the maps with red crosses which are numbered in the continuum panel, are shown in the lower panels.

In the text

Fig. 5 Comparison of near-IR Ks images from Bik et al. (2012, greyscale image, 0.9′′ resolution) with VLA 22.5 GHz radio continuum observations from Tieftrunk et al. (1997, orange contours) and free–free-subtracted NOEMA continuum at 220 GHz (magenta contours). The contour levels are the same as in Fig. 3. The cyan cross indicates the approximate geometrical centre of the UCH II region used as the centre for the radial analysis in Fig. 4, while the white circle indicates the source attributed by Bik et al. (2012) as powering W3 C and the red circle indicates the nearest source to the centre of the free–free emission, which is much redder.

In the text

Fig. 6 Maps of integrated intensity (moment 0, left, integrated between −55.3 and −29.8 km s−1), intensity-weighted velocity (moment 1, middle), and intensity-weighted FWHM (moment 2, right) for combined NOEMA and 30 m data for H2CO 30,3 −20,2, focusing on the filament. The grey contours indicate 5σ for the NOEMA 1.4 mm continuum with uniform weighting, as in the right-hand panel of Fig. 2. The cyan contours show the moment 0 of the H2CO 32,2 −22,1 transition, starting at 5σ (σ = 0.019 Jy beam−1 km s−1) and increasing to 25σ in 5σ steps. The mean velocity near the hot core is approximately −44.5 km s−1. The moment1 and 2 maps are calculated using channels above 4σ and only for pixels where the integrated intensity is also above 4σmom0, with grey used to indicate those regions that do not meet these requirements.

In the text

Fig. 7 XCLASS fitting for H2CO in the filament. The panels show: integrated intensity (between −55.3 and −29.8 km s−1) in the H2CO 30,3 −20,2 (top-left) and 32,2−22,1 (bottom-left) transitions; example spectra (black) and  XCLASS model fits (red) for H2CO 30,3 −20,2 (top-middle-left) and 32,2−22,1 (bottom-middle-left); column density in p-H2CO (top-middle-right); excitation temperature (top-right); peak velocity (bottom-middle-right); FWHM (bottom-right). The data for the two lines were resampled to a common spherical beam with a FWHM of 0.8′′, as shown in the top-right panel, and converted to the TMB scale. The grey contours in the maps show the continuum smoothed to the same resolution as the line data, starting at 5σ (σ = 0.81 mJy beam−1) and increasing to 45σ in 10σ steps. The spectra are extracted at the continuum peak (coordinates in Table 2), which is indicated in the left-hand panels with a red marker. A scale-bar is shown in the lower-left panel. The model results are only displayed for pixels where the integrated intensity for H2CO 30,3 −20,2 is above 5σmom0, with grey used to indicate those regions that do not meet this requirement. The white dashed contours in the map of excitation temperature indicate 3σmom0 for H2CO 32,2 −22,1.

In the text

Fig. 8 Calculated physical properties for the filament. The panels show: mass-per-pixel (left); H2 column density (middle); and volume number density (right) assuming a depth along the line-of-sight equal to the widest part of the filament (~8000 AU). The volumedensity is therefore likely an underestimate. The values are only calculated where the continuum is above 3σ (σ = 0.81 mJy beam−1). The blackcontours show the continuum smoothed to the same resolution as the line data, starting at 5σ and increasing to 45σ in 10σ steps. The beam and a scale-bar are shown in the left and middle panels respectively.

In the text

Fig. 9 LBT LUCI Ks (top-left, data from Bik et al. 2012) and Spitzer IRAC 4.5 μm (bottom-left) continuum maps, as well as integrated emission in the blue (top) and red (bottom) line-wings of combined NOEMA and 30 m maps of 13CO 2−1 (middle) and SO 6−5 (right). The NOEMA continuum emission with natural weighting is shown as green (left) or grey (middle and right) contours from 5 to 30σ in 5σ intervals (σ = 0.95 mJy beam−1). The velocity integration windows for the 13CO and SO maps are indicated in the right-hand panels. The purple arrows indicate the tentatively identified direction of the outflow, with dashed lines used to highlight features which may trace the outflow cavity walls (see text for details).

In the text

Fig. 10 XCLASS fitting for CH3CN in the hot core. The panels show: integrated intensity in the CH3CN 124 −114 transition (top-left); column density of CH3CN (top-middle); excitation temperature (top-right); peak velocity (bottom-middle); FWHM (bottom-right); example spectra (black) and  XCLASS model fits (bottom-left) for CH3 CN (red), CHCN (cyan), and HNCO (blue). As k = 0−3 of CH3 CN were notincluded in the fit, the resulting fit solution for those transitions is shown with a dashed line. The data were imaged with an intermediate weighting between uniform and natural in order to provide a good balance between sensitivity and resolution, resulting in a beam size of 0.67′′ × 0.56′′ (PA = −83.2°), as shown in the top-left panel with a scale-bar, and converted to the Tmb scale. The black contours in the maps show the continuum at the same resolution as the line data, starting at 5σ and increasing to 45σ in 10σ steps (σ = 0.48 mJy beam−1). The spectra are extracted at the continuum peak, which is indicated in the top-left panel with a red square. All maps are clipped to show only results where the integrated intensity of the CH3 CN 124 −114 is above 5σ (σ = 1.27 K kms). The purple arrow indicates the suggested direction of the outflow as in Fig. 9, while the green dashed line is perpendicular to this axis.

In the text

Fig. 11 Intensity-weighted velocity (moment 1, top, third, and fifth rows) and intensity-weighted FWHM (moment 2, second, fourth, and sixth rows) maps of key species in the high velocity-resolution data, along with their spectra at the continuum peak, as indicated with a red box on the moment 1 maps and a cyan box on the moment 2 maps. The lines are shown from left to right and top to bottom in order of increasing transition upper energy level (Eup). The NOEMA continuum emission with similar beam-size to the data are shown as grey contours from 5 to 45σ in 10σ intervals (σ = 0.43 mJy beam−1) in the maps, which are clipped at 5σ of the corresponding integrated intensity map. The purple arrows indicate the suggested direction of the outflow as in Fig. 9.

In the text

Fig. 12 Moment 1 (left) and position-velocity (PV) maps for molecular gas around W3 IRS4. From top to bottom, each row shows the moment 1 map and PV cuts along three different axes (middle-left, middle-right and right), for H2CO 32,2 −22,1, OCS 18-17, HNCO 101,9−91,8, CH3 CN 123 −113, and CH3CN 126 −116. The three PV cuts, taken along the outflow axis (1), rotated by 45° (2), and perpendicular to it (3), are shown on the moment 1 maps with grey lines numbered in red, where each cut is averaged over a width similar to the beam perpendicular to the cut. The 5σ contour of the data in each PV diagram is indicated with a white solid contour, while white dashed horizontal and vertical lines indicate the vLSR and zero offset from the continuum peak. The spatial and spectral resolution of the PV diagrams is indicated with a white cross in the lower-left corner of PV diagram 1 in each row. The black dashed lines indicate the results of the same cut through the moment 1 map, while the red dashed lines indicate what would be expected for Keplerian motion with a central mass of 10 and 20 M⊙. The black contours in the moment 1 maps show the continuum, as in Fig. 11.

In the text

Fig. 13 Temperature profile of H2CO as a function of offset from the MYSO and hot core W3 IRS4 along a slice which passes through the filament. Left: the position of the slice is indicated as a grey line on the temperature map from the  XCLASS fitting of H2CO (see Sect. 3.2.1 and Fig. 7). The temperature at each position is calculated as an average over the beam FWHM of the map (shown in the lower-left corner, 0.8′′), indicated with the grey dashed lines. This profile for H2CO (black), and that for CH3CN (red), are shown in the middle-left panel. The middle-right panel shows a zoom-in focusing on IRS4, with the profile shown in logarithmic units in the right-hand panel. The green dashed lines in both panels indicate three boundaries, between which we perform power-law fits to the temperature profile. The solid magenta line shows the fit over the whole region, while the cyan and blue lines show the fits to only the inner and outer regions. These latter two fits are shown as solid lines in the regions where they are fit and dashed lines where they are extrapolated to the same region as the magenta line. The results of these fits are given in Table 7. The vertical dashed black lines in the middle-left and right-hand panels indicate the beam size. The black contours in the maps show the continuum at the same resolution as the H2CO data, starting at 5σ and increasing to 45σ in 10σ steps (σ = 0.81 mJy beam−1), while the purple arrow in the middle-right panel indicates the suggested direction of the outflow as in Fig. 9.

In the text

	Fig. 14 WideX spectrum centred on the position of the hot core for both the combined NOEMA and 30 m data imaged with the SDI algorithm (cyan) and the NOEMA only data imaged with the Clark algorithm (black).
In the text

Fig. 15 Integrated intensity maps of the strongest unblended line of the principle detected molecular species, as well as the continuum (bottom-right), imaged with natural weighting and resampled to a circular 1.2′′ beam. The combined NOEMA and 30 m data, imaged using this SDI algorithm, were used for those species indicated with a • before the species label, while the NOEMA only data, imaged with the Clark algorithm, were used for all others. The white contour indicates 5σ on the continuum (σ = 1.56 mJy beam−1). In the continuum panel (bottom-right), the black dashed circle indicates the ring of continuum emission associated with W3 C, the black circle indicates the position of W3 IRS4 and the black crosses indicate the two continuum peaks associated with cold cores, as shown in Fig. 2. WideX spectra are fitted in (see Sect. 4.2) at the locations of the cold cores, IRS4 and a position in the outflow with a knot of continuum emission, indicated by the green circle.

In the text

Fig. 16 XCLASS fitting results for the Clark-imaged WideX data at the hot core position. Results for the outflow and two cold core positions, as well as SDI-imaging for the hot core, are shown in Figs. C.2–C.8. The data are shown in black and the combined fit in red, while the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen). The data is cut off below the peak intensities for 13CO and SO in order to improve the visibility of the fainter lines. Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey band on the data spectrum, which extends beyond the range of the data on the right and left, indicates 3σ (σ = 0.044 K).

In the text

Fig. 17 Abundance distribution of H2CO (left), DCN (middle), and CH3CN (right). The black contours in the maps show the continuum smoothed to the same resolution as the line data, starting at 5σ and increasing to 45σ in 10σ steps (σ = 0.81 mJy beam−1). The beam is shown in each panel and a scale-bar is shown in the middle panel. Abundances are only calculated where both the continuum emission and integrated intensity of the H2CO 303 −202 and DCN 3−2 lines are above 3σ and where CH3CN 124 −114 is above 5σ. The uncertainty in the abundances is conservatively estimated to be around a factor of two (see text for details).

In the text

Fig. 18 Abundance vs. temperature from the maps shown in Fig. 17 (dots) and the fits to the WideX spectra discussed in Sect. 4.2 (connected crosses). The different colours denote different species (see legend), using the same colours as in Fig. 16. The dashed lines indicate a rolling median of map points within ± 10 K, measured every 10 K as long as at least 10 points are included. The downward vertical arrows show the effect of increasing N[H2] at the outflow position (which is likely underestimated, leading to the molecular abundances being overestimated) such that the 12C/13C ratio would be 76. The uncertainty in the abundances is conservatively estimated to be around a factor of two (see text for details).

In the text

Fig. A.1 IRAM 30 m spectra at 11′′ resolution towards the hot core position of W3 IRS4, resampled to 3 km s−1 channels. The location of the main detected transitions within these frequency ranges are labelled and indicated with red dashed vertical lines (see Table C.1 for details). The location of SiO 5−4, which is not detected, is indicated in blue. The frequency axis has been shifted to the rest velocity using a source velocity of vLSR = −44.5 km s−1.

In the text

Fig. A.2 Overview of the IRAM 30 m data for selected lines, one per row. For each transition the panels show a spectrum at the source centre (left), as well as maps of the integrated intensity (moment 0, middle-left, between −87.5 and 2.2 km s−1), intensity-weighted velocity (moment 1, middle-right), and intensity-weighted line FWHM (moment 2, right). The position of the spectrum is indicated on the left-hand panel by a black square. The contours in the three maps indicate the NOEMA 1.4 mm continuum with uniform weighting, as in the right-hand panel of Fig. 2. The beam-size and a scale bar are given in the middle-left panel. For the middle-left and middle-right panels, values are only calculated where the spectrum has a peak intensity above 3σ, with grey used to indicate those regions that do not meet this requirement.

In the text

	Fig. A.2 continued.
In the text

Fig. B.1 Integrated intensity (0th moment, integrating between −55.3 and −29.8 km s−1) maps of H2CO 30,3 −20,2 imaged with the SDI algorithm (see text for more details) and uniform weighting for NOEMA only (left) andmerged NOEMA and 30 m data (right). The black-solid and white-dashed contours indicate 5σ and −5σ, where σ = 12.5 and 12.6 mJy beam−1 km s−1 for the NOEMA only and the merged data respectively. The black circle, dashed circle and crosses indicate the positions of W3 IRS4,W3 C and the two cold cores, as indicated in Fig. 2.

In the text

Fig. B.2 Integrated intensity (0th moment) maps of OCS 18−17 (left, between −55.3 and −29.8 km s−1) and C18O 2−1 (middle and right, between −81.8 and −0.8 km s−1) for data images with the Clark (top) and SDI (bottom) clean algorithms. All are imaged with uniform weighting. The OCS map is from only NOEMA data while the C18O is merged NOEMA and 30 m data. The C18O maps are without (middle) and with (right) primary beam correction, where in the latter case the data have been truncated at the FWHM of the primary beam. The black contours indicate 5σ, where σ = 13.8 (OCS, Clark), 12.6 (OCS, SDI), 19.0 (C18O before primary-beam correction, Clark), and 20.0 mJy beam−1 km s−1 (C18O before primary-beam correction, SDI). The black circle, dashed circle and crosses indicate the positions of W3 IRS4, W3 C and the two cold cores, as indicated in Fig. 2.

In the text

Fig. B.3 Comparison of the integrated intensity of SO 65−54 between −81.8 and −0.8 km s−1 with different imaging choices. The left and middle-left columns show results with only the NOEMA visibilities while the middle-right and right columns show the results with the 30 m short-spacing information included. The maps in the left and middle-right columns were imaged with the Clark algorithm while those in the middle-left and right columns were imaged with the SDI algorithm. The top row was imaged with uniform weighting, the middle row using a GILDAS-robust parameter of 3 to give an intermediate beam-size, while the bottom row using natural weighting. The intensity scale is chosen to show the extended emission, with the peak intensity for the hot core being between ~240 and ~790 K km s−1, depending on the method.

In the text

Fig. C.1 XCLASS fitting for the IRAM 30 m H2CO data. The panels show: integrated intensity (between −87.5 and 2.2 km s−1) in the H2CO 30,3 −20,2 (top-left) and 32,2−22,1 (bottom-left) transitions; column density in p-H2CO (top-middle-left); excitation temperature (top-middle-right); peak velocity (bottom-middle-left); FWHM (bottom-middle-right); example spectra (black) and  XCLASS model fits (red) for H2CO 30,3 −20,2 (top-right) and 32,2−22,1 (bottom-right). The beam is shown in the top-left panel. The black contours in the maps show the NOEMA 1.4 mm continuum smoothed to a FWHM of 0.8′′, as used for the analysis of the combined H2CO data in Sect. 3.2, starting at 5σ and increasing to 45σ in 10σ steps. The spectra are extracted at the position covering the hot-core, which is indicated in the right-hand panels with a red square. A scale-bar is shown in the lower-left panel.

In the text

Fig. C.2 XCLASS fitting results for the SDI-imaged WideX data at the hot core position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey band on the data spectrum indicates 3σ (σ = 0.092 K).

In the text

Fig. C.3 XCLASS fitting results for the Clark-imaged WideX data at the outflow position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey bands on the data spectrum indicate 3σ (σ=0.045 K).

In the text

Fig. C.4 XCLASS fitting results for the SDI-imaged WideX data at the outflow position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey bands on the data spectrum indicate 3σ (σ=0.084 K).

In the text

Fig. C.5 XCLASS fitting results for the Clark-imaged WideX data at the cold core 1 position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey bands on the data spectrum indicate 3σ (σ=0.052 K).

In the text

Fig. C.6 XCLASS fitting results for the SDI-imaged WideX data at the cold core 1 position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey bands on the data spectrum indicate 3σ (σ=0.085 K).

In the text

Fig. C.7 XCLASS fitting results for the Clark-imaged WideX data at the cold core 2 position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey bands on the data spectrum indicate 3σ (σ=0.051 K).

In the text

Fig. C.8 XCLASS fitting results for the SDI-imaged WideX data at the cold core 2 position. The data are shown in black and the combined fit in red, with the fits for individual species are shown in different colours below, in some cases scaled for greater visibility, with common colours used for different chemical families (blue for organics, yellow and orange for sulphur-bearing, red for nitrogen-bearing, purple for species bearing both nitrogen and oxygen, dark green for hydro-carbons). Vertical dashed lines are shown at the position of the brightest transition of each species to guide the eye. The grey bands on the data spectrum indicate 3σ (σ=0.086 K).

In the text