A case study of the IRAM NOEMA large program CORE

Context. The fragmentation mode of high-mass molecular clumps and the properties of the central rotating structures surrounding the most luminous objects have yet to be comprehensively characterised. Aims. We study the fragmentation and kinematics of the high-mass star-forming region W3(H2O), as part of the IRAM NOrthern Extended Millimeter Array (NOEMA) large programme CORE. Methods. Using the IRAM NOEMA and the IRAM 30m telescope, the CORE survey has obtained high-resolution observations of 20 well-known highly luminous star-forming regions in the 1.37 mm wavelength regime in both line and dust continuum emission. Results. We present the spectral line set-up of the CORE survey and a case study for W3(H2O). At ∼0′′ .35 (700 AU at 2.0 kpc) resolution, the W3(H2O) clump fragments into two cores (west and east), separated by ∼2300 AU. Velocity shifts of a few km s−1 are observed in the dense-gas tracer, CH3CN, across both cores, consistent with rotation and perpendicular to the directions of two bipolar outflows, one emanating from each core. The kinematics of the rotating structure about W3(H2O) W shows signs of differential rotation of material, possibly in a disk-like object. The observed rotational signature around W3(H2O) E may be due to a disk-like object, an unresolved binary (or multiple) system, or a combination of both. We fit the emission of CH3CN (12K − 11K) K = 4 − 6 and derive a gas temperature map with a median temperature of ∼165 K across W3(H2O). We create a Toomre Q map to study the stability of the rotating structures against gravitational instability. The rotating structures appear to be Toomre unstable close to their outer boundaries, with a possibility of further fragmentation in the differentially rotating core, W3(H2O) W. Rapid cooling in the Toomre unstable regions supports the fragmentation scenario. Conclusions. Combining millimetre dust continuum and spectral line data toward the famous high-mass star-forming region W3(H2O), we identify core fragmentation on large scales, and indications for possible disk fragmentation on smaller spatial scales.


Introduction
Fundamental questions pertaining to the fragmentation of highmass clumps and the accretion processes that result in the birth of the most massive stars (M 8 M ) still remain unanswered. This is in part due to the clustered nature of star formation and the typically large distances involved. For a long time the existence of high-mass stars had been puzzling as it was thought that the expected intense radiation pressure would prevent the accretion of enough material onto the protostar (e.g. Kahn 1974;Wolfire & Cassinelli 1987). More recently, two-and three-dimensional (magneto)hydrodynamical simulations of collapsing cores have validated the need for accretion disks in the formation of very massive stars, analogous to low-mass star formation (e.g. Yorke & Sonnhalter 2002;Krumholz et al. 2009;Peters et al. 2010;Kuiper et al. 2010Kuiper et al. , 2011Kuiper & Yorke 2013;Klassen et al. 2016). Furthermore, different fragmentation processes can contribute to the final stellar mass distribution within a single region, including fragmentation from clouds down to core scales (e.g. Bontemps et al. 2010;Palau et al. 2013, Based on observations from an IRAM large program. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). 2015; Beuther et al. 2018; see review by Motte et al. 2017) and disk fragmentation at smaller spatial scales (e.g. Matsumoto & Hanawa 2003; see review by Kratter & Lodato 2016).
In the disk-mediated accretion scenario, the non-isotropic treatment of the radiation field reduces the effect of radiation pressure in the radial direction, such that radiation can escape through the poles along the disk rotation axis, while the disk is shielded due to the high densities. Observationally, the existence of these disks is expected due to ubiquitous observations of collimated outflows (e.g. Beuther et al. 2002;Fallscheer et al. 2009;Leurini et al. 2011;Frank et al. 2014;Maud et al. 2015), which have also been predicted by theoretical models (e.g. Pudritz et al. 2007). Although some accretion disks in differential, Keplerianlike rotation about B-type (proto)stars have been found in recent years (e.g. Carrasco-González et al. 2012;Sánchez-Monge et al. 2013;Beltrán et al. 2014; see reviews by Cesaroni et al. 2007, andde Wit 2016), the existence of such rotating structures around the most massive, O-type protostars is still elusive; only a few cases have been reported to date (Johnston et al. 2015;Ilee et al. 2016;Cesaroni et al. 2017).
As higher resolution observations are becoming more accessible, thus allowing structures to be resolved on scales <1000 AU, it is important to determine whether disks around intermediate to high-mass stars (OB-type) are ubiquitous, and if so, to characterise their properties. What is the typical extent of these disks? Are they in differential rotation about a centrally dominating protostar, similar to their low-mass counterparts, and if so, over what range of radii? Is there any scale where a core stops fragmenting? 1 At what scales do we see the fragmentation of disks? Are close binary/multiple systems an outcome of disk fragmentation (as suggested by e.g. Meyer et al. (2018))? If so, stability analyses of these high-mass rotating cores and disks are needed to shed light on fragmentation at disk scales. These questions can only be answered with a statistical approach for a large sample of high-mass star-forming regions.
We have undertaken a large programme at IRAM, called CORE (Beuther et al. 2018), making use of the IRAM NOrthern Extended Millimeter Array (NOEMA, formerly Plateau de Bure Interferometer) at 1.37 mm in both line and continuum emission to study the early phases of star formation for a sample of 20 highly luminous (L > 10 4 L ) star-forming regions at high angular resolution (∼0.4 ) to analyse their fragmentation and characterise the properties of possible rotating structures. Additionally, observations with the IRAM 30m telescope are included to complement the interferometric data, allowing us to understand the role of the environment by studying high-mass star formation at scales larger than those covered by the interferometer. Observations in the 1.3 mm wavelength regime of the CORE project began in June 2014 and finished in January 2017, consisting of a total of more than 400 hours of observations with NOEMA. The sample selection criteria and initial results from the observed level of fragmentation in the full sample are presented in Beuther et al. (2018), and details of the 30m observations and the merging of single-dish observations with interferometric observations can be found in Mottram et al. (in prep.). In this work, we describe our spectral set-up and present a case study of one of the most promising star-forming clouds in our sample, W3(H 2 O).
W3(H 2 O), also known as the Turner-Welch object, resides in the W3 high-mass star-forming region and was initially identified through observations of the dense-gas tracer HCN at 88.6 GHz (Turner & Welch 1984). It is located ∼0.05 pc (5 ) east of the well-known ultra compact H ii region (UCHII) W3(OH). The name W3(H 2 O) stems from the existence of water masers in the vicinity of the source (Dreher & Welch 1981), allowing for an accurate distance measurement of 2.0 kpc for this region (Hachisuka et al. 2006;cf. Xu et al. 2006). The relative proper motions of these masers are further explained by an outflow model oriented in the east-west direction (Hachisuka et al. 2006). A continuum source elongated in the east-west direction and spanning the same extent as the water maser outflow has been observed in subarcsecond VLA observations in the radio regime with a spectral index of -0.6, providing evidence for synchrotron emission Wilner et al. 1999). This source of synchrotron emission has been characterised by a jetlike model due to its morphology, and the point symmetry of its wiggly bent structure about the centre hints at the possibility of jet precession. Moreover, Shchekinov & Sobolev (2004) attribute this radio emission to a circumstellar jet or wind ionised by the embedded (proto)star at this position. Additional radio continuum sources have been detected in the vicinity of the synchrotron jet, the closest of which is to the west of the elongated structure and has a spectral index of 0.9 (Wilner et al. 1999;Chen 1 Here a core is defined as a gravitationally bound region that forms single or multiple stars, following Williams et al. 2000.  Wyrowski et al. (1999) in the 1.36 mm band allowed for the detection of three continuum peaks in thermal dust emission, one of which peaks on the position of the water maser outflow and synchrotron jet, and another on the position of the radio continuum source with positive spectral index, confirming the existence of a second source at this position. The detection of two bipolar molecular (CO) outflows further supports the protobinary scenario, suggesting that W3(H 2 O) may be harbouring at least two rotating structures ). The two cores within W3(H 2 O) have individual luminosities of the order of 2 × 10 4 L , suggesting two 15 M stars of spectral type B0 2 .
In this paper, we study the fragmentation properties of W3(H 2 O) and the kinematics of the rotating structures within it. We use this source as a test-bed for an expanded study in a forthcoming paper which will focus on the kinematic properties of a larger sample within our survey. The structure of the paper is as follows. Section 2 presents our spectral line set-up within the CORE survey with the details of our observations and data reduction for W3(H 2 O). The observational results are described in Section 3. The kinematics, temperature, and stability analysis of W3(H 2 O) is presented in Section 4. The main findings are summarised in Section 5.

NOEMA observations
Observations of W3(H 2 O) at 1.37 mm were made between October 2014 and March 2016 in the A-, B-, and D-array configurations of NOEMA in track-sharing mode with W3 IRS4. The compact D-array observations were made with six antennas, while seven antennas were used for the more extended Aand B-array observations. Baselines in the range of 19 − 760 m were covered, therefore the NOEMA observations are not sensitive to structures larger than 12 (0.1 pc) at 220 GHz. On-source observations were taken in roughly 20-minute increments distributed over an observing run and interleaved with observations of various calibration sources. The phase centre for the observations of W3(H 2 O) is α(J2000) = 02 h 27 m 03 s .87, δ(J2000) = 61 • 52 24 . 5. A summary of the observations can be found in Table 1.
The full CORE sample of 20 regions has been observed simultaneously with a narrow-and a wide-band correlator. The wide-band correlator, WideX, has four units, each with 1.8 GHz bandwidth, covering two overlapping ranges in frequency in both horizontal and vertical polarisations (H and V) with a fixed spectral resolution of 1.95 MHz (∼2.7 km s −1 at 219 GHz). 2 The luminosity and spectral type calculations are described in detail in Section 4.3.
The full coverage of the WideX correlator is shown in Fig. 1 with bright lines marked. The narrow-band correlator has eight units, each with 80 MHz bandwidth and a spectral resolution of 0.312 MHz (∼0.43 km s −1 ), placed in the 1.37 mm wavelength regime. The frequency coverage of the correlator bands are listed in Table 2. The narrow-band correlator can only process the signal from six antennas; therefore, in cases where the sources were observed with more than six antennas, the correlator automatically accepts the signal from the antennas that yield the best uvcoverage. Important lines covered by the narrow-band receiver are listed in Table 3 and presented in Fig. 2 for the pixel at the phase centre toward W3(H 2 O).
Data reduction and imaging were performed with the clic and mapping program of the gildas 3 software package developed by the IRAM and Observatoire de Grenoble. The continuum was extracted by identifying line-free channels in the range 217 078.6−220 859.5 MHz covered by all four spectral units of the WideX correlator. As we are interested in achieving the highest possible angular resolution, we CLEANed the cubes using the CLARK algorithm (Clark 1980) with a uniform weighting (robust parameter of 0.1) 4 yielding a synthesised beam size of 0 . 43 × 0 . 32, PA=86 • , and an rms noise of 3.2 mJy beam −1 for the continuum emission using the combined set of observations in the A-, B-, and D-array configurations (hereafter ABD). We also imaged the data from just the A-and B-array configurations (hereafter AB), as well as the A-array only, for which the synthesised beam sizes and rms noise values are summarised in Table 4.
Continuum subtraction for the lines was performed in the uvplane by subtracting the emission in the line-free channels in the spectral unit in which the line was observed. Due to line contamination in spectral unit L03, we used the continuum from spectral unit L02 to remove the continuum from spectral unit L03, under the assumption that there is no significant spectral slope between the two adjacent spectral windows. For the WideX continuum subtraction, we subtracted the continuum obtained from line-free channels in all four spectral units. All narrow-band spectra were resampled to a spectral resolution of 0.5 km s −1 and when imaged with the CLARK algorithm and uniform weighting have a negligibly smaller synthesised beam than the continuum images. The average rms noise of the line images in the ABD configuration is 11.2 mJy beam −1 km s −1 . The synthesised beam size and the average rms noise of the line data for all imaged combinations of array configurations are listed in Table 4. Notes. Rest frequencies and upper energy levels have been obtained from the Cologne Database for Molecular Spectroscopy (CDMS) (Müller et al. 2001; with the exception of those for HCOOCH 3 transitions, which were acquired from the Jet Propulsion Laboratory (Pickett et al. 1998).

IRAM 30m observations
Observations of W3(H 2 O) with the 30m telescope were obtained on March 13, 2015, centred on the same position as the phase centre of the interferometric observations. We used the Eight MIxer Receiver (EMIR) covering the range 213 − 236 GHz, reaching a spectral resolution of 0.3 km s −1 . In this work, we have merged the NOEMA observations of 13 CO with the singledish observations using the mapping software and CLEANed the merged cube with the Steer-Dewdney-Ito (SDI) method (Steer et al. 1984) in order to recover more of the extended features to study molecular outflows. Further details of the 30m observations and data reduction, as well as the merging process, can be found in Mottram et al. (in prep.). The resulting merged image has an angular resolution of 1 . 14 × 0 . 92, PA=49 • , and an rms noise of 8.4 mJy beam −1 km s −1 . We also make use of our single-dish 13 CO data which have been reduced and converted to brightness temperatures for a detailed outflow analysis (see Section 3.3).

Observational results
In the following, we present our detailed analysis for W3(H 2 O), and when applicable, we also showcase our observational results for W3(OH). Our analysis mainly uses the continuum and CH 3 CN spectral line emission. Maps of the other lines are shown in Appendix A.    Figure 3 shows the 1.37 mm (219 GHz) continuum emission map of W3(H 2 O) and W3(OH) in the ABD configuration. At this wavelength, the continuum emission in our field of view is dominated in W3(OH) by free-free emission, while the emission in W3(H 2 O) is due to dust (Wyrowski et al. 1999). In the following, we focus on the fragmentation and kinematics of the younger region, W3(H 2 O). Figure 4 shows a comparison of the continuum emission maps of W3(H 2 O) obtained by imaging the ABD, AB, and Aarray only observations. The integrated flux within the 6σ contours is 1220 mJy, 656 mJy, and 364 mJy for the ABD, AB, and A-array images, respectively. The fragmentation of W3(H 2 O) into two cores, separated by ∼2300 AU, is best seen in the AB image at 700 AU scales. The two cores are labelled W3(H 2 O) E and W3(H 2 O) W, and their peak continuum positions are depicted by stars in Fig. 4, marking the positions of embedded (proto)stars. The peak position of W3(H 2 O) E is α(J2000) = 02 h 27 m 04 s .73, δ(J2000) = 61 • 52 24 . 66, and that of W3(H 2 O) W is α(J2000) = 02 h 27 m 04 s .57, δ(J2000) = 61 • 52 24 . 59. The approximate separation boundary between W3(H 2 O) E and W is marked with a vertical dashed line. The integrated flux within 6σ contours and the separation boundary are 377 mJy and 279 mJy for W3(H 2 O) E and W3(H 2 O) W, respectively. Furthermore, there is an additional emission peak to the north-west of W3(H 2 O) at an offset of −1 . 6, 1 . 6 (α(J2000) = 02 h 27 m 04 s .37, δ(J2000) = 61 • 52 26 . 35), which is best seen in the ABD image as it has the best sensitivity. This is most likely a site for the formation of lower mass stars.

Continuum emission
Radiative transfer models by Chen et al. (2006) for W3(H 2 O) in the 1.4 mm wavelength regime show that the averaged dust optical depth is less than 0.09; therefore, we assume the thermal dust emission to be optically thin in our observations. The positions of continuum peaks A and C from Wyrowski et al. (1999) coincide well with the continuum peaks W3(H 2 O) E and W3(H 2 O) W in our observations to within a synthesised beam (see Fig. 4). Wyrowski et al. had attributed the third peak in their observations between the other two cores to an interplay of high column density and low temperatures in the central region, and concluded that W3(H 2 O) is harbouring two cores at the positions of continuum peaks A and C. Our double-peaked continuum image in AB, with a better spatial resolution than that of Wyrowski et al. (1999) by a factor of 3.3, supports this interpretation. Furthermore, peak B coincides well with the north-east extension of W3(H 2 O) W in our observations, best seen in the middle panel of Fig. 4.

Line emission
W3(H 2 O) is one of the most chemically rich sources in our sample (see Fig. 1) with detections of sulphur-bearing species such as 33 SO and 34 SO 2 , complex species such as HCOOCH 3 , and vibrationally excited lines of HC 3 N, among many others. Figure 5 shows integrated intensity (zeroth moment) and intensityweighted peak velocity (first moment) maps of CH 3 CN (12 3 − 11 3 ) for W3(H 2 O) and W3(OH). The zeroth moment map confirms the fragmentation of W3(H 2 O) into two cores. The moment maps of most lines covered by the narrow-band receiver (see Table 3 and Fig. 2) are presented in Appendix A. All moment maps have been created inside regions where the signal-tonoise ratio is greater than 5σ. The integrated intensity maps of most tracers for W3(H 2 O) also show two peaks coincident with the locations of W3(H 2 O) E and W3(H 2 O) W. While the continuum emission is stronger for W3(H 2 O) E, some dense gas tracers (e.g. CH 3 CN, HC 3 N) show stronger line emission towards The bottom panels in Figure 5 show the intensity-weighted peak velocity (first moment) map of the region in CH 3 CN (12 3 − 11 3 ). We chose to do our kinematic analyses on this transition as it is the strongest unblended line in the methyl cyanide (CH 3 CN) K-ladder. There is a clear velocity gradient in the east-west direction across W3(H 2 O), and in the NW-SE direction across W3(OH). The systemic velocities of both clumps are determined by averaging the spectra of CH 3 CN (12 3 − 11 3 ) over a 4 × 4 area centred on each source and fitting a Gaussian line to the resulting averaged spectrum. In this way, W3(H 2 O) and W3(OH) have average velocities of -49.1 and -45.0 km s −1 , respectively.
The velocity gradient across W3(H 2 O) is detected in most of the high spectral resolution lines in our survey (see Fig. A.2 in Appendix A) and spans ∼6000 AU in size, corresponding to an amplitude of 170 km s −1 pc −1 . The velocity gradient resolved in W3(OH) has an amplitude of ∼100 km s −1 pc −1 and is roughly perpendicular to the motion of the ionised gas in the east-west direction as traced by the H92α line in observations of Keto et al. (1995), and is perpendicular to the direction of the 'champagne flow' observed to the north-east at radio frequencies (Keto et al. 1995;Wilner et al. 1999). Thus, in W3(OH) we seem to be witnessing the large-scale motion of the remnant molecular gas. The interpretation of line emission for W3(OH) is complex because most of the continuum emission at 1.37 mm is due to freefree emission which affects the appearance of molecular lines. In practice, free-free emission at 1.3 mm reduces the molecular line emission which is reflected by the reduced integrated line emission towards the peaks of W3(OH). As this is beyond the scope of this paper, we refrain from further analysis of W3(OH).

Outflow structure
In Figure 6 we show integrated intensity (zeroth moment) maps of outflow-tracing molecules ( 12 CO and 13 CO) for the redshifted and blueshifted gas. The minimum intensity below which a pixel is not considered in the creation of the moment maps is based on 5σ rms noise level in emission-free channels. The single-dish CO (2-1) map (see top left panel of Fig. 6) shows the existence of a bipolar outflow in the overall vicinity of W3(H 2 O) but also encompassing W3(OH) in an approximately NE-SW direction. The 13 CO (2-1) single-dish map (see top right panel of Fig. 6) shows less extended emission than the 12 CO as the 13 C isotopologue is a factor of ∼76 less abundant (Henkel et al. 1982), and the map highlights the general NE-SW direction of the outflow even better.
The integrated intensity map of 12 CO from SMA interferometric data of Zapata et al. (2011) allows for the detection of two bipolar outflows (see bottom left panel of Fig. 6). The outflow emanating from W3(H 2 O) E has its blueshifted side to the south-west and its redshifted lobe to the north-east, while the second outflow emanating from W3(H 2 O) W has its blueshifted side extending to the north-east with its redshifted side to the south-west . The difference between the position angles of the two outflows (in the plane of the sky) is 25 • . Furthermore, the resulting zeroth moment map of 13 CO emission from our combined NOEMA and 30m single-dish observations, presented in the bottom right panel of Fig. 6, confirms the findings of Zapata et al. (2011) with regards to the directions and origin of the redshifted outflow lobe from W3(H 2 O) W and the origins of the blueshifted outflow lobe from W3(H 2 O) E. However, we miss much of the emission that is detected in the 12 CO SMA interferometric data, mainly due to the lower abundance of the 13 C isotopologue, and thus its lower sensitivity to the outflowing gas. The same coloured arrows obtained from Zapata et al. are redrawn in a zoom panel inside the bottom right panel of this figure, highlighting that the two outflows are in fact emanating from different positions.
In Figure 7, we show how the centimetre emission aligns with the millimetre continuum emission. The directions of bipolar molecular outflows emanating from the cores are shown by red and blue arrows, and the positions of water masers shown by yellow triangles. The elongated radio source centred on W3(H 2 O) E has a spectral index of -0.6 and is a source of synchrotron emission, best described by a jet-like model due to its morphology Wilner et al. 1999). The radio source centred on W3(H 2 O) W, however, has a rising spectral index, possibly due to a circumstellar wind being ionised by an embedded (proto)stellar source at this position. Furthermore, the discrepancy between the direction of the synchrotron jet-like object and the molecular outflow could be due to the existence of multiple objects in this core, unresolved by our observations. Interestingly, it has been shown that synchrotron radiation can be produced not only via jets, but also through the acceleration of relativistic electrons in the interaction of disk material with a stellar wind (Shchekinov & Sobolev 2004), providing an alternative explanation for the maser and cm emission. Top left: Intensity map of CO (2-1) emission from IRAM 30m integrated over the velocity range from −40 to −20 km s −1 for the redshifted and from −80 to −55 km s −1 for the blueshifted gas. Contours start at 5σ and increase in steps of 5σ (1σ = 1.9 K km s −1 ). Top right: Intensity map of 13 CO (2-1) emission from IRAM 30m integrated over the velocity range from −40 to −20 km s −1 for the redshifted and from −80 to −55 km s −1 for the blueshifted gas. Contours start at 3σ and increase in steps of 2σ (1σ = 0.42 K km s −1 ). Bottom left: Intensity map of CO (2-1) emission obtained by Zapata et al. (2011) with the SMA integrated over the velocity range from −45 to −25 for the redshifted and from −75 to −55 for the blueshifted gas. Contours start at 5σ and increase in steps of 5σ (1σ = 1.1 mJy beam −1 km s −1 ). The blue and red arrows highlight the positions and directions of the CO outflows emanating from each source. Bottom right: Intensity map of 13 CO (2-1) emission from IRAM 30m merged with NOEMA observations integrated over the velocity range from −41 to −30 for the redshifted and from −70 to −53 for the blueshifted gas. Contours start at 3σ and increase in steps of 3σ (1σ = 0.03 mJy beam −1 km s −1 for redshifted and 0.17 mJy beam −1 km s −1 for blueshifted sides). The zoom panel in the top left corner highlights the launching positions of each outflow, with the red and blue arrows showing the directions of the two bipolar outflows from Zapata et al. (2011). The rms noise used in drawing the contours of the integrated intensity maps have been determined by first creating the maps without any constraints on the minimum emission level (threshold of 0) and calculating the noise in an emission-free part of the resulting map.

Dense gas kinematics
The kinematics of W3(H 2 O) can be further analysed by looking at position-velocity (PV) diagrams for various transitions of dense gas and potentially disk-tracing molecules such as CH 3 CN and HCOOCH 3 (e.g. Beuther et al. 2005). In the following, we divide our focus between the large-scale kinematics of the entire W3(H 2 O) region, where we put forward arguments for the observed velocity gradients in CH 3 CN and HCOOCH 3 being due to rotation instead of infall (Section 4. tion 4.1.2). Although the alignment of the elongated radio emission with the water masers can be described by an outflow model oriented in the east-west direction, the detection of two molecular outflows emanating from the positions of the two continuum peaks, roughly perpendicular to the observed velocity gradient in CH 3 CN at large scale and perpendicular to the observed velocity gradients on smaller scales (see Section 4.1.2), make it unlikely that the motions in CH 3 CN would be due to expansion or outflow. Furthermore, CH 3 CN and HCOOCH 3 line profiles do not show the broad components typically seen in emission from expanding gas, and these species are nevertheless too complex to exist in an ionised jet. We therefore conclude that the observed velocity gradient is most likely due to rotation, which we will assume for the remainder of this paper.

Large-scale kinematics
The PV plots of CH 3 CN and HCOOCH 3 for W3(H 2 O) are shown in Fig. 8 for a cut in the direction of the velocity gradient going through the continuum peaks (dashed line in bottom left panel of Fig. 5) obtained from the NOEMA observations. The white curves in the top left panel correspond to gas in Keplerian rotation with V rot ∝ R −1/2 about a 10, 25, and 50 M central object. These white curves are not fits to the PV diagram, but are merely drawn to guide the eye. It is clear that the gas is not in Keplerian rotation; however, higher velocity gas is observed closer to the centre of W3(H 2 O) which can be a signpost for differential rotation. Ohashi et al. (1997) created models for comparison to their interferometric data of the low-mass protostar L1527, using a thin disk with 2000 AU extent configured edge-on and present PV diagrams for cases with various degrees of infall and rotation. In the case of infall-only motions, their PV plots are axisymmetric with two peaks offset symmetrically in the velocity axis. With the addition of rotation, the peaks become blueshifted and red- shifted away from the central positions such that in the case of a pure Keplerian rotation and in the absence of infall we would recover the classical butterfly-shaped rotation curve. Comparing our PV plots to the Ohashi et al. (1997) scenarios, we observe minimal contributions in the top left and bottom right quadrants of our plots where one would expect some emission in the cases including infall motions. Therefore, we do not detect infall motions in our interferometric data probably because the infalling envelope is too diffuse and filtered out. Furthermore, models by Tobin et al. (2012) for spherical rotating collapse and filamentary rotating collapse showed similar results to that of Ohashi et al. (1997), confirming that the absence of infall results in the lack of emission in those quadrants and not the projection or source morphology. Moreover, the PV plots of HCOOCH 3 , which is a less abundant species than CH 3 CN, are more representative of rigid body-like rotation. Chen et al. (2006) generated a binary model with radiative transfer post-processing of methyl cyanide for each source within W3(H 2 O), showing that the high-velocity deviations from solid body rotation in their PV plots could be a result of two spatially unresolved cores (on similar scales to our observations) with a small radial velocity difference. The detection of a smooth velocity gradient in all tracers across W3(H 2 O) E and W suggests that the observed angular velocity difference cannot be solely due to binary motion of the cores, but that there is also a large-scale toroidal structure encompassing and rotating about the two resolved cores. Such circumbinary toroidal structures have indeed been observed previously at lower angular resolutions in other sources (e.g. Beltrán et al. 2005;Beuther et al. 2007).

Small-scale kinematics
Although we see a smooth velocity gradient in CH 3 CN across the entire W3(H 2 O) structure, the existence of at least two cores within it as presented by our continuum results, and two collimated outflows requires further analysis. Imaging the data in the most extended configuration and therefore reaching our highest resolution, we can filter out the large-scale envelope to analyse the kinematics of gas around each of the cores. Figure 9 shows the first moment map of CH 3 CN (12 3 − 11 3 ) for the cores to the east and west using the A-array observations exclusively. The maps have been scaled and masked to highlight the velocity structure of each core. The mean line-of-sight velocity difference between the two main cores, W3(H 2 O) W and E, is a few km s −1 . Velocity differences of a few km s −1 are observed across each core, approximately perpendicular to the directions of the bipolar molecular outflows emanating from each core . This indicates that the small-scale gradients across each individual core are most likely due to rotation. The two line-ofsight velocity gradients are of the order of ∼1000 km s −1 pc −1 , depending on the choice for the extent of the gradient, much faster than the rotational motion of 170 km s −1 pc −1 for the whole of W3(H 2 O) from the ABD data. Furthermore, the directions of the velocity gradients observed for each core are slightly inclined with respect to the overall east-west motion of the gas on larger scales (Fig. 5), suggesting that the general rotation seen around the two cores may be inherited from the large-scale rotation. Meanwhile, the direction of the blueshifted (redshifted) outflow emanating from W3(H 2 O) W is almost in the opposite direction of the blueshifted (redshifted) outflow ejected from W3(H 2 O) E. This implies that the inclination angles of the two rotating structures with respect to the plane of the sky are likely different. Figure 10 shows PV diagrams for W3(H 2 O) E (left) and W3(H 2 O) W (right) corresponding to cuts in the directions of rotation as depicted by dashed lines in Figure 9. Based on the PV plots, the LSR velocities of W3(H 2 O) E and W3(H 2 O) W are estimated to be -51 km s −1 and -47 km s −1 , respectively. White curves correspond to gas in Keplerian rotation about a 5, 10, and 15 M centrally dominated object. As on larger scales, these PV diagrams do not show the symmetric four-quadrant shape expected if the velocity gradients were due to infall.
The PV plot for W3(H 2 O) W contains contributions from W3(H 2 O) E due to the angle and extent of the cut, hence there is added emission in quadrants toward positive offset, making it difficult to infer rotational signatures pertaining to the blueshifted emission of W3(H 2 O) W. The redshifted rotational signatures seen in the quadrants toward negative offsets, however, show signatures of increased gas velocities closer to the centre of the core. Such a trend in the PV plot implies differential rotation of material, possibly within a disk-like object.
Position-velocity plots for W3(H 2 O) E have a lower signalto-noise ratio than W3(H 2 O) W as line emission is typically weaker for this fragment even though its 1.37 mm dust continuum peak is stronger. The linearity of the rotation curves in Fig. 9. Intensity-weighted peak velocity (first moment) map of CH 3 CN (12 3 − 11 3 ) using only the A-array observations and masked out to show contributions from W3(H 2 O) W (top) and W3(H 2 O) E (bottom). The solid contours correspond to the 1.37 mm continuum in the A-array only observations and start at 6σ and increase in steps of 3σ (1σ = 2.5 mJy beam −1 ). The dashed lines correspond to the cuts made for the PV plots (Fig. 10). The blue and red arrows show the directions of bipolar molecular outflows (Fig. 6). A scale bar and the synthesised beam (0 . 39 × 0 . 28, PA=88 • ) are shown in the bottom. The velocity ranges are different for the two cores. CH 3 CN does not reveal Keplerian signatures, but is more consistent with rigid body-like rotation. In order to increase the signal-to-noise ratio, we stacked the PV plots of CH 3 CN (12 K − 11 K ) K = 2 − 5 and show this stacked PV plot in the bottom right panel of Fig. 10 for W3(H 2 O) E. Stacking is a reasonable strategy as the variation in average linewidths of these transitions around this core is below our spectral resolution; therefore, assuming these lines to be probing the same surface is valid. In the stacked image, the 4σ contour reveals a high-velocity feature close to the centre of the core toward positive offsets. As this feature has an extent of the order of our spatial resolution, it is unclear whether the rotation observed for W3(H 2 O) E is due to a disk-like object, an unresolved binary (or multiple) system, or a combination of the two.

Temperature distribution
As a symmetric top molecule, CH 3 CN is an excellent thermometer of hot molecular gas (e.g. Loren & Mundy 1984;Zhang et al. 1998) since its relative populations in different K-levels are dominated by collisions. Our high spectral resolution observations of CH 3  their isotopologues which have upper energy levels in the range 70-325 K with 0.5 km s −1 spectral resolution (see Table 3).
We made use of the eXtended casa Line Analysis Software Suite (xclass 5 , Möller et al. 2017) to model the observed spectra under the assumption of Local Thermodynamical Equilibrium (LTE), which is typically valid for CH 3 CN in such high-density environments (see Section 4.3). In summary, xclass solves the radiative transfer equation in one dimension for a set of initial conditions (source size, column density, temperature, linewidth, and peak velocity) provided by the user, and through a χ 2 minimisation routine changes the initial conditions within ranges that have been provided by the user to obtain the best fit to the observed spectra. The details of our xclass modelling are summarised in Appendix B.
In Figure 11, we present our results of pixel-by-pixel xclass modelling of CH 3 CN (12 K − 11 K ) K = 4 − 6, including CH 3 13 CN (12 K − 11 K ) K = 0 − 3, in AB configuration which produces rotational temperature, column density, peak velocity, and linewidth maps for CH 3 CN. The column density map is double peaked, similar to the continuum emission, although the column density peaks are slightly offset to the north-west by a synthesised beam. This offset is consistent with the offset found between the continuum peaks and the integrated intensity maps of CH 3 CN (12 3 − 11 3 ) (see Fig. 5) and most high-density trac-5 https://xclass.astro.uni-koeln.de ers (see Fig. A.1 in Appendix A). The median CH 3 CN column density is 1.4 × 10 15 cm −2 . The velocity gradient derived in this way is consistent with the large-scale east-west velocity gradient observed in the first moment map of CH 3 CN, confirming that our modelling strategy captures this accurately, regardless of its origin. The linewidths are also larger by a few km s −1 for W3(H 2 O) W, also seen in the second moment maps. The median rotational temperature of W3(H 2 O) is ∼165 K and the temperature structure does not follow any particular pattern within either core. Two high-temperature features can be seen which may be associated with regions carved out by the molecular outflows allowing a deeper look into the cores, or regions that may have been additionally heated by the outflow. Nevertheless, the compactness of these features, which are of the same order as the size of our synthesised beam, along with decreased signal-tonoise ratio at the edges of the map, prevent us from making firm conclusions in this regard.

Mass estimates
The combined bolometric luminosity of W3(H 2 O) and W3(OH), determined from fitting the SED from the near-IR to sub-mm, is 8.   Table 4). Bottom right: observed spectrum of a given pixel drawn in black and overlaid with the resulting fit for CH 3 CN (12 K − 11 K ) K = 4 − 6 in red and the marginally detected where F ν is the flux density of the free-free radio continuum emission at frequency ν, T e is the electron temperature, and d is the distance. The observed integrated radio flux at 15 GHz for W3(OH) is 2.53 Jy (Kurtz et al. 1994), and assuming a typical electron temperature of 10 4 K, this results in a value of Q 0 = 1.2×10 48 photons s −1 . Interpolating from Table 1 of Davies et al. (2011) and the relationships between spectral type and photospheric temperature from Martins et al. (2005) and Boehm-Vitense (1981) for O and B stars, respectively, this ionising photon flux approximately corresponds to an O8.5 main sequence star, with M ≈ 20 M and L ≈ 4.4 × 10 4 L . This leaves a total bolometric luminosity of 3.9 × 10 4 L , which we assume to be evenly distributed between the two cores within W3(H 2 O), 1.95 × 10 4 L each; using the same tables this would correspond to two 15 M stars of spectral type B0. The above estimates are based on Davies et al. (2011) for zero age main sequence (ZAMS) stars. High-mass protostars growing by accretion resemble ZAMS stars in terms of their luminosities and temperatures when core nuclear burning dominates other sources of luminosity such as accretion and envelope burning. When this occurs depends primarily on when the protostar gains sufficient mass, but also on the accretion rate. Stellar structure calculations suggest that this occurs at 5 to 10 M for accretion rates of 10 −5 to 10 −4 M yr −1 , respectively (Norberg & Maeder 2000;Behrend & Maeder 2001;Keto & Wood 2006). The mass estimates further assume that all the emitted energy has been produced within the (proto)stars, ignoring contributions from episodic accretion to the luminosity. Therefore, the 15 M estimates can be taken as upper limits, in agreement with calculations of Chen et al. (2006) who find a minimum binary mass of 22 M for the protostars within W3(H 2 O). Furthermore, assuming the gas to be in gravito-centrifugal rotation around the two individual cores, our PV plots (see Fig. 10) suggest that the range of 5 − 15 M is a reasonable estimate for the protostellar masses.
As dust is typically optically thin in the 1.3 mm wavelength regime and proven to be for this region in particular (Chen et al. 2006), we use the prescription by Hildebrand (1983) to convert the flux density, F ν , of the continuum observations to a mass. In the form presented by Schuller et al. (2009), where R is the gas-to-dust mass ratio, B ν (T D ) is the Planck function at a dust temperature of T D , and κ ν is the dust absorption coefficient. We adopt a gas-to-dust mass ratio of R = 150 (Draine 2011) and a value of κ ν = 0.9 cm 2 g −1 for the dust absorption coefficient from Ossenkopf & Henning (1994), corresponding to thin ice mantles after 10 5 years of coagulation at a density of 10 6 cm −3 .
High-mass cores such as the ones we are studying typically have densities high enough to thermalise the methyl cyanide lines. To check this, using the spontaneous decay rate of CH 3 CN (12 4 − 11 4 ) obtained from the LAMDA database 6 , 7.65 × 10 −4 s −1 , and the collision rate of 2.05 × 10 −10 cm 3 s −1 with H 2 at 140 K (Green 1986), we calculate the simplified twolevel critical density of this line to be n crit ≈ 3.7 × 10 6 cm −3 . The effective density, once the radiation field is taken into account, is typically an order of magnitude lower (Evans 1999;Shirley 2015). Following Schuller et al. (2009), the H 2 column density is calculated via where S ν is the peak intensity, θ B is the beam solid angle, µ is the mean molecular weight assumed to be equal to 2.8, and m H is the mass of the hydrogen atom. At a temperature of 140 K and using our continuum intensity of a given position near the centre, we can estimate the H 2 column density to be 4.5×10 24 cm −2 . This can be converted to a volume density of n H 2 > 7×10 7 cm −3 , assuming the maximum extent of the third dimension to be the plane-of-sky size of the clump (∼4000 AU). Therefore, since the density of molecular hydrogen is much higher than the critical density of the lines, the LTE assumption is valid and the rotational temperature map of CH 3 CN obtained from our modelling can be assumed to be tracing the gas kinetic temperature. Using Eq. 2, our continuum map with its unit converted to Jy pixel −1 and the temperature map obtained from xclass, assuming dust and gas temperatures are coupled, we obtain a pixel-by-pixel mass density map (see Fig. C.1 in Appendix C). Summing over the pixels in our mass density map in the ABD observations, the total mass for W3(H 2 O) is calculated to be ∼26.8 M , with 15.4 M contributed from the core to the east, and 11.4 M from the core to the west. Similarly, using the AB observations, we obtain a total mass of ∼11.4 M for W3(H 2 O), with a core mass of 6.7 M and 4.7 M from the cores to the east and west, respectively. The effect of spatial filtering of the interferometer is clear in these mass estimates as the exclusion of the D-array data removes more than half of the mass.
Comparing our NOEMA observations to SCUBA-2 850 µm single-dish observations, about 25% of the flux is filtered out by the interferometer in our ABD observations (assuming a ν −3.5 frequency-relation, Beuther et al. 2018), implying that our core mass estimates of 15.4 and 11.4 M are lower limits, but nevertheless reasonable. Masses estimated in this manner have contributions from the cores and the disk-like structures and not from any embedded (proto)stars.

Toomre stability
For a differentially rotating disk, the shear force and gas pressure can provide added stability against gravitational collapse. This idea was originally introduced by Safronov (1960) and further quantified by Toomre (1964) for a disk of stars, and has since been used in various applications ranging from planet formation to galaxy dynamics. We investigate the stability of the rotating structures in W3(H 2 O), assuming that they are disks in gravitocentrifugal equilibrium, against axisymmetric instabilities using the Toomre Q parameter, where c s is the sound speed, and Ω is the epicyclic frequency of the disk which is equivalent to its angular velocity. The surface 6 http://home.strw.leidenuniv.nl/~moldata/ density of the disk, Σ, is calculated by multiplying the column density (Eq. 3) by the mean molecular weight and mass of the hydrogen atom (µm H ) to convert the number column density to a mass surface density. A thin disk becomes unstable against axisymmetric gravitational instabilities if Q < 1.
Having obtained a temperature map representative of the kinetic temperature, we are able to calculate the Toomre Q parameter pixel by pixel. In particular, the temperature is used in the calculation of B ν (T D ) and the sound speed, where γ is the adiabatic index with a value of 5/3. We assume two 10 M (proto)stellar objects to be present, one at the position of each of the two continuum peaks (see AB image in Fig. 4). The angular velocity maps are created by adding up mass within concentric circles starting at the position of each core and going outwards. In this way, we incorporate the mass of the rotating structures, calculated via Eq. 2, and the mass of the central objects. The angular velocity of a disk in gravitocentrifugal equilibrium at a radius r is where the mass of the central object, M * , is 10 M and M disk (r) is the gas mass enclosed within r. Given the radii involved, this means that in practice, most parts of the maps are dominated by M disk rather than M * . In Figure 12 we present our Toomre Q map of W3(H 2 O), which is created by stitching together the Toomre Q maps of the two individual cores. The stitching boundary is shown by a solid vertical line and the positions of the two central objects corresponding to our continuum peaks in the AB image (see Fig. 4) are indicated by stars. While the Toomre Q calculations were performed on the AB-array data, we have drawn the continuum contours from the A-array only observations. The boundary where Q = 1 is shown by a dashed line.
The most significant factor stabilising the disk against Toomre instability is the high gas temperatures and the fast differential rotation of material closest to the (proto)star; therefore, we find the highest Toomre Q values closest to the presumed locations of the (proto)stars depicted by stars in Fig. 12 as expected. We find low Toomre Q values at the periphery of W3(H 2 O) E and W, suggesting that both rotating structures are only marginally stable against axisymmetric instabilities. Moreover, W3(H 2 O) W shows low Toomre Q values at the positions of additional continuum peaks in the highest resolution A-array data (solid contours in Fig. 12), suggesting further fragmentation at these positions may be possible. As the dimensions of the observed continuum peaks are smaller than the size of our synthesised beam, the validity of this hypothesis needs to be confirmed with higher resolution observations (at 345 GHz) to be taken with NOEMA.
A Toomre stability analysis performed by Chen et al. (2016) for the high-mass protostar IRAS 20126+4104 showed that the massive disk around that protostar is hot and stable to fragmentation, with Toomre Q > 2.8 at all radii. Hydrodynamic simulations of collapsing protostellar cores by Klassen et al. (2016) find massive accretion disks which evolve to become asymmetrically Toomre-unstable, creating large spiral arms with an increased rate of accretion of material onto the central object. They find that the self-gravity of these spiral arm segments is lower than the tidal forces from the star, causing the spiral arms in their simulations to remain stable against collapse. However, they find that the Toomre conditions combined with the cooling of the disk (Gammie 2001) would potentially yield the formation of a binary companion. Similarly, simulations by Kuiper et al. (2011) see such spiral arms in a massive disk, and simulations by Krumholz et al. (2009) and more recently Meyer et al. (2017Meyer et al. ( , 2018 see disk fragmentation on even smaller spatial scales, of the order of hundreds of AU. Therefore, disk fragmentation scenarios in the high-mass regime are valid from a theoretical stand-point. Among the assumptions that go into this Toomre analysis is the thin-disk approximation. According to Gammie (2001), the critical value of Q for an isothermal disk of finite thickness should be 0.676 (as derived from Goldreich & Lynden-Bell 1965). Even with such a low critical value, the rotating structures around each core would still show signs of instability in their outskirts. To get an estimate on the importance of disk thickness, following Gammie (2001), the instability condition for a disk with scale height, H c s /Ω, can be written as Using our temperature and angular velocity maps, we create a map of the H/r ratio. In this way, we calculate an average value of H/r ≈ 0.3. Assuming the mass of each (proto)star is roughly 10 M , we calculate that for the candidate disks to be unstable against gravitational collapse, the mass of the candidate disk must be M disk 3 M . Since the mass estimates for the cores obtained from the continuum maps using Eq. 2 are 15 and 11 M , the masses are large enough to result in the rotating structures being Toomre-unstable. Therefore, the thickness of the candidate disks does not affect the general interpretation of our results. The angular velocity maps used in the calculation of the Toomre Q map have been created assuming a face-on geometry for the candidate disks. If the objects are inclined, we would be overestimating the surface density of the candidate disks and the line-of-sight angular velocities. Since these two parameters have the opposite effect on the Toomre Q value (see Eq. 4), our results are still reasonable. Furthermore, W3(H 2 O) E may harbour multiple (proto)stars within it (see Section 4.1) such that the possibility of rotational contributions from multi-body rotation would require a different prescription for the epicyclic frequency than what we have assumed.
The estimation of 10 M for the central (proto)stars have been deduced from our own and previous mass estimations (see Section 4.3), and by examining the emission in the PV plots with Keplerian curves that are just 'grazing' the emission. The latter method has the caveat that the observed PV plot is the result of the convolution with the instrumental beam and the line width. Therefore, the real emission is less extended in space and velocity than we observe, and as a result we may be overestimating the mass of the (proto)stars. Performing our Toomre calculations, but with the assumption of 5 M (proto)stars, we obtain median Toomre Q values of 1.3 compared to a value of 1.7 in the 10 M calculations. Therefore, lower protostellar masses would lead to the decrease in the Toomre Q value, and thus larger unstable regions, making the disk fragmentation scenario even more plausible. On the other hand, as the mass estimates using the bolometric luminosities suggest the central (proto)stellar masses to be in the 15 M regime, the median Toomre Q value would be 2.0, with unstable regions still present in the outskirts. The Toomre Q maps created assuming two 5 and 15 M (proto) stars at the position of the continuum peaks can be found in Appendix D. The change of ±5 M does not strongly affect the Toomre Q maps because we add the mass of the candidate disks, estimated from the continuum emission to be 15 and 11 M , in rings (see Eq. 6) and they dominate over an addition or subtraction of 5 M , especially in the outermost regions. This, together with the root dependence on mass lead to uncertainties in M * to only have a small effect on the Toomre Q results.
We are currently conducting simulated observations from 3D hydrodynamic simulations to investigate the reliability of our methods in even more depth (Ahmadi et al. in prep.). Our preliminary results on the effect of inclination and (proto)stellar mass on the calculation of Toomre Q maps using hydrodynamic simulations suggest that the Toomre Q value is not strongly dependent on these parameters. A forthcoming paper is dedicated to this topic.

Cooling
For a marginally stable disk (Q ∼ 1), the locally unstable regions compress as they start to collapse, providing compressional heating in these regions. This increase in the local thermal pressure can counteract the local collapse and put the disk back into a state of equilibrium such that the region will not fragment unless heat can be dissipated on a short timescale. Through numerical simulations, Gammie (2001) and Johnson & Gammie (2003) have introduced a cooling parameter, β, to study the effect of cooling for marginally stable thin disks in gravito-centrifugal equilibrium. In their prescription where the cooling time, t cooling , is assumed to be constant and a function of the surface density and internal energy of the disk. If the local collapse needs more than a few orbits, due to a long cooling timescale, a locally collapsing region within the disk can be ripped apart by shear. Therefore, there is a critical β value below which a locally collapsing region would be rapidly cooling, and if it is sufficiently self-gravitating, it would be prone to fragmentation. Conversely, values above this critical β value would put a marginally stable disk into a self-regulating scenario such that heating acts as a stabilising force and directly counteracts the cooling rate. The Gammie cooling criterion is relevant only if the disk is marginally stable, with regions that are locally unstable against axisymmetric gravitational instabilities. A Toomre-stable disk will not fragment regardless of how quickly it may be cooling.
This critical value is determined through numerical simulations and varies depending on the heating and cooling recipes and convergence issues. For an overview of these issues, see Kratter & Lodato (2016) and Baehr et al. (2017). Regardless of the specific numerical simulation issues, it is generally assumed that the critical β value ranges between 1 − 5.
Following the prescription of Klassen et al. (2016), where E z is the column internal energy integrated in the zdirection (along the axis of rotation) and F grey is the radiative flux away from the disk adopting a grey-body approximation. The internal energy of an ideal gas is calculated via where ρ is the volume density and c v is the specific heat capacity defined as with N A as Avogadro's number and M u the molar mass. By replacing the volume density with the surface density, the column internal energy is then Knowing the cooling time, which is of the order of months, and the rotation rate at each pixel, we construct a local β map. Figure 13 shows the β map for the two disks presumed in gravitocentrifugal equilibrium about the positions of peak continuum emission as depicted by stars. The median β value is 1.8 × 10 −4 , and the maximum value found in the outskirts is of the order of 10 −2 , which is much lower than the critical value of 1-5. This finding is consistent with that of Klassen et al. (2016) and Meyer et al. (2018) who also find β 1 in their simulation of collapsing protostellar cores with a mass of 100 M at later times. The rotating bodies within W3(H 2 O) are therefore able to cool rapidly and any local collapse induced by the gravitational instabilities in the disks will lead to fragmentation. Rapid cooling in this context means that the cooling timescale is comparable to the local dynamical timescale.

Conclusions and summary
Our IRAM large programme, CORE, aims to study fragmentation and the kinematics of a sample of 20 high-mass star-forming regions. We have performed a case-study for one of the sources in our sample, the prototypical hot core W3(H 2 O). In this paper we present details of the spectral line set-up for our NOEMA observations in the 1.37 mm band which covers transitions of important tracers of dense gas and disks (e.g. CH 3 CN, HCOOCH 3 , CH 3 OH), inner envelopes (e.g. H 2 CO), and outflows (e.g. 13 CO, Fig. 13. Map of β cooling parameter obtained by assuming two disklike structures in gravito-centrifugal rotation about the positions of peak continuum emission as depicted by stars. The solid contours correspond to our continuum data in the most extended configuration, starting at 6σ and increasing in steps of 3σ (1σ = 2.5 mJy beam −1 ). The solid vertical line corresponds to the stitching boundary. The dashed line corresponds to Toomre Q = 1. Regions outside of the 6σ mm continuum emission contour in the AB configuration are masked out. SO). We cover the range of 217 078.6 − 220 859.5 MHz with a spectral resolution of 2.7 km s −1 , and eight narrower bands centred on specific transitions to provide higher spectral resolution of 0.5 km s −1 .
With the aim of studying the fragmentation and kinematic properties of W3(H 2 O), the following is a summary of our findings: -At an angular resolution of ∼0 . 35 (∼700 AU at 2 kpc), W3(H 2 O) fragments into two cores, which we refer to as W3(H 2 O) E and W3(H 2 O) W, separated by ∼2300 AU, as seen in both line and thermal dust continuum emission; -Based on the integrated dust continuum emission, W3(H 2 O) has a total mass of ∼26.8 M , with 15.4 M contributed from W3(H 2 O) E, and 11.4 M from W3(H 2 O) W; -On large scales, there is a clear velocity gradient in the eastwest direction across W3(H 2 O), spanning ∼6000 AU, attributed to a combination of circumbinary and circumstellar motions. On smaller scales, velocity gradients of a few km s −1 are observed across each of the cores, perpendicular to the directions of two bipolar molecular outflows, one emanating from each core, as traced by 12 CO and 13 CO. The direction of motion of the gas around the individual cores deviates little from the overall larger scale rotation of W3(H 2 O), suggesting that these motions seen around the cores, which we interpret as rotation, may be inherited from the largescale rotation; -The kinematics of the rotating structure about W3(H 2 O) W shows differential rotation of material about a (proto)star as deduced from the redshifted part of its PV plot, suggesting that the rotating structure may be a disk-like object. The radio source with a rising spectrum at this position can be attributed to a circumstellar jet or wind ionised by the embedded (proto)star at this position; -The PV plots of the rotating structure about W3(H 2 O) E is inconclusive on whether the observed rotation is due to a disk-like rotating object, an unresolved binary (or multiple) system, or a combination of the two;   Top: Spectrum of a given pixel (in black) along with the xclass fit for CH 3 CN (12 K − 11 K ) K = 0 − 6 and CH 3 13 CN (12 K − 11 K ) K = 0 − 3 lines (in solid red). The corresponding fit parameters are provided in the panel. Bottom: Spectrum of the same pixel as above (in black) with the xclass fit for CH 3 CN (12 K − 11 K ) K = 4 − 6 and CH 3 13 CN (12 K − 11 K ) K = 0 − 3 (in solid red). The dashed red line corresponds to the predicted spectrum for the CH 3 CN lines that were not used in the fitting process. This shows that the exclusion of the low-K lines in the fitting process allows xclass to provide a better fit for the optically thinner high-K lines. The bright line detected between the K = 4 and 5 components is C 2 H 5 CN.
proto)stars (bottom) at the positions of peak continuum emission as depicted by the two stars. The solid contours correspond to our continuum data in the most extended configuration, starting at 6σ and increasing in steps of 3σ (1σ = 2.5 mJy beam −1 ). The solid vertical line corresponds to the stitching boundary. The dashed line corresponds to Q = 1. Regions outside of the 6σ mm continuum emission contour in the AB configuration are masked out.