From clump to disc scales in W3 IRS4 A case study of the IRAM NOEMA large programme CORE

Context. High-mass star formation typically takes place in a crowded environment, with a higher likelihood of young forming stars affecting and being affected by their surroundings and neighbours, as well as links between different physical scales affecting the outcome. However, observational studies are often focused on either clump or disc scales exclusively. Aims. We explore the physical and chemical links between clump and disc scales in the high-mass star formation region W3IRS4, a region that contains a number of different evolutionary phases in the high-mass star formation process, as a case-study for what can be achieved as part of the IRAM NOrthern Extended Millimeter Array (NOEMA) large programme named CORE: “Fragmentation and disc formation in high-mass star formation”. Methods. We present 1.4mm continuum and molecular line observations with the IRAM NOEMA interferometer and 30m telescope, which together probe spatial scales from ∼ 0.3 − 20 ′′ (600 − 40000AU or 0.003 − 0.2pc at 2kpc, the distance to W3). As part of our analysis, we used XCLASS to constrain the temperature, column density, velocity, and line-width of the molecular emission lines. Results. The W3IRS4 region includes a cold ﬁlament and cold cores, a massive young stellar object (MYSO) embedded in a hot core, and a more evolved ultra-compact (UC)H II region, with some degree of interaction between all components of the region that affects their evolution. A large velocity gradient is seen in the ﬁlament, suggesting infall of material towards the hot core at a rate of 10 − 3 − 10 − 4 M ⊙ yr − 1 , while the swept up gas ring in the photodissociation region around the UCH II region may be squeezing the hot core from the other side. There are no clear indications of a disc around the MYSO down to the resolution of the observations (600AU). A total of 21 molecules are detected, with the abundances and abundance ratios indicating that many molecules were formed in the ice mantles of dust grains at cooler temperatures, below the freeze-out temperature of CO ( . 35K). This contrasts with the current bulk temperature of ∼ 50K, which was obtained from H 2 CO. Conclusions. CORE observations allow us to comprehensively link the different structures in the W3IRS4 region for the ﬁrst time. Our results argue that the dynamics and environment around the MYSO W3IRS4 have a signiﬁcant impact on its evolution. This context would be missing if only high resolution or continuum observations were available.

M o t t r a m , J. C., B e u t h er, H ., Ah m a di, A., Kla a s s e n , P. D., B elt r á n, M . T., C s e n g e ri, T., F e n g , S., Gi e s er, C., H e n ni n g, Th., Joh n s t o n, K. G., Kuip er, R., Le u ri ni, S., Linz, H ., Lo n g m o r e , S. N., Lu m s d e n , S., M a u d , L. T., M o s c a d elli, L., P al a u, A., P e t e r s , T., P u d ri tz, R. E., R a g a n , S. Thi s v e r sio n is b ei n g m a d e a v ail a bl e in a c c o r d a n c e wit h p u blis h e r p olici e s. S e e h t t p://o r c a . cf. a c. u k/ p olici e s. h t ml fo r u s a g e p olici e s. Co py ri g h t a n d m o r al ri g h t s fo r p u blic a tio n s m a d e a v ail a bl e in ORCA a r e r e t ai n e d by t h e c o py ri g h t h ol d e r s .

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 mmwavelength 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" ) is GMA (a) 150 10 10 6 −10 7 6 GMC (b) 50 100 10 5 1,2,4,6 Cloud 10 500 10 4 1,2,3,4,6 Clump 1 5000 10 3 1,2,3,4,6,8 Filament 1−10 × 10 4 1−30 × ℓ 5,6 0.05−0.25 Cores 0.1 10 4 −10 6 10−10 3 1,2,3,5 Embedded 5−500 × 10 −5 10 6 −10 9 0.1−100 7 Discs (10-1000 AU) Notes. 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 2 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 2 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  (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/H 2 O, 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) K s 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.
(for example Tieftrunk et al. 1997Tieftrunk et al. , 1998Wang 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 2 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 4 L ⊙ from SED fitting 1 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 A&A 636, A118 (2020)  Notes. (a) Taken as the peak of the continuum emission after resampling to a circular beam of 0.8 ′′ , as used in much of the analysis in Sect. 3. (b) Approximate geometric centre of the ring of NOEMA continuum emission.
( 40 ′′ , 80 000 AU or 0.4 pc) to the NE and N in VLA NH 3 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 with surrounding hot core, and an UCH II region) next to each other in one region, making it both an efficient and interesting 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 shortspacing information from IRAM 30 m observations for 12 CO and 13 CO, so could not explore the variation of temperature and abundance in the region, and their limited velocity and spatial resolution (1.2 km s −1 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 where discs 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. 2017Meyer et al. , 2018Kö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 higher spectral resolution, while also including shortspacing 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 observations and 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.

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.

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 2 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 −1 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 process the signals of six antennas. During the observations, the rest-frame of W3 IRS4 was taken to be −42.8 km s −1 , based on the IRAM 30 m 12 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 the 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 2 . For line data, the continuum was subtracted using linefree channels in the same band, while a combined continuum

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 −1 at 217 GHz and 0.25 km s −1 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 −1 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.

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 highand low-resolution bands (see Table 4) are imaged with a velocity resolution of 0.5 and 3.0 km s −1 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.

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 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 freefree 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 ( 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 4 cm −3 , considerably less than the 10 6 cm −3 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 0 ) can also be estimated from the radio continuum flux (again following Sánchez-Monge et al. 2013) as: and is calculated to be 2 × 10 47 photons s −1 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 44 photons s −1 . 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 4 L ⊙ . As the combined bolometric luminosity for W3 IRS4 and W3 C is ∼4.5 × 10 4 L ⊙ (see Sect. 1), this would leave ∼2.3 × 10 4 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 4 L ⊙ , leaving only ∼9 × 10 3 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 . 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 13 CO 2−1 (top-left), C 18 O 2−1 (top-centre), c-C 3 H 2 6−5 (middle-left), and H 2 CO 3 0,3 −2 0,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-C 3 H 2 and H 2 CO 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.
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. 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, 13 CO 2−1, C 18 O 2−1, c-C 3 H 2 6−5, and H 2 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 −1 ) is from −81.8 to −0.8 km s −1 for 13 CO and C 18 O, and −48.8 to −28.8 km s −1 for c-C 3 H 2 in order to avoid contamination at the position of IRS4 from 33 SO. For H 2 CO, which was observed at high spectralresolution (0.5 km s −1 ), −55.3 to −29.8 km s −1 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 3 H 2 .

Molecular line emission
A118, page 7 of 44 A&A 636, A118 (2020) These gas tracers, along with SO 6 5 −5 4 and possibly HC 3 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 3 H 2 peaks closest to the UCH II region at 3.6 ′′ (7200 AU), followed by 13 CO, C 18 O, and H 2 CO at 4.9, 5.3, and 6.1 ′′ (9800, 10 600, and 12 200 AU) respectively. Only 13 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 2 CO spectra showing signs of two velocity components, one around −44.5 km s −1 , the lsr of IRS4, and another around −42.5 km s −1 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 −1 , offset by ∼20 km s −1 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 13 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 3 H 2 peaks closer to 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). 13 CO and C 18 O then peak next, followed by H 2 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 2 CO slightly closer to the ionising source than C 18 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 18 O has a lower relative intensity near the centre of the UCH II region than 13 CO (see lower-right panel of Fig. 4). This could be due to the peak 13 CO intensity being decreased at the edge of the shell due to higher optical depth with respect to the centre and C 18 O. However, it is also consistent with isotopeselective photodissociation by FUV photons, where 13 CO is able to more effectively self-shield than C 18 O at smaller radii due to its higher abundance (see for example Glassgold et al. 1985;Röllig & Ossenkopf 2013;Shimajiri et al. 2014).

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 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. spectroscopy of W3 Main using LUCI 1 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. , intensity-weighted velocity (moment 1, middle), and intensity-weighted FWHM (moment 2, right) for combined NOEMA and 30 m data for H 2 CO 3 0,3 −2 0,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 H 2 CO 3 2,2 −2 2,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 moment 1 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.

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 3 observations (Tieftrunk et al. 1998). This structure is also detected in a number of the molecular gas tracers, including H 2 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 2 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 −1 which starts in 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 largescale velocity shift is even seen the 30 m data for H 2 CO, SO, C 18 O, and, to a lesser extent, 13 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.

XCLASS fitting of H 2 CO
The two H 2 CO transitions included in the high spectralresolution 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, columndensity of p−H 2 CO, peak velocity, and FWHM for each pixel, the maps of which are shown in Fig. 7. Assuming that H 2 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 2 CO] = 5.8 × 10 13 cm −2 ,v = −45.3 km s −1 , and FWHM = 2.9 km s −1 . 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 03 −2 02 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 3 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 2 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 temperature of 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 2 O (∼100 K, for example Fraser et al. 2001). This affects the chemistry in both the gas and solid phases, by in the H 2 CO 3 0,3 −2 0,2 (top-left) and 3 2,2 −2 2,1 (bottom-left) transitions; example spectra (black) and XCLASS model fits (red) for H 2 CO 3 0,3 −2 0,2 (top-middle-left) and 3 2,2 −2 2,1 (bottom-middle-left); column density in p-H 2 CO (top-middle-right); excitation temperature (top-right); peak velocity (bottom-middleright); 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 T MB 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 H 2 CO 3 0,3 −2 0,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 H 2 CO 3 2,2 −2 2,1 . Notes. (a) Calculated over a circular region with a radius 1 ′′ , centred at the coordinates given in Table 2.
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).

Mass
If we assume that the dust is optically thin and in thermal equilibrium with the gas, then the gas temperature derived from H 2 CO and the dust continuum emission can be used to calculate maps of N[H 2 ] and mass following the prescription by Hildebrand (1983) in the forms presented by Schuller et al. (2009): and 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 , and κ ν = 0.9 cm 2 g −1 from Ossenkopf & Henning (1994) for a wavelength of 1.4 mm assuming coagulation for 10 5 yr at n = 10 6 cm −3 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 2 column density were calculated using the temperature map derived from H 2 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 −4 M ⊙ pixel −1 and 2.6 × 10 22 cm −2 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 Cold core 1 5.4 7.8 6.5 1. properties). The median column density is 2.8 × 10 23 cm −2 with a maximum value of 1.8 × 10 24 cm −2 . These values for the masses and H 2 column densities are likely underestimated due to the large amount of missing continuum flux towards W3 IRS4 (∼78%, Beuther et al. 2018).

Volume density
The third panel in Fig. 8 shows an estimate for the volume number density, n, calculated from N[H 2 ] 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 2 ] is likely underestimated due to the missing flux in the continuum map. The median value of 2 × 10 6 cm −3 is well above the effective critical density for H 2 CO under the conditions in the region (n eff ∼10 5 cm −3 for 10 K, ∼10 4 cm −3 for 100 K, Evans 1999;Gerner et al. 2014;Shirley 2015), so the assumption of LTE for the XCLASS modelling, and therefore that the gas and dust are in thermal equilibrium, is sound.

Stability
For regions within the filament, the thermal Jeans mass can be measured using: 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: Using the median temperature of 50 K, the filament has a value of M l,crit ∼114 M ⊙ pc −1 . 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 −1 . 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 2 CO is 2.9 km s −1 , while the thermal contribution to the FWHM at the median temperature of 50 K is 0.28 km s −1 , 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 A118, page 11 of 44 A&A 636, A118 (2020) filament (see for example Fiege & Pudritz 2000;Kirk et al. 2015). If the thermal sound-speed is replaced with the observed FWHM/ √ 8 ln (2), thus assuming that all the non-thermal contribution to the line-width provides turbulent support against gravity, then Eq. (5) becomes: and Eq. (6) becomes: Using the median FWHM for the filament from H 2 CO of 2.8 km s −1 , M l,crit,nonthermal = 658 M ⊙ pc −1 , so the filament could be stable against radial collapse if supported by nonthermal 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: which is of the order of 10 4 yr for the three regions considered (see Table 6).

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 −1 , 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 5 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 3 observations of 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 −1 pc −1 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 bẏ Treating the filament as a solid body results in mass flow rates along the axis of the filament of 3.3 × 10 −4 M ⊙ yr −1 and 1.0 × 10 −3 M ⊙ yr −1 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 −4 M ⊙ yr −1 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 higherdensity 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 2 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 4 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.

Comparison with single-dish-only results
It is helpful to quantify how these results from resolved observations compare to what one would obtain with only singledish observations. A fit to the IRAM 30 m data for H 2 CO (see Fig. C.1) gives T ex = 67.5 K, N[p-H 2 CO] = 9.6 × 10 13 cm −2 , v = −44.2 km s −1 , and FWHM = 4.2 km s −1 for the central position. The values of T ex , N[p-H 2 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 −1 ) 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 beamaveraged H 2 column density similarly increases from 6.2 × 10 22 to 8.6 × 10 22 cm −2 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 2 column density.

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 and submm sources in W3. The combined bolometric luminosity for the W3 IRS4 region from single-dish measurements is ∼4.5 × 10 4 L ⊙ (see Sect. 1) and a contribution of ∼2.2 × 10 4 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. 3 and 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 H II 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 12 CO, 13 CO and C 18 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 2 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 linewings of 13 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 −1 from the v LSR of IRS4, which is approximately −44.5 km s −1 from H 2 CO (see Table 5), CH 3 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 the cavity (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 13 CO primarily traces transverse motions perpendicular to the direction of the outflow. The region of higher emission in the map of red-shifted 13 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).

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 3 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 3 CN and CH 3 13 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 3 CN 12 6 −11 6 . Since CH 3 13 CN may be present but is not clearly detected, we assume a ratio of CH 3 13 CN to CH 3 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 2 O). Following Ahmadi et al. (2018), we only fit the K = 4-6 transitions of CH 3 CN as this gives more consistent results.
Unlike H 2 CO, the peak CH 3 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,   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 13 CO 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 13 CO 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). 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 2 O (see for example Ahmadi et al. 2018). Towards IRS4, some transitions show complicated or unclear velocity structures (for example CH 3 OH, HC 3 N, CH 3 CN 12 2 −11 2 ), while in others there is a shift in the axis of the velocity gradient from more NW-SE at lower energies (for example H 2 CO) to N-S (for example HNCO, CH 3 CN 12 3 −11 3 , and CH 3 CN 12 6 −11 6 ). 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 2 CO 3 2,2 −2 2,1 , OCS 18-17, HNCO 10 1,9 −9 1,8 , CH 3 CN 12 3 −11 3 , and CH 3 CN 12 6 −11 6 . 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 −1 ) 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 13 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 3 CN 12 3 −11 3 , 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)   3 CN (cyan), and HNCO (blue). As k = 0−3 of CH 3 CN were not included 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 T mb 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 CH 3 CN 12 4 −11 4 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. 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 3 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.

Temperature structure
Returning to the XCLASS results, the median temperature obtained from fitting CH 3 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 2 O and other solidphase 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 3 CN (see Fig. 10) and H 2 CO over the small region around IRS4 where both are well detected (middle-left panel). The uncertainties increase in the H 2 CO temperature map for values well above the upper energy level of the higher excited transition (68 K). In addition, CH 3 CN does not show the high-temperature spot to the SE of the hot core that is seen in the H 2 CO map, though this is where the CH 3 CN intensity  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 (E up ). 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 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 H 2 CO 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.
becomes lower so the temperature from CH 3 CN becomes more uncertain due to lower sensitivity. The hot spot in H 2 CO could also be due to chemical effects or enhancement in the outflow, and CH 3 CN likely traces denser conditions, on average, along the line of sight than H 2 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 2 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 2 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.

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 2 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   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 of a 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 2 CCO, t−HCOOH, NH 2 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 −1 channels of the WideX data.

Spatial intensity distributions
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 18 O, and c−C 3 H 2 , 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 18 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 13 CO near IRS4 is offset slightly towards the SE, similar to H 2 CO, so the 13 CO emission near IRS4 may be outflow related, which would explain why it doesn't show the same distribution as C 18 O. c−C 3 H 2 , 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 18 O, c−C 3 H 2 does not show much enhancement at the position of hot core. The absence of an enhancement of c−C 3 H 2 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 33 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 3 H 2 is often detected in UV-illuminated gas in photodissociation regions (PDRs) and H II regions, together with C 2 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 3 H 2 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 3 H 2 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 (bottomright), 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.   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 13 CO 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).

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 2 ], using the temperature derived from H 2 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 2 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 13 CO, C 18 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.
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 2 , 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 temperature and 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 18 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 12 CO, so in order to be able to explore the 12 C/ 13 C ratio using the observed 13 CO abundances these are ratioed to the 12 C abundance for the solar system from Asplund et al. (2009). Since 12 C is less abundant than 16 O in the solar abundance distribution, we make the simple assumptions that the availability of 12 C sets the abundance for 12 CO, that 12 CO is the main source of 12 C in the W3 IRS4 region, and that the cause of the change in 12 C/ 13 C with Galactocentric radius is due to variation of 13 C rather than 12 C.
The global 12 C/ 13 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 12 C to 13 CO and H 2 CO to H 13 2 CO for the hot core give a ratio of <40. The determination of N[H 2 ] for the hot core cancels out for the ratio of H 2 CO to H 13 2 CO, since the same value comes into the calculation of both abundances, but does not for the ratio of solar 12 C to 13 CO since it only enters the calculation for X[ 13 CO]. However, the resonable agreement between the two measured ratios suggests that the value of N[H 2 ] used for the hot core position is relatively good and therefore that the derived abundances are reasonable.   Asplund et al. (2009). (b) The solar 12 C abundance with respect to H 2 is taken to be 1.33 × 10 −4 (Asplund et al. 2009).
The lower 12 C/ 13 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 3 13 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 12 C/ 13 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 3 OH and CH 3 CN compared to an expected ratio of ∼60-70, Palau et al. 2017) and G31.41+0.31 ( 10 in CH 3 CN and CH 3 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, 13 C can be enhanced in CO at low temperatures due to the reaction 13 C + + 12 CO −→ 13 CO + 12 C + + 35K, and depleted in species which form through reactions with C + (for example Langer et al. 1984). Low-temperature gas-phase enhancement of 13 CO relative to 12 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 2 CO can be formed at warmer temperatures in the gas-phase via or at cold temperatures on grain surfaces via sequential hydrogenation of CO a process that also leads to the formation of CH 3 OH (for example Walsh et al. 2014): The similar enhancement of 13 C in CO and H 2 CO, with the ratio in CO slightly lower than in H 2 CO as predicted by Langer et al. (1984), is consistent with a cold temperature origin for H 2 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 3 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 3 OH in W3 IRS4. The ratio of A-CH 3 OH ν t = 1 to A-CH 3 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 13 CO/C 18 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 18 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 12 C to 13 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 12 C to 13 CO abundance ratio at the outflow position (∼11) could be due to N[H 2 ] 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 13 CO/C 18 O ratio than is observed. If the 12 C/ 13 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 2 ] 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 2 are similar to those obtained by Esplugues et al. (2013)  various regions around Orion KL. As in that study, our value for 34 SO 2 / 33 SO 2 = 5.8 is close to the solar abundance ratio of 5.6. The extremely low value for SO/ 33 SO is likely due to optical depth effects in the main isotopologue, but it seems less likely that SO 2 is suffering from such issues. It is possible that the isotopologues of SO 2 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, H 13 2 CO, HCOOH, and CH 3 CN towards IRAS 20126+4104. However, the values for H 2 CCO, CH 3 OH, CH 3 OCHO, and CH 3 COCH 3 are approximately an order of magnitude higher in IRAS 20126 (Palau et al. 2017), and both CH 3 CN and CH 3 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.

Spatial abundance distributions
The maps of N[H 2 ] and N[p-H 2 CO] can be used to produce a map of the abundance of H 2 CO by assuming a value for the ortho-to-para ratio for H 2 CO. If H 2 CO is produced in the gasphase, 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 2 CO/H 13 2 CO 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 2 CO] varies between 2.4 × 10 −9 and 1.1 × 10 −7 around the hot core IRS4 with a median of 1.0 × 10 −8 . The median abundance over the rest of the filament is 7.4 × 10 −10 , though these are upper limits since the H 2 CO data contain the large-scale information but the continuum (used for N[H 2 ]) 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 2 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 2 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 2 CO of 5−7 × 10 −10 , 5−7 × 10 −10 , 2−9 × 10 −9 , and 0.7−5 × 10 −9 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 ν 4 scaling, and use this to calculate N[H 2 ] = 6.2 × 10 22 cm −2 using the temperature from the fit to the IRAM 30 m H 2 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 2 CO] = 4.1 × 10 −9 , 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 2 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 2 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 2 CO. The values are broadly within the range predicted by the models of   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[H 2 ] at the outflow position (which is likely underestimated, leading to the molecular abundances being overestimated) such that the 12 C/ 13 C ratio would be 76. The uncertainty in the abundances is conservatively estimated to be around a factor of two (see text for details). 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 3 CN varies in the map between 9.5 × 10 −10 and 9.6 × 10 −9 with a median of 2.7 × 10 −9 , 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 3 OH (∼10 −2 ) 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). 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 2 ] (see Sect. 4.2 and the downward arrows in Fig. 17), most species show an increasing abundance with increasing temperature. The abundance of H 2 CO rises more quickly below ∼90−100 K than above it, and CH 3 CN, CH 3 OH, HNCO, c−C 3 H 2 , and SO all show 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).

Variation of abundance with temperature
In the case of H 2 CO, CH 3 CN, CH 3 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-thermal desorption 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 2 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 3 H 2 and SO may be 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 13 CO and C 18 O depend on the correction needed for N[H 2 ] 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 3 N show modest variation with temperature at best, suggesting that their abundances are maintained through gasphase 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 A118, page 25 of 44 A&A 636, A118 (2020)  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.

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 combination of 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 largerscale 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.

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 C likely 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 4 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 −4 −10 −3 M ⊙ yr −1 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 2 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 5 −10 6 yr. Since the timescales for highmass prestellar cores has been observationally constrained to be very short (10 4 −10 5 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 4 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.

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 2 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.

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 . 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 13 CO, C 18 O, H 2 CO, c-C 3 H 2 , and SO. The radial intensity distribution of these species suggest that isotopeselective 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 −1 , median temperature of ∼50 K, median volume density of 2 × 10 6 cm −3 , 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 −4 −10 −3 M ⊙ yr −1 (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 3 N, which therefore likely have a gas-phase origin. The most chemically complex species detected is acetone (CH 3 COCH 3 ), while the spatial distribution of c−C 3 H 2 stands out from the other species, likely tracing UV-illuminated gas, as suggested for PDRs and outflows from low-mass Class 0 protostars. A118, page 27 of 44 A&A 636, A118 (2020) 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 2 CO) will provide fruitful grounds for further studies when expanded to the to full CORE sample.   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 v LSR = −44.5 km s −1 .
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).

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 (T * A ) 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 4 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 −1 ) 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 −1 ). 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 12 CO and 13 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-positionvelocity 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.1 shows 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

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 −1 . 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. 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.

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 18 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 −5 4 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 −5 4 transition  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 2 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.

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 Notes. (a) Taken from the Cologne Database for Molecular Spectroscopy (CDMS, Müller et al. 2001Müller et al. , 2005. (b) Not detected towards W3 IRS4. (c) Blend of two lines with different quantum numbers but the same E up at the same frequency. The highest A ul is given.
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 2 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 3 CN and HNCO). Once a good fit was achieved for a given molecule, the solution was held constant dring the fitting of subsequent molecules.  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).

C.2. Detected transitions
A118  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).
A118 . 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).
A118  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).
A118 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).