© T. Nony et al. 2020

1 Introduction

It is now common knowledge that gas accretion during star formation is not a continuous process (see, e.g., Audard et al. 2014). The temporal variability of accretion has recently been statistically studied for young low-mass stars in the final phase of star formation. Variability timescales from weeks to years have been found in the light curves of young stars observed in the near-IR wavelength range (e.g., Parks et al. 2014; Cody & Hillenbrand 2018). Events of much greater variability on larger timescales have also been found in the light curves of FU Orionis-type stars. Luminosity variability could probe a wide variety of processes with different timescales. When associated with variable accretion rates, this variability could trace instabilities developing within the disk, tidal interactions between the disk and a companion, or even inflowing gas streams arising from the stellar environment (e.g., Zhu et al. 2009; Lodato & Clarke 2004; Alves et al. 2019).

Unfortunately, much less is known about the variability of accretion during the initial phase of star formation, whereas it is during this phase that most of the final stellar mass accumulates from the gas mass reservoir. Jets and outflows develop as protostellar accretion occurs, and they may evacuate the angular momentum of the collapsing protostellar core. The morphology and kinematics of molecular outflows could thus provide fossil records of the accretion history of a protostar (see, e.g., reviews by Arce et al. 2007; Bally 2016). The link between the outflow power and the accretion rate is far from direct, however. It could even be loose due to the still-unknown physical mechanism that is at the origin of jets and outflows, and because of the density and kinematic structure of the protostar’s environment. Despite these caveats, outflow studies have already given very valuable constraints on the protostellar accretion mechanisms. They revealed that accretion rates are stronger at the beginning of the protostellar accretion phase of low-mass stars (e.g., Bontemps et al. 1996) and that they are stronger in high-mass (>8 M_⊙) than low-mass protostars (e.g., Beuther et al. 2002; Duarte-Cabral et al. 2013; Maud et al. 2015). Detailed studies of molecular outflows driven by low-mass stars also suggested that the accretion could be episodic (e.g., Cernicharo & Reipurth 1996; Arce & Goodman 2001; Plunkett et al. 2015; Jhan & Lee 2016). A statistical study of bipolar outflows driven by high-mass protostars of Cygnus X also favored a scenario with intermittent accretion rates (e.g., Duarte-Cabral et al. 2013).

Several theoretical models have been proposed to explain how molecular outflows are formed and driven in the cloud (see the review by Frank et al. 2014). The two most compelling ones invoke either a wind launched from the whole surface of the protostellar disk (disk-wind model by, e.g., Konigl & Pudritz 2000), or a jet + wind initiated at the innermost region of this disk (X-wind model by, e.g., Shu et al. 2000). Whatever the launching mechanism, these winds could easily recollimate along magnetic-field lines and become jet-like. Outflows could then develop from this collimated jet when it crosses the envelope, creates internal shocks and a terminal bow shock, and sweeps the ambientgas. A wide-angle wind would itself drive shells through the surrounding envelope. Given the large variety of length, opening angle, mass, and velocity structure observed for protostellar outflows, both X-wind and disk-wind physical mechanisms may be acting. Each model predicts different relations in diagrams, such as the position-velocity (PV), mass-velocity, and position-magnetic field diagrams, which could thus be used to differentiate between them (e.g., Arce et al. 2007). A couple of models have shown PV diagrams associated with episodic ejection (e.g., Rohde et al. 2019; Vorobyov et al. 2018), but they strongly depend on the structure adopted for the protostellar envelope.

From an observational point of view, young, low-mass protostars referred to as Class 0s have molecular outflows extending up to a few ~0.1 pc, with collimation angles of only a few degrees, and velocities up to 100 km s⁻¹ (Bally 2016). These molecular outflows are often composed of a collimated jet at high velocity, referred to as the “extremely high-velocity” component, plus a shell-like wind tracing gas at lower speeds, referred to as the “low-velocity” or sometimes “standard high-velocity” component. Besides, collimated jets frequently appear as chains of knots interpreted as internal shocks produced byepisodic variations in the protostar mass-loss rate and/or the ejection velocity. Position velocity diagrams of these jets often display velocity spurs referred to as “Hubble laws”, since the gas velocity of these ejection phenomena increases with the distance to the protostar (case of, e.g., the HH211 Class 0, Gueth & Guilloteau 1999). A series of these velocity spurs composes a jagged profile sometimes called “Hubble wedge” (Arce & Goodman 2001). These complex PV diagrams have been investigated in about a dozen low-mass Class 0 protostars so far including L1448-C(N) (Bachiller et al. 1990), HH111 (Cernicharo & Reipurth 1996), L1551 (Bachiller et al. 1994), L1157 (Gueth et al. 1998), IRAS 04166(Tafalla et al. 2004), HH212 (Lee et al. 2006), HH46-47 (Arce et al. 2013), and CARMA-7 in Serpens South (Plunkett et al. 2015).

As for high-mass protostars, current outflow studies often are limited by their angular resolution (Beuther et al. 2007; Maud et al. 2015; Cunningham et al. 2016) or are focused on only a handful of objects (Duarte-Cabral et al. 2013; Cheng et al. 2019). Surveys with high angular resolution in high-mass star-forming regions are thus unavoidable when it comes to characterizing high-mass protostellar outflows and finally constraining the mass-loss or even accretion history of high-mass protostars. Interferometric mosaics at (sub)millimeter wavelengths are currently the best way to image a large and homogeneous sample of outflows driven by protostars forming within the same cloud. W43-MM1 has been imaged with ALMA and represents one of the richest protocluster forming massive stars (Motte et al. 2018a). It is located at 5.5 kpc from the Sun (Zhang et al. 2014) at the tip of the Galactic bar (Nguyen Luong et al. 2011). This 6 pc² cloud ridge has a 2 × 10⁴ M_⊙ mass and qualifies as a “mini-starburst,” because its star formation activity is reminiscent of that of starburst galaxies (SFR ~ 6000 M_⊙ Myr⁻¹, Motte et al. 2003; Louvet et al. 2014). Motte et al. (2018a) detected 131 cores in W43-MM1, and among them 13 high-mass cores. With masses from 16 to 100 M_⊙ within their ~2000 AU diameters, these cores are expected to form high-mass stars in the near future.

Here, we present an ALMA mosaic of the W43-MM1 cloud at ~230 GHz, which reveals a cluster of molecular ouflows driven by high- to intermediate-mass protostars. From observations presented in Sect. 2, we identify 46 molecular outflow lobes and characterize their morphologies (see Sect. 3). In Sect. 4, we discuss the ejection events found in molecular jets and interpret them as reflecting some accretion bursts, and in Sect. 5 we summarize our findings.

2 Observations

Observations were carried out in Cycle 2 between July 2014 and June 2015 (project #2013.1.01365.S), with ALMA 12-m and 7-m (ACA) arrays and baselines ranging from 7.6 to 1045 m. W43-MM1 was imaged with a 78′′ × 53′′ (2.1 pc × 1.4 pc) mosaic composed of 33 fields with the 12-m array and 11 fields with ACA. In the 12-m configuration, the primary beam is 26.7′′ (45.8′′ with ACA), and the maximum detectable scale is ~ 12′′ (~ 21′′ with ACA). Data were reduced in CASA 4.7, applying manual and self-calibration scripts. We ran the CLEAN algorithm with CASA 5.4 using a robust weight of 0.5 and the multiscale (with 0, 1, 3, 9, and 17 beams) option in order to minimize interferometric artefacts due to missing short spacings. Line data were imaged using the merged (12 m + 7 m) configurations, while the continuum was imaged with the 12-m and 7-m arrays separately. The 12-m-only data provide the best sensitivity image chosen to identify compact cores (Motte et al. 2018a), while the 7-m-only map is the best to probe the column density background of these cores. Parameters of the four line spectral windows used in this study are presented in Table 1. The 1.3-mm continuum images themselves have 0.44′′ and 5.3′′ angular resolutions and 1σ rms noise levels of ~0.13 and ~9.1 mJy beam⁻¹ (0.015 and 0.007 K).

Fig. 1 Cluster of molecular outflows in the W43-MM1 massive protocluster. Contours of the integrated blue- and red-shifted emission of the CO(2–1) (a) and SiO(5–4) (b) line wings are overlaid on the 1.3 mm continuum emission of the 7 m array (gray scale, with a 5.3′′ beam and 10 mJy beam−1 rms). (a) The CO line is integrated over 42−85 km s−1 for blue lobes and 111−154 km s−1 for red lobes. Contours are 15, 35 to 235 by steps of 50 in units of σblue = 33 mJy beam−1 km s−1 and σred = 42 mJy beam−1 km s−1. (b) The SiO line is integrated over 43−93 km s−1 for blue lobes and 103−153 km s−1 for red lobes.Contours are 8, 18, 28 to 203 in units of σ = 23 mJy beam−1 km s−1. Ellipses in the lower-left corner represent the angular resolutions of each image, and a scale bar is shown on the top panel to indicate the size in physical units. Boxes outline the zooms-ins of Figs. 3, 4.

3 Results and analysis

Using the CO(2–1) and SiO(5–4) lines¹, we conducted a systematic search for outflows driven by the cores that were detected at 0.44′′ angular resolution (or 2400 AU at 5.5 kpc) in the 1.3-mm continuum image by Motte et al. (2018a).

3.1 Discovery of a rich cluster of outflows

Figures 1a,b provide an overview of the cluster of molecular outflows found in W43-MM1. The blue- and red-shifted wings of the CO(2–1) and SiO(5–4) lines centered at the cloud velocity at rest, V_LSR ~ 98 km s⁻¹, were integrated over the entire velocity range that we observed for the line wings. In the following, we took the same velocity range for all outflows: 42−85 km s⁻¹ and 111−154 km s⁻¹ for CO(2–1); 43−93 km s⁻¹ and 103−153 km s⁻¹ for SiO(5–4) (see Fig. 2). These velocity ranges are different because CO lines present effects of self-absorption and missing short spacings close to their central velocities (see Fig. 2a). The CO(2–1) line also exhibits emission and absorption features associated with low-density clouds located along the line of sight of W43, which were identified by Nguyen Luong et al. (2011) and Carlhoff et al. (2013). These features have a negligible impact on our wide-band integrations, except for the absorption observed at ~80 km s⁻¹, which we excluded by ignoring a handful of channels.

In Figs. 1a,b, we discovered several tens of outflow lobes, whose existence was confirmed by investigating the velocity cubes of both CO(2–1) and SiO(5–4) lines (see Sect. 3.3). Table 2 lists the 27 cores that drive these outflows, gives their core masses, M_core, their velocity at rest, V_LSR, and characterizes each of their associated outflow lobes. Table 2 separates cores according to their location in W43-MM1: 20 are located within its center and seven in its south-western part (Center and SW in Fig. 1b).

A channel map analysis revealed that most of the CO outflows of W43-MM1 protostars present two distinct velocity components: a collimated jet-like structure at high velocity (|ΔV | > 30 km s⁻¹ from the cloud velocity at rest, V_LSR ~ 98 km s⁻¹) plus a broader and more intense outflow at lower velocity. Molecular jets generally exhibit a series of knots that correspond to bumps in the CO(2–1) and SiO(5–4) spectra (see Figs. 2b–d). Similar features have been observed in the spectra of outflows driven by low-mass protostars (e.g., Santiago-Garcia et al. 2009; Lefloch et al. 2015). We therefore integrated the wings of the CO line in two separated ranges, hereinafter called the low- and high-velocity (LV and HV) ranges, recalling the standard and extremely high-velocity ranges used by, for example, Tafalla et al. (2004). For the blue- and red-shifted wings, they correspond to ΔV = −16 ; − 10 km s⁻¹ (resp. 10 ; 21 km s⁻¹) and ΔV = −56 ; − 34 km s⁻¹ (resp. 30 ; 60 km s⁻¹) from the V_LSR. Figures 3 and 4 present the outflows driven by the protostars listed in Table 2 and highlight their ejection direction and length. Isolated contours in Figs. 3 and 4 outline diffuse gas observed along the same line of sight as W43-MM1, but not strictly associated with it. These diffuse structures could lie within the outskirts of the W43 complex (Nguyen Luong et al. 2011).

Fig. 2 Typical CO(2–1) (a and b) and SiO(5–4) (c and d) lines toward molecular outflows, here that of core #8. (a, c): lines averaged over the complete extent of the outflow, thus including both the blue and red lobes. (b, d): lines averaged over one high-velocity knot, here knot R2 of core #8. Line absorption around the velocity at rest, indicated by a dashed line, is due to missing short spacings and optical thickness. The transition between the low- and high-velocity components of the molecular outflows is indicated in panels b and d with a blue dotted line.

3.2 Association with dust cores

The ¹³CS(5–4) line was chosen to estimate the cores’ velocity because it traces high-density gas (critical density of ~ 10⁶ cm⁻³). The ¹³CS emission is strong in the central region of W43-MM1 but weaker in its less dense regions like its southwestern and northern parts. It is clearly detected toward cores or their close surrounding in the central region (see Fig. A.1a), while only half of the cores are detected in the south-western part of W43-MM1 (see Fig. A.1b). We used the velocities of a Gaussian fitted to the ¹³CS(5–4) line integrated over the core, except for cores #1, #2, and #4 for which we used the O¹³ CS(18–17) line. This line traces extremely dense and warm gas, and is thus less confused than ¹³CS(5–4) by the complex environment of the most central part of W43-MM1, while still being strong enough for a fit (Pouteau et al., in prep.). For the cores detected neither in ¹³CS nor in O¹³CS(18–17), we assumed that their velocity at rest would be close to that of their neighboring cores and used 95 km s⁻¹ for undetected cores in the south-westernmost part, and 96 km s⁻¹ for core #31. In W43-MM1, the cores’ velocities at rest, listed in Table 2, range from 95 to 102 km s⁻¹ with a gradient from north-east to south-west and a median of ~97.7 km s⁻¹ over the whole cloud, and ~98.6 km s⁻¹ in the central region (see Fig. A.1). Given that the ¹³CS(5–4) lines were fit without subtracting any emission from the core’s foreground and background, their velocities at rest, V_LSR, are estimated to be uncertain by up to 1 km s⁻¹. These V_LSR values and the gradient are consistent with the ones measured at large scales (28′′ or ~0.7 pc) by Nguyen Luong et al. (2013).

The association between cores and molecular outflows was mostly done using the high-velocity jet component of CO(2–1), which is the most collimated and the easiest to distinguish from the cloud emission (see Figs. 3a–d). However, the low-velocity component was necessary to identify one third of the 46 outflow lobes because their high-velocity component is either missing or too weak. SiO(5–4) was used as a secondary outflow tracer, as in Nony et al. (2018), since it is detected along all but four of the CO outflow lobes, excluding outflows from cores #11, #23, #31 (see Table 2, Fig. 9 and compare Figs. 1a and b, 3, 4 and A.2). The association between cores and outflows is considered reliable for 19 cores and more tentative for eight cores (see Table 2). The latter consist of outflows developing in confused environments, with several lobes overlapping (cores #4, #10, #13, #14, red lobes of cores #11 and #14); outflows only detected in CO (cores #23, #11, red lobe of core #31); and/or outflows only detected at low velocities (ΔV_max ≤ 20 km s⁻¹, cores #4, #13, #31, blue lobes of cores #11 and #36).

3.3 Main characteristics of outflow lobes

Among the main characteristics of outflows, one can estimate their maximum velocities and lengths. For this purpose, we investigated the channel maps of CO(2–1) using a 5σ = 12.5 mJy beam⁻¹ threshold. For each outflow, we measured its maximum velocity projected along the line of sight and relative to the core velocity at rest, ΔV_max = |V_max − V_LSR|, and its maximum length projected onthe plane of the sky, L_max (see Table 2). In a handful of cases, these values are lower limits because their estimates are confused by other outflows or line-of-sight features that dominate the CO wings (see Figs. 3, 4 and 2a). We estimated that these velocity and size estimates are uncertain by up to 2 km s⁻¹ (twice the precision of our velocity measurements) and 0.2′′ (about half our resolution element, 0.005 pc at 5.5 kpc), respectively. Measured velocities range from 11 to 101 km s⁻¹, with a median value of ~47 km s⁻¹, and measured lengths range from 0.02 to 0.4 pc.

Figures 3b–d and 4a,b display the directions of the molecular jets associated with each lobe of the 28 detected outflows. Molecular jets are generally well collimated, keeping a constant width close to our resolution limit, of about 0.04 pc (or 8400 AU) all along their length. The longest molecular jets, ≳0.1 pc without anycorrection for the plane-of-the-sky projection frequently show direction variations, which are outlined by broken segments and arrows in Figs. 3 and 4. Table 2 lists the position angles of molecular jets, defined as the angle of the line passing through the furthest knot of the molecular jet before it undergoes any jet deflection. We thus characterized the closest, and therefore the youngest ejections from the core, connecting at best the blue and red lobes. In most of the cases, this comes down to taking the core center as reference. However a couple of jets are offset by ~0.3′′ with respectto the core center (cores #8 and #9 especially, see Figs. 3b and d). The angle uncertainty relies on the precision of the knots’ location, corresponding to our ability to connect one beam element at a typical distance of ~ 2′′. It is generally 0.5°−1°, but reaches 4° for the shortest lobes.

The majority (18 out of 28) of the molecular outflows we discovered are bipolar (see Table 2). Three of the monopolar outflows could have their blue or red counterparts confused by other outflows (#10 and #14) or disregarded at low velocity (#4). In total, we thus evaluate the number of monopolar outflows to correspond to about 25% of the complete sample of outflows in W43-MM1.

For most of the molecular outflow lobes, the low-velocity (|ΔV |~ 10−20 km s⁻¹) emission presents a wider morphology, but is similar to that of the high-velocity jet (see flows driven by cores #3, #7, #8, #9, and #26 in Figs. 3b–d and 4a).

For the longest outflows, the observed cavity, however tends not to be symmetric around the jet axis (case of cores #7 and #8, see Figs. 3b,c). In marked contrast, core #15 developed outflow lobes with wider-open angle, 30°, cavities, and varying velocities when projected on the line of sight (see Fig. 4a). To give more general trends is, however, difficult because our observation limitations (effect of self-absorption and missing short spacings) preclude the detection of the outflow cavities developing close to the systemic velocity (|ΔV | < 10 km s⁻¹).

Fig. 3 Molecular outflows driven by cores in the central region of the W43-MM1 ridge. Outflows detected at high velocity (a) are shown in detail with their low- and high-velocity components in zooms (b), (c), and (d). Contours of the integrated blue- and red-shifted emission of the CO(2–1) line wings are overlaid on the 1.3-mm continuum emission of the 12-m array (gray scale). For blue lobes, emission is integrated over 42− 64 km s−1 (HV, blue contours) and 82−88 km s−1 (LV, cyan contours). For red lobes, emission is integrated over 108−119 km s−1 (LV, orange contours) and 128−158 km s−1 (HV, red contours). For blue lobes, contours are 7, 15, 30 to 230 by steps of 40 (HV), and 7, 15 to 120 by steps of 15 (LV), in units of σ = 20 mJy beam−1 km s−1. For red lobes, contours are 10, 20 to 160 by steps of 20, in units of σLV, R = 37 mJy beam−1 km s−1 (LV) and 7, 15, 30 to 280 by steps of 50, in units of σ (HV). Green ellipses locate the W43-MM1 cores, arrows (sometimes broken) indicate the direction of their outflows, dashed line outlines the best developed outflow cavities. Knots of the blue lobe of core #8 are labelled B1 to B7 in zoom (b) (the first knot B0 is not visible on this integrated map). Ellipses in the lower corners represent the angular resolutions of the images, and a scale bar is shown in (a) and (b) to indicate the size in physical units (same physical units in c, d and Fig. 4 than in b). Some isolated features, like those outlined by the blue contours below the Z2 box in a, correspond to diffuse gas not associated with the W43-MM1 outflows.

3.4 Cloud environment crossed by outflow lobes

We characterized the cloud density within which each outflow lobe developed by estimating the mean column density, or, equivalently, the visual extinction of the cloud over the complete extent of each outflow lobe. The (filtered) column density along any given line of sight may be derived from I_7m, the 1.3-mm flux density measured in the 7 m-array continuum image in a 5.3′′ -beam, by: $$\begin{array}{ll}{N}_{{\text{H}}_2,\text{f}}\hfill & =\frac{I_{\text{7m}}}{{\Omega }_{\text{beam}}\mu m_{\text{H}}{\kappa }_{\text{1}.\text{3mm}}{B}_{\text{1}.\text{3mm}}(T_{\text{dust}})}\hfill \\ \hfill & =9.4\times {10}^{23}{\text{cm}}^{-2}\left(\frac{I_{\text{7m}}}{\text{Jy\hspace{0.17em}}{\text{beam}}^{-1}}\right),\hfill \end{array}$$(1)

where Ω_beam is the beam solid angle, μ = 2.8 is the mean molecular weight, m_H is the mass of atomic hydrogen, κ_1.3mm is the dust opacity per unit mass column density at 1.3 mm, and B_1.3mm(T_dust) is the Planck function for a dust temperatureT_dust. Because W43-MM1 is a warm and dense cloud, we assumed T_dust = 23 K and κ_1.3mm = 0.01 cm g⁻¹ (see arguments in Motte et al. 2018a). Since, due to interferometer filtering, scales larger than ~21′′ are missing, we added a Av = 130 mag level measured at the outskirt of the W43-MM1 cloud ridge on the Herschel column density image, with 25′′ resolution, by Nguyen Luong et al. (2013).

Equation (1) then leads to the following equation for the visual extinction, A_V : $$A_{\text{V}}=\left(\frac{N_{{\text{H}}_2}}{{10}^{21}{\text{cm}}^{-2}}\right)\approx 130 \text{mag+940} \text{mag}\times \left(\frac{I_{\text{7m}}}{\text{Jy\hspace{0.17em}}{\text{beam}}^{-1}}\right).$$(2)

A_V values span more than one order of magnitude, from 140 to 2600 mag (see Table 2). Since the median core environments should have rather homogeneous temperatures and dust opacities, we estimate the relative uncertainty of our A_V values to be about 30% or less. Similarly, the large-scale filtered cloud structures have A_V values uncertain by up to 30%, A_V = 130 ± 40 mag. As for the absolute values of the visual extinction, they could be uncertain by up to a factor of two due to opacities that are generally poorly constrained.

Fig. 4 Molecular outflows driven by cores of south-western part of W43-MM1. Contours correspond to the integrated blue- and red-shifted emission of the CO(2–1) line shown with the same color code as in Fig. 3. They are overlaid on the 1.3-mm continuum emission of the 12-m array (gray scale). Forblue lobes, contours are 5, 15, 30 by steps of 15 (HV) and 10, 20, 40, 60 (LV), in units of σ = 20 mJy beam−1 km s−1. For red lobes, contours are 5, 15, 25, 35, in units of σLV, R = 37 mJy beam−1 km s−1 (LV) and 5, 15 to 75 by steps of 15, in units of σ (HV). Same physical units and same convention as in Fig. 3 for lines, arrows, and ellipses. Isolated features, like those outlined by the cyan contours at the top of panel a, correspond to diffuse gas not associated with the W43-MM1 outflows.

3.5 Knots in jets

Most of the molecular jets observed both in CO(2–1) and SiO(5–4) integrated maps consist of a continuous emission plus a succession of knots, which are emission peaks aligned along the outflow axis (see Figs. 3b–d, 4a,b and Figs. 9a,b, A.2a–c). In order to validate the velocity coherence of these knots initially identified in space (see, e.g., Figs. 3b–d and 4a,b) and to characterize their velocity and intensity structure, we built Position-Velocity (PV) diagrams along the axis of each jet as defined in Sect. 3.3 and displayed in Figs. 3 and 4. In order to smooth local variations and gain in signal-to-noise level, we averaged the PV diagrams over 0.6′′, about half of the mean width of molecular jets detected in W43-MM1. A smaller averaging width was taken when needed to avoid confusion from neighboring outflows. In Figs. 3 and 4, directions of molecular jets are indicated by arrows that follow the jet deflections when they exist. For the few molecular jets with strong deflection, the PV diagram is the juxtaposition of the ones built along each segment of the jet.

Figures 5a–c, 6a–f and 7a,b show the PV diagrams of the outflows that display several knots and are not too confused by other outflows. For a given core with a bipolar outflow, we chose to give negative position offsets to its red lobe and positive offsets to its blue lobe. The PV diagrams present a very complex structure that varies from core to core and often, for a single core, from its red lobe to its blue lobe. As a general trend, the PV diagram of an outflow lobe is constituted of an intense emission at very low velocity plus several linear features/fingers that develop from these low velocities up to the highest velocities of the outflow. The low-velocity part of the outflow is not properly constrained by present observations as it corresponds to absorbed zones at the cloud velocity at rest (see Sect. 3.1). As for the fingers, most of them have velocities that increase approximately linearly with their positional offsets, and thus their distance from the core. This distance-to-velocity relation, generally qualified as “Hubble law” (see, e.g., review by Arce et al. 2007), is outlined by segments in Figs. 5–7. The best examples of fingers with Hubble law gas distribution can be found in theoutflow lobes of core #8 and core #7 (see Figs. 5a and c). The tips, at high velocity, of these linear features are labelled with “R” or “B” plus their ordered number, because they correspond to the spatial knots observed along the jets (see Figs. 3b,c and 4a,b). We stress that this numbering does not indicate the pairs of (red and blue) knots that would have been simultaneously ejected by the protostar.

In total,we discovered 86 knots in the integrated maps and PV diagrams of CO(2–1), contributing to 38 out of 46 outflow lobes drivenby 22 cores in W43-MM1. The mean number of jet knots per outflow lobe with at least one knot thus is ~2.3. Interestingly enough, outflows that do not display any jet knot are low-velocity (<15 km s⁻¹, for cores #31 and #36) or are confused by the overlap with other outflows (cores #4, #11, and #14) (see Table 2 and Figs. 3b,c).

Table 3 lists each of the detected knots with the distance to their associated protostellar core or the connecting point between red and blue lobes, r_knot, and their velocity relative to the core velocity at rest, ΔV_knot = |V_knot − V_LSR|. Velocities are measured on the high-velocity vertical segments associated to each knot, either at its local maximum or at half the segment when no clear maximum is observed (see Figs. 5–7). The uncertainties of knot velocities and positions were estimated to be ~3 km s⁻¹ and ~0.2′′ (~0.005 pc), respectively. Table 3 also lists the mean velocity offset of knots detected in each lobe, calculated as $\Delta V_{\text{lobe}}={\overline{\Delta V_{\text{knot}}}}^{\text{lobe}}$. The uncertainty on ΔV_lobe is taken to be the dispersion between all the (1 to 8) measurements of V_knot with a minimum value of 3 km s⁻¹. While the dispersion is relatively small, typically ~4 − 5 km s⁻¹ (see, e.g, Figs. 5b or 6a), four lobes have a dispersion larger than 10 km s⁻¹, among which we count the blue lobes of cores #8 and #15. The first one is the longest lobe in the sample and could display velocity wiggling like those found by Lee et al. (2015) or Choi et al. (2017). The second atypical lobe propagates in a complex environment. The velocity offsets measured for the blue and red lobes of the 15 bipolar outflows of Table 3 are similar to within ~15% or ~5 km s⁻¹, with the highest asymmetry measured for the high-mass cores #3, #7, and #16. Over the complete sample of 46 lobes, the blue and red lobes have the same median value, ${\overline{\Delta V_{\text{knot}}}}^{\text{blue\hspace{0.17em}lobe}}\simeq {\overline{\Delta V_{\text{knot}}}}^{\text{red\hspace{0.17em}lobe}}~42 \text{km\hspace{0.17em}}{\text{s}}^{-1}$.

Fig. 5 PV diagrams along molecular outflows of cores #8 (a), #3 (b), and #7 (c). Position offsets, r on the X-axis, are given with respect to the core center or the connecting point between the red and blue lobes. Negative and positive position offsets correspond to the red and blue lobes, respectively. Each figure has the same physical scales on the X-axis and on the Y -axis. Contours are 5, 10, 22 to 67 by 15 steps in units of σCO = 2.5 mJy beam−1. Additional 16σCO contours are drawn in (a) and (c). Vertical and horizontal segments in red locate the detected knots (at rknot in Table 3) and outline their associated velocity range from which Vknot was measured. Segments with arrows in (a) and (b) indicate the knots whose velocity range is limited by absorption features. Fingers, with Hubble-law gas distributions connecting the knots with lower velocity structures, are indicated with continuous black lines, gas layers possibly associated with lateral forward shocks are indicated with dotted red lines. The PV area confused by another outflow is hatched in (b).

Fig. 6 PV diagrams along the molecular outflows of cores #9 (a), #49 (b), #2 (c), #39 (d), #44 (e), and #67 (f). Same convention of position offsets, segments, and lines as in Fig. 5. Contours are 5, 10, 22 to 67 by steps of 15 in (a), 4, 8, 15 to 65 by steps of 10 in (b) and (d–f), 5, 10, 22 to 122 by steps of 25 in (c), in units of σCO =2.5 mJy beam−1. Additional 2σCO contours aredrawn in dotted line in (b) and (f) to highlight weaker emission. PV areas confused by another outflow in (a) and hot core lines in (c) are hatched.

Fig. 7 PV diagrams along the outflows of cores #15 (a) and #26 (b). Same convention of position offsets, segments, and lines as in Fig. 5. Contours are 4, 8, 15 to 65 by steps of 10 in units of σCO =2.5 mJy beam−1. PV areas with confusion between the blue lobe of core #15 and the red lobe of core #26 are hatched.

4 Discussion

With 46 outflow lobes detected, corresponding to 36 bipolar and 10 monopolar outflows, W43-MM1 is one of the richest protocluster found to date in terms of outflows per area. In particular, it is two times richer than the nearby protocluster called NGC 1333 (Plunkett et al. 2013), over a similar spatial area, ~1 pc. Our discovery is due to a combination of cluster characteristics – with a close-packed cluster of cores (Motte et al. 2018a) – suitable angular resolution (2000 AU) to disentangle outflow lobes, and the detailed study in velocity cubes (see Sect. 3.1) allowed by the 1.3 km s⁻¹ resolution. This sample of 46 outflow lobes, ejected by 27 individual cores from a single protocluster at an early stage of evolution, uniquely provides us with the opportunity to make meaningful statistical studies of outflow properties. We hereinafter investigate the influence of the core mass and cloud environment on the main characteristics of outflows (see Sect. 4.1) and discuss the history of ejection events in W43-MM1 (see Sect. 4.2).

4.1 Outflow characteristics in a dynamical environment

The protostars that drive the outflow lobes listed in Table 2 have core masses spanning two orders of magnitude, 1.4 M_⊙ to 102 M_⊙ (Motte et al. 2018a). They could be the precursors of stars with masses ranging from ~0.5 M_⊙ to ~80 M_⊙, assuming a gas-to-star conversion factor of 30−80% (Alves et al. 2007; Bontemps et al. 2010; Louvet et al. 2014). Since the W43-MM1 protostellar population is young (Motte et al. 2003, 2018a), in agreement with the strong outflow collimation found in Sect. 3.3, we expect protostellar outflows to develop in a relatively pristine cloud, devoid of stellar feedback effects. The cloud regions crossed by outflows have different densities, with visual extinctions varying from 140 to 2600 mag (see Table 2).

4.1.1 Outflow maximum length and velocity

Figures 8a,b compare the maximum length of molecular outflows with the mass of the launching cores and the extinction level, which is used as a proxy for the density in which the outflows propagate. While no relation is seen between the outflow length and the core mass in Fig. 8a, a clear anti-correlation of the length of outflow lobes projected on the plane of the sky with the visual extinction appears in Fig. 8b. A linear regression gives the relation $L_{\text{max}}\propto A_{\text{V}}^{\alpha }$ with α =−0.5 ± 0.1, with a Pearson correlation coefficient of −0.70 and a p-value of 7 × 10⁻⁸. This anti-correlation is consistent with the asymmetry in size of the blue and red lobes driven by cores #7 and #8: 0.08 pc vs. 0.18 pc and0.44 pc vs. 0.13 pc, respectively. Their longest outflow lobe indeed propagates away from the cluster, and their shortest lobe toward the inner part of the cluster (see Figs. 3b, c).

We also compared the maximal velocity, ΔV_max, of the molecular jets with both the mass of the driving core and the background cloud extinction (see Figs. A.3a,b). The large scatters of the observed points do not suggest any correlation. However, we remind the reader that the unknown inclination of the individual lobes introduces errors on the length and velocity measurements. As shown for dynamical times in Sect. 4.2.3, this should introduce a scatter when plotting observed length and velocitycharacteristics for the complete sample of 46 outflow lobes. This scatter could hide real physical trends, like those expected between the mass of the launching core and the outflow characteristics, in particular the correlation observed with the outflow momentum flux (e.g., Bontemps et al. 1996; Duarte-Cabral et al. 2013; Maud et al. 2015).

Regardless the correlations or noncorrelations discussed above, they are all obviously limited by our ability to detect protostellar outflows in the central part of W43-MM1. This ridge contains extremely dense, 10⁴ −10⁷ cm⁻³ (see Sect. 3.4, Nguyen Luong et al. 2013; Louvet et al. 2014, 2016), material that should slow down the propagation of molecular outflows. While no global slowing down is observed in Fig. A.3b, the inspection of the PV diagrams of the asymmetric outflows of cores #8 and #7 possibly give some insights (see Figs. 5a and c). For core #8, the velocities of knots, where they remain observable, are not much impacted, with the velocity offsets, ΔV_knot, of R3 and R4 similar to that of R1 (see Fig. 5a), and to the mean velocity of the blue knots, ${\overline{\Delta V_{\text{knot}}}}^{\text{blue}\text{lobe}}$ (see Table 3). In contrast, the series of nested shock structures, seen in the PV diagram of Fig. 5a as a continuous emission plateau over ΔV ≃ 20−40 km s⁻¹ connecting knots, is less extended in the X-axis for the red lobe (containing R1 to R3) than for its blue counterpart (containing B0 to B3): ~0.1 pc vs. ~0.19 pc. It also getsmore intense at low velocity, reaching I_CO = 0.5 Jy beam⁻¹ vs. I_CO = 0.3 Jy beam⁻¹. These last two elements suggest quenching/containment of the outflowing gas by the denser gas of the globally infalling cloud structure inside which outflows develop. A similar trend is observed in Fig. 5c, where the complete red lobe of core #7 is two times less extended than its blue lobe.

Besides the increased pressure expected for an increased density, the complex kinematics of W43-MM1 should be considered. The W43-MM1 ridge is indeed a highly dynamical cloud with 5− 10 km s⁻¹ gradients (Nguyen Luong et al. 2013; Louvet et al. 2016) and inflowing gas at several tens of km s⁻¹ (Louvet 2015). These complex gas dynamics result in a blurring effect along the velocity axis that prevents the identification of molecular outflows with velocities lower than ~10 km s⁻¹. Other processes, such as the global infall that was reported in W43-MM1 by Motte et al. (2005) and Cortes et al. (2010), could also slow down or even disperse the gas contained within protostellar outflows. This global infall is probably neither spherical nor continuous in time, but may rather consist of sporadic accretion flows along filaments, thus following privileged directions (e.g., Schneider et al. 2010; Peretto et al. 2013; Motte et al. 2018b). The most massive protostars, #1 and #4, are the best sites to which this scenario could apply since they are located at the center of gravity of the infalling protocluster (Herpin et al. 2012; Cortes et al. 2010) and they only developed small, low-velocity monopolar (red) outflows (see Fig. 3c). Interestingly, both their missing lobes and the largest gas inflow observed toward the clump hosting cores #1 and #4 are from the blue-shifted part of the spectrum (Pouteau et al., in prep.). Therefore, one can expect this inflow of material to partly quench and disperse the outflowing gas, thus preventing the development of well-organized blue lobes. The general lack of cavity shells developing around each molecular jet could also be explained by the dispersal of these moderate-velocity outflow components by infalling gas. Numerical simulations of protostellar outflows interacting with a dense and dynamical environment are definitively needed in order to go beyond the qualitative hints given above. They could use the shape in PV diagrams of both the high-velocity part of the plateau and the low-velocity arch observed in, for example, Fig. 5a, to distinguish the passive quenching of density from the active dispersion of a complex gas inflow.

Fig. 8 Maximal outflow length observed projected on the plane of the sky, Lmax, vs. the mass of the launching core, Mcore (a) and the visual extinction of the cloud background crossed by the outflow, AV (b). For the longest outflows, measurements at their tip (in parenthesis in Table 2) are used instead of that averaged over their complete extent. Filled and empty squares pinpoint outflow lobes developing in the central and south-western parts of W43-MM1, respectively. Blue and red colors show theblue- and red-shifted lobes, respectively. The anti-correlation found in (b) (best fit line $L_{\text{max}}\propto A_{\text{V}}^{\alpha }$ with α = −0.5 ± 0.1) suggests that molecular outflows are shorter in denser environments.

Fig. 9 Molecular outflows observed in SiO(5–4) compared with their CO(2–1) emission in the central region (a, Z1 field ofFig. 3a) and in the surrounding of core #67 (b). The SiO line is integrated over 43− 93 km s−1 (continuous blue contours) and 103−153 km s−1 (continuous red contours), with contours of 6, 11, 18 to 198 by steps of 15 (a) and 5, 10, 20 (in b), in units of σSiO = 23 mJy beam−1 km s−1. The CO line is integrated over 21−64 km s−1 (dotted blue contours), with contours 7, 15, 30 to 270 by steps of 40 (a) and 5, 12, 19 (b), in units of σCO,B= 32 mJy beam−1 km s−1; and 128−158 km s−1 (dotted red contours), with contours of 7, 15, 30 to 280 by steps of 50 (a) and 5, 15, 30, 45 (b), in units of σCO,R= 22 mJy beam−1 km s−1. Contours are overlaid on the 1.3-mm continuum emission of the 12-m array (gray scale). Green ellipses locate the W43-MM1 cores, arrows (sometimes broken) indicate the direction of their outflows. Knots of the molecular outflow of core #8 arelabelled B1 to B7.

4.1.2 Outflow chemistry through CO versus SiO tracers

Another argument for the large impact of the cloud kinematics on W43-MM1 outflows arises from the comparison of SiO(5–4) and CO(2–1) tracers. The SiO has indeed proven to give insights on the cloud chemistry and, as a result, on the cloud dynamics (e.g., Louvet et al. 2016).

In total, 42 out of the 46 outflow lobes (i.e., 91%) are detected both in CO(2–1) and SiO(5–4) (see Figs. 9a,b and A.2a–c). Such a high detection rate of the molecular outflows in SiO is frequent in high-mass star-formingregions (see, however, Widmann et al. 2016). Low-angular resolution surveys of the SiO(2–1), (3–2) or (5–4) emission lines toward 50–400 high-mass molecular clumps indeed reported detection rates as high as ~60–90% (López-Sepulcre et al. 2011; Csengeri et al. 2016; Li et al. 2019). These rates are in marked contrast with the rare detection of the SiO(5–4) emission line toward low-mass star-forming regions. For instance, less than 30% of the ~20 best-known outflows arising from low-mass protostars were detected by Chandler & Richer (1997) and Gibb et al. (2004). This discrepancy is most likely due to the excitation conditions of the SiO(5–4) emission line, which has a high critical density of (5− 10) × 10⁶ cm⁻³. Such high densities are more easily reached in high-mass than in low-mass star-forming regions.

Beyond the detection rate of the molecular outflows in the SiO(5–4) emission, the important new result here is that the morphologies of SiO outflows are very similar to those of the CO outflows (see Figs. 9a and A.2a–c). The high angular resolution studies by Duarte-Cabral et al. (2014) and Zhang et al. (2015) reported a (less striking) similarity for 10–20 outflow lobes in Cygnus X and IRDC G28.34+0.06P1, respectively.This excellent correlation between CO and SiO probably stems from the peculiar chemistry of SiO in W43. We recall that the most efficient path for the SiO formation is a gas phase reaction between O₂ and Si (Le Picard et al. 2001). Getting Si in the gas phase generally requires high-velocity shocks that erode the core of the dust grain (≥25 km s⁻¹, Schilke et al. 1997; Gusdorf et al. 2008). In W43-MM1, Louvet et al. (2016), however, showed that a wide-spread SiO emission was associated with low-velocity shocks (~10 km s⁻¹). They explained this atypical association by the presence of a fraction of the total Si abundance (1 to 10%) in the form of SiO located in the mantle of the dust grain. Such a hypothesis can be explained if a first event of high-velocity shocks formed the SiO that got depleted onto the mantle of dust grain. This proposed peculiarity would permit that SiO emits over a large range of shock velocities within the outflows of W43-MM1, and not only where high-velocity (>25 km s⁻¹) shocks are present.

In this respect, the case of core #67, the most isolated core of the W43-MM1 protocluster, is enlightening. It drives a bipolar outflow that is nicely traced in CO(2–1), whereas SiO(5–4) is only detected at the terminal bow shocks, plus in a few locations along the red-shifted lobe (see Fig. 9b). Core #67 is located in a background cloud with visual extinction four times lower than the median A_v = 600 mag value. We speculate that along the outflow of core #67, the volume density is not enough to reach the critical density of the SiO(5–4) emission line, except at the terminal bow shocks. The blue lobe of core #8 shows another good example of the impact of the cloud density on the SiO detection in outflows. While the first part of the jet develops in gas column densities typical of the central region, A_v ≃ 400 mag, its second part is characterized by a lower background level, A_v ≃ 150 mag. The first knots, labeled B0 to B2, propagating in high-density gas, are strongly detected in SiO (>20σ_SiO, see Fig. 9a), while the knots further away, labeled B3 to B6, are not (or barely) detected (≤6σ_SiO). Finally, the furthest knot, labelled B7, is itself detected; it probably corresponds to the terminal bow shock of the outflow lobe, which presents a stronger density.

4.1.3 Angle of ejection versus the direction of outflow lobes

A last argument for the influence of cloud kinematics on W43-MM1 molecular outflows comes from the study of outflow propagation directions. As shown in Figs. 3b–d and 4a,b, we observed changes of directions for the longest molecular jets. Most of them are consistent with “wiggling”, periodic oscillation of the atomic and/or molecular jet position (see outflows ejected by cores #8, #3, #9, and #39 in Figs. 3b–d and Fig. 4a). This pattern could be a signature of jet precession, associated with the misalignment of the rotation and magnetic-field axes or orbital motion due to the presence of close multiple protostars (e.g., Terquem et al. 1999).

When the separation of protostellar multiples is large, one could expect multiples to have formed from the fragmentation of the envelope rather than from that of the disk. Their outflow would then have a greater probability of being misaligned, which would facilitate the detection of multiple outflows driven by multiple protostars hosted in a single core. The only clear-cut case of ejections from multiple systems is observed for core #22. It drives a bipolar outflow along the NE-SW direction (#22a in Table 2) and a single red lobe (#22b), almost perpendicular to #22a (see Fig. 3c). This small detection number of multiple outflows (1 out of 27 cores) could be understood if secondary outflows tend to be weaker and/or of lower velocity, or if multiples generally have small separations. The first interpretation would point to uneven mass binaries, which is amiss for a cluster dominated by high-mass protostars (Duchêne & Kraus 2013). The second interpretation recalls results obtained for young Class 0s (Maury et al. 2010) and agrees with the youth of the W43-MM1 protostars, inferred from their low luminosity (Motte et al. 2003, 2018a) and strong outflow collimation (see Sect. 3.3 and Figs. 3 and 4).

In a few cases, long portions of the molecular jets are deflected, a pattern sometimes inconsistent with jet wiggling. In particular, the red lobes of cores #7 and #67 show, at their mid-length, a bend of ~30° and ~10°, respectively (see Figs. 3c and 4b). For core #67, the gradual bend of its red lobe is in marked contrast with the blue counter-jet that stays straight and narrow along its whole length. In the case of core #7, one cannot directly compare the red bent lobe to its blue counterpart, because the latter develops in an environment confused by other outflows. Since these two bends are obvious and occur locally, three interpretations are possible. These outflows could have encountered an obstacle, like, for instance, a dense core, although none were yet detected at the location of these bends. These deflections could also trace density inhomogeneities in the surrounding medium, assuming molecular outflow consists of entrained gas loosely associated with the jet – as gas forming the cavity shell. Otherwise, if outflows, and especially jet-like outflows, are constituted of both ejected and dragged molecular gas, outflow deflection would indicate that outflows propagate in a gas stream with varying velocities. This third interpretation seems to be the most likely since the W43-MM1 cloud consists of several layers of gas, which all host cores and move at several km s⁻¹ relative to each other, and thus relative to some of the cores (Louvet 2015). Similar outflow bends have been observed for a few protostars, including IRAS 4A and IRAS 18162-2048 (Choi 2005; Fernández-López et al. 2013).

Another deflection could also explain the peculiar behavior of the blue lobe of core #8. At high velocities (ΔV > 50 km s⁻¹), its blue jet is not aligned with the core center and points toward core #11, located 0.8′′ away. At lower velocities however, it connects with core #8, following a direction similar to the last ejection(first knot) of the red outflow lobe. This strong global deflection, with an angle of ~40° and the velocity gradient (~3500 AU over 60 km s⁻¹) observed perpendicular to the jet suggest that the blue outflow lobe of core #8 strongly interacted with the surrounding gas, at the location of this bend, ~0.08 pc east of the core. We cannot exclude that similar deflections happened in other molecular outflows (especially those of cores #3 and #9) since they are generally smaller in size and/or more confused by nearby outflows. Further investigations of the cloud structure and kinematics are necessary in order to come to an accurate conclusion regarding the main reason for the outflow deflections observed here.

The three elements presented in Sects. 4.1.1–4.1.3 argue in favor of a large impact of the dynamical environment of molecular outflows on their characteristics. While the maximum velocity of outflows may not strongly be influenced, their length, chemistry, and direction angle are.

4.2 History of outflow ejections and protostellar accretion

Molecular outflows have been widely used to trace the history of protostellar accretion. Here, we measure the dynamical timescale of the 28 protostellar outflows developing in W43-MM1 (Sect. 4.2.1) and use PV diagrams (Sects. 4.2.2–4.2.3) to constrain the history of outflow ejection and protostellar accretion in this region (Sect. 4.2.4).

4.2.1 Dynamical timescales

Determining protostellar ages is far from trivial. Estimations are obtained from various indirect measurements such as the protostar location in envelope mass vs. luminosity diagrams, spectral energy distribution (SED) slopes, chemical ages, or the dynamical timescale of outflows. Here, we estimated the dynamical timescale of outflows using two different methods. The first one classically estimates it as the ratio between the maximum length and maximum velocity of outflows, t_dyn = L_max∕ΔV_max. Since these maximum values, listed in Table 2, are independently measured in the CO cube, length and velocity do not necessarily correspond to the same gas element and, even less so, to the furthest knot corresponding to the oldest ejection event. We therefore adopted as second estimations of t_dyn the maximum value of the dynamical timescales of knots. They were estimated as the ratio of the distance from the knot to the core and the velocity (offset) of the outflow lobe in which it formed, t_knot = r_knot∕ΔV_lobe (see Table 3 and Sect. 3.5). The uncertainty of this dynamical timescale depends on the dispersion of V_knot and on the angular resolution used to identify knots in each PV diagram, here half of the beam (~0.2′′ or ~0.005 pc). We consider that this second method provides the most reliable estimates of the ejection timescale. The two timescale estimates given above are consistent within the uncertainties for 19 out of the 23 cores or multiples (see Fig. 10a). The dynamical timescales of outflows in W43-MM1, before any correction of projection effects, range from 400- to 10⁴ yr, with a median value of ~1500 yr.

To estimate the real values of dynamical timescales, one needs to correct for projection effects, a process which is hazardous when the inclination angle remains unknown (see however Sect. 4.2.3). Measured velocities and lengths are in fact velocities projected along the line of sight, ΔV_max = ΔV_real × cos(i), where i is the angle to the line of sight, and lengths projected on the plane of the sky L_max = L_real × sin(i). The dynamical timescales corrected for the projection angle would then be t_{dyn, corr} = t_dyn∕tan(i). Assuming arandom distribution of inclination angles, the mean value is given by $\overline{i}={\displaystyle {\int }_0^{\pi /2}i}\sin (i)\text{d}i=1\text{rad}={57.3}^{°}$. This leads to dynamical timescales that, once statistically corrected for homogeneous projection effects, range from t_dyn,corr = 280 to 6300 yr.

Figures 10a,b display the dynamical timescale of protostellar outflows as a function of the mass of the launching core and as a function of the visual extinction of the background cloud. While no correlation is observed with the core’s mass (see Fig. 10a), an anti-correlation is found between the dynamical timescale of outflows and the visual extinction of the background they cross (see Fig. 10b). As shown in Figs. 10a,b, the dynamical timescale of outflows in the central part of W43-MM1 (filled symbols) seems on average – but not systematically – a couple of times shorter than the dynamical timescale of protostellar outflows in the south-westernmost part of W43-MM1 (empty symbols). At face value and, assuming that this dynamical time is representative of the complete ejection history and thus of the true protostellar lifetime, this would imply that cores in the central region are two times younger than those in the south-west. This would suggest thatcloud and subsequent star formation would have started in the south-western part of W43-MM1 first. This fits in the framework of recent cloud formation models (e.g., Lee & Hennebelle 2016; Vázquez-Semadeni et al. 2019), where several cloud and star formation events sum up to form the final population of stars. That said, the anti-correlation observed in Fig. 10b is reminiscent of the one found in Fig. 8b between L_max and A_V. We remind the reader that outflows in the central part of W43-MM1 generally arise from massive protostars and develop in dense environments (Motte et al. 2018a, and Table 2). Therefore, the trend observed here could be completely environmental (see Sect. 4.1). In any case, we stress that outflow dynamical timescales should only be considered as rough estimates of the protostellar lifetimes. The first reason is that outflows very likely get missed or dispersed in the cloud, if they do not simply fall below our detection limit. A second important reason is that our angular resolution does not allow us to measure the velocity of the atomic and/or molecular jet directly ejected by the protostar, but only the velocity of the entrained molecular layer that could be a factor of a few to ten times slower.

Fig. 10 Outflow dynamical timescale corrected for a mean projection angle of 57.3°, tdyn, corr, vs. mass of the launching core, Mcore (a) and the visual extinction of the cloud background crossed by the outflow, AV (b). The dynamical timescale of a given core is evaluated in two ways: as the oldest ejection event (triangles) and as the Lmax ∕ΔVmax ratio (squares). Both dynamical timescales are plotted in (a), showing that they are generally consistent within errorbars. A single estimate is used in (b), which corresponds to the oldest ejection event (triangles) when it can be measured and the Lmax ∕ΔVmax ratio otherwise (squares). No significant correlation is observed in (a). The anti-correlation between tdyn and AV in (b), $t_{\text{dyn}}\propto A_{\text{V}}^{\alpha }$ with α = −0.7 ± 0.1 (magenta line) and a Pearson correlation coefficient of −0.67, is reminiscent of that between Lmax and AV in Fig. 8b.

4.2.2 PV diagrams as open books on ejection variability

To go beyond the global characterization of outflow lobes and study the evolution of the ejection process, PV diagrams are among the most powerful tools. The PV diagrams of Figs. 5a–c, 6a–f, and 7a,b display fingers that generally have velocities that increase approximately linearly, with their positional offsets, and thus their distance, from the core. This is nicely illustrated by the handful of linear features found in the outflow lobes of core #8 (connected to knots R2, B1, B2 in Fig. 5a), of core #7 (connected to knots R2, R3, B1 in Fig. 5c), of core #49 (see Fig. 6b), and of core #44 (see Fig. 6e). Interestingly, the slope of these fingers tends to be steeper or even vertical close to the core and shallower afterward (see B1, B2, and B3 in Fig. 5a and R1, R2, and R3 in Fig. 5c). In a ballistic scenario, this behavior could indicate that the molecular ejecta plus the gas layer at its immediate interface, both lying within our spatial resolution unit, are initially all located at short distances from the core with various velocities. Then, after propagation, the ejecta, roughly corresponding to the knot, reaches further distances than the entrained gas layer that now draws a finger in the PV diagram. Similar changes of slopes have been predicted by models of variable ejections (e.g., Rohde et al. 2019). The distance-velocity relations of these fingers are frequently referred to as Hubble laws. Here, we interpret them as the sum of (1) molecular gas ejected by the protostar and ballistically propagating, and (2) ambient gas dragged and accelerated by the jet (see also Arce et al. 2007; Rohde et al. 2019). In our case, fingers with Hubble-law gas distributions represent traces of enhanced ejection events that add to a more continuous average ejection level. At the tip of these fingers observed in PV diagrams (see, particularly, Figs. 5a–c and Figs. 6b,c,f), vertical segments appear that correspond to knots. These knots, also called internal shocks, mark enhanced velocity and density spots of the forward shock where the ejecta and/or its accompanying accelerated gas catches up with slower material from previous ejections (e.g., Jhan & Lee 2016). In some PV diagrams, a structure presents a vertical or strongly declining slope that connects the knot to the V_LSR, reminiscent of the features seen in simulations (e.g., Lee et al. 2001). This gas layer could correspond to the lateral forward shocks.

In our general case, fingers with increasing velocities, knots, and layers with decreasing velocities create triangular-shape structures in the PV diagrams (see Figs. 5a,c, 6a,b,e,f). Several of these triangular shape structures add up and form nested shocks that appear as a continuous emission plateau in the PV diagrams (see, e.g., Figs. 5a,c). At large offsets, most of the low- to intermediate-velocity (< 30 km s⁻¹) emission disappears from PV diagrams, leaving the high-velocity vertical segments relatively isolated (case of B5 to B7 in Fig. 5a). For a couple of knots, their highest velocity appears upstream, drawing a “slow head, fast tail” profile (see B2 and B3 of core #67 in Fig. 6f). Santiago-Garcia et al. (2009) and more recently Jhan & Lee (2016) attributed these “inverse-Hubble laws” to projection effects of internal shocks, in agreement with the pulsed jet simulations of Stone & Norman (1993).

Most of the PV diagrams observed for outflow lobes of W43-MM1 protostars display several knots (see, e.g., Figs. 5a–c and 6f). This is probably the signature of episodic ejection, with variation of the velocity and/or of the mass of the gas ejected, which could be related to episodic accretion. When observed with lower spatial and/or spectral resolutions or lower sensitivity, only the high-velocity part of the PV diagram stands out. These velocity peaks, corresponding to internal shocks, are then the only structures that can clearly be identified. When several velocity peaks appear, PV diagrams present a jagged profile sometimes called a “Hubble wedge” (Arce & Goodman 2001; Stojimirović et al. 2006; Arce et al. 2013; Plunkett et al. 2015), which should resemble a smoothed version of Fig. 5a, for example.

4.2.3 Time span between consecutive ejection events and projection effects

Despite the fact that most of the lobes in the central region seem limited by the environment (see Sect. 4.1), 22 lobes from 14 cores or multiples display two or more knots. For those, we measured the timescale difference between knots, Δt (see Table 3), suggesting it could be an estimate, within the limit of projection effects, of the time between two successive ejection events. The distribution of the 48 Δt values between these 72 knots varies from 290 to 4300 yr, with a median of ~780 yr (before correction for projection effect). Since we used the same velocity to calculate all t_knot values for each lobe, the Δt measurements mostly relate to the distance between consecutive knots, Δr, through the equation Δt = Δr∕ΔV_lobe. Uncertainties of about 200 yr were calculated using those set for ΔV_lobe (see Table 3) and a constant uncertainty of 0.4′′ for Δr. Figure 11 displays the distribution of timescale differences between knots. Our Δt sample should be complete as long as two nearby knots can be separated in the PV diagrams. This will be the case when their relative distance exceeds the spatial resolution of the CO cube, FWHM = 0.013 pc, and their velocity offset exceeds the first quartile from which Δt values are measured in our sample, ΔV_lobe ~ 33 km s⁻¹ (see Sect. 3.3 and Table 3). When corrected for an average projection angle of 57.3°, the median time span between ejection events and the 75% completeness limit of our sample are ~500 and ~250 yr. The histogram of Fig. 11 clearly peaks at ~500 yr, well above the completeness level, suggesting that a typical timescale between two ejecta exists in W43-MM1.

We evaluated the projection effects of Δt measurements by considering two extreme configurations presented in Figs. 12a,b. Firstly, we assumed that ejections in W43-MM1 have the same periodicity throughout the whole cluster and that the observed dispersion is entirely due to projection effects. We remind the reader that measured distances and velocities are projected on the plane of the sky and on the line of sight, respectively: Δr = Δr_corr × sin(i) and ΔV_lobe = ΔV_lobe,corr × cos(i). In order to look for strong projection effects, we plotted the mean distance between knots of a given lobe, ${\overline{\Delta r}}^{\text{lobe}}$, against the mean velocity of knots relative to the V_LSR, ΔV_lobe (see Fig. 12a). We assumed that the most probable inclination angle, 57.3°, is associated to lobes in the median mode of our sample, meaning those that present a median distance between knots of ${\overline{\Delta r}}^{\text{lobe}}=$ 1.5′′ (or 8300 AU) and a median velocity of ΔV_lobe = 45 km s⁻¹. A configuration with a single periodicity implies inclination angles varying from 29° to 81° to cover the full dispersion in our sample. In Fig. 12a, we plotted the ${\overline{\Delta r}}^{\text{lobe}}$ vs. ΔV_lobe linear relations expected for the most probable inclination angle of i = 57.3° and for the first and last quartile angles, i = 41° and i = 75°. Most lobes of our sample lie within this 50% percentile zone, with striking outliers being lobes driven by core #39 (see Fig. 12a). The low-velocity,long outflow lobes (see Figs. 4a and 6d) driven by core #39 could be lying close to the plane of the sky, i ≳ 80°. If, however, the low-velocity outflow is intrinsic to the core, the latter could be more evolved than the typical cores of W43-MM1. On the other hand, the high-velocity, short outflow lobes driven by core #59, and to a lesser extent the blue lobe of core #16 (see Figs. 3c and 4a), argue for them to lie closer to the line of sight than the average, i ≲ 30°. Given that these cores lie in a confused environment, their peculiarity needs to be confirmed to reinforce this interpretation.

In the second configuration, we once again assumed that each lobe is seen with the same inclination angle, 57.3°, and thus that the observed dispersion of the timescale differences between knots, Δt, is caused by an intrinsic dispersion of time variability in outflows. In that case, the timescale difference between two knots, corrected for projection effects, directly is Δt_corr = Δt∕tan(57.3°) and thus has a median value of ~500 yr (see Fig. 11). The distribution of timescale differences clearly peaks at ~500 yr, but ten Δt_corr values are above ~1000 yr. Four out of the five outliers above ~1200 yr correspond to knots of the #39 outflow,which would thus correspond to the most particular outflow of our sample. In this configuration, timescale variations Δt_corr = 400−800 yr can account for the dispersion of half of the sample (see Fig. 12b).

The other methods we used to infer the inclination angle of observed outflows use the morphology of the ejection in integrated images and/or PV diagrams. In our sample, this is most likely the case of cores #15, #3, and #67. The bipolar outflow of core #15 indeed displays a unique morphology in our sample, with an ejection both blue- and red-shifted, both in its eastern and western lobes (see Fig. 4a). The western jet and the western and eastern cavities contain both blue- and red-shifted gas, with the blue-shifted emission more pronounced east of the core, and the red-shifted emission west of it. The outflow cavity, though not centered on the jet direction has a maximum opening angle from its axis of θ_m = 16°. This configuration suggests that the outflow of core #15 is observed in the third geometrical configuration of Cabrit & Bertout (1986), a bipolar lobe seen close to the plane of the sky: i ≥ 90 − θ_m ≃ 74°. As for core #3, its PV diagram presents curved fingers (see Fig. 5b) that resemble the arc-like structures observed by Lee et al. (2015). This study suggests that carved gas layers could come from sideways ejections from an internal bow shock observed with large inclination angles, which means close to the plane of the sky, i ~ 90°. With similar arguments and following the interpretation of Santiago-Garcia et al. (2009) for the IRAS 04166 protostar, core #67, which presents internal shocks with “slow-head, fast-tail” profiles (see Sect. 4.2.2), would be seen with an inclination angle of about 45°.

Since the three indirect methods (statistical, morphology, and PV diagrams arguments) do not converge toward the same objects that should lie on the plane of the sky or on the line of sight, we chose not to individually correct outflow characteristics, but to keep the most probable angle, 57.3°, for homogeneously deprojecting timescales.

Fig. 11 Distribution of timescale differences between successive ejection events. The top and bottom axes show the timescales before and after the correction by a homogeneous inclination angle of 57.3°, Δtcorr, and Δt, respectively. The 75% completeness limit, ~250 yr after deprojection, is represented by a dashed line. Outflows display a variability with a typical timescale, corrected for projection effects, of $\overline{\Delta t_{\text{corr}}}~500$ yr.

Fig. 12 Mean projected distance between knots, ${\overline{\Delta r}}^{\text{lobe}}$, vs. the mean velocity of knots, ΔVlobe, for 21 outflow lobes with at least 2 knots. Error bars represent the ± 1σ dispersion. Measurements and core numbers are indicated in red for outflow lobes discussed below. (a): blue lines locate the measurements for the mean inclination angle of a random distribution of outflows (57.3°) and for the inclination angles of the first and last quartiles (41°, 76°), assuming a constant variability timescale of Δtcorr = 500 yr. (b): blue lines locate the measurements for the median timescale difference between knots, Δtcorr = 500 yr, and the first and last quartiles (Δtcorr = 400 and 800 yr), assuming a constant inclination angle of 57.3°.

4.2.4 Episodic ejection constraining variable accretion

Our study revealed the existence of a typical timescale for the molecular outflow episodicity in W43-MM1, $\overline{\Delta t_{\text{corr}}}~500{\pm }_{100}^{300}$ yr, with a dispersion that can partly be attributed to projection effects (see Figs. 11 and 12 and Sect. 4.2.3). If this timescale is not induced by propagation instabilities along the outflow, but is instead associated with ejection events, it could characterize episodic accretion.

Constraining the origin, within the close environment of the outflowing protostar, of the instabilities creating this episocity requires the comparison with dedicated models. Unfortunately, very few models have so far attempted to link the ejection and accretion processes using their temporal variability. Among the most recent, the model of Vorobyov et al. (2018) describes the accretion instabilities in the disk of low-mass protostars, like the one whose outflow variability was studied by Plunkett et al. (2015). This model revealed luminosity bursts with a bi-modal distribution constituted by one major and one minor mode with 10³ −10⁴ and 10−100 yr timescales, respectively. The first mode would be associated with the accretion of disk cloudlets that remain unique entities, the second one would correspond to cloudlets that undergo tidal torques and fragment into a series of smaller cloudlets. The Vorobyov et al. (2018) model is likely not suitable for describing the intermittent accretion expected for intermediate- to high-mass protostars in W43-MM1. Indeed, this model does not consider irregular feeding of the disk from the envelope and possibly its environment, while intermittent inflows are expected from a dynamical cloud like W43-MM1 (e.g., Louvet et al. 2016). Moreover, as mentioned by Vorobyov et al. (2018), accretion from a more massive disk, like the ones surrounding more massive protostars, would lead to more intense luminosity bursts and possibly longer delays between bursts. We therefore refrain from comparing the episodicity timescale we obtained in W43-MM1 with the timescales found for the minor and major modes of accretion variability modeled by Vorobyov et al. (2018).

Discrete emission peaks have long been observed along the axis of outflows driven by low-mass protostars (see the reviews of, e.g., Arce et al. 2007; Bally 2016). Very few published mm-wave studies have, however, quantified the time span between knots, even taking into account studies done 20 yr ago (Bachiller et al. 1990; Cernicharo & Reipurth 1996; Arce & Goodman 2001) and more recent studies (Arce et al. 2013; Plunkett et al. 2015; Chen et al. 2016). This handful of measurements also has low statistical value since the figures were obtained for a single source, and they are, except for Plunkett et al. (2015), limited by the small number of observed ejecta, typically 3–5 (22 in Plunkett et al. 2015). In marked contrast, the present study measured a statistically robust characteristic timescale using 72 knots, detected in a cluster of outflows driven by low- to high-mass protostars. The timescale measured here, ~500${\pm }_{100}^{300}$ yr, is broadly consistent with the 300−1000 yr range of values obtained by the handful of studies mentioned above. If we limit our comparison to the statistically robust study of Plunkett et al. (2015), our estimated timescale is slightly larger: ~500${\pm }_{100}^{300}$ vs. ~310 ± 150 yr². A word of caution needs to be added regarding the relevance of direct comparisons with the previously published timescales. While some measurements are computed using velocities typical of atomic jets (200–300 km s⁻¹, Arce & Goodman 2001; Arce et al. 2013), our study and others assume the measured, probably lower, velocity of knots (10–50 km s⁻¹, Cernicharo & Reipurth 1996; Plunkett et al. 2015). As for the inclination angles, Cernicharo & Reipurth (1996), Arce et al. (2013), and Chen et al. (2016) had a precise measurement derived from proper motions but others (Bachiller et al. 1990; Arce & Goodman 2001; Plunkett et al. 2015) only roughly estimated it based on strong, and often unverifiable, assumptions. Our study itself uses a sample of 22 lobes to circumvent the bias associated with unknown inclination angles. Interestingly, outflows observed at increasing angular resolutions generally present several spatial – and thus temporal – characteristic scales. L1448 was, for instance, studied with 3000 and 250 AU resolutions by Bachiller et al. (1990) and Hirano et al. (2010), revealing timescales of ~400³ and ~20 yr, which could correspond to several modes of accretion episodicity. Our timescale estimate is also reminiscent of the 100–200 yr timescales found for jet wiggling (e.g., Choi 2001; Lee et al. 2015) but it is not obvious that similar instabilities would be at the origin of both accretion episodicity and jet wiggling. In addition to the absolute value of the measured variability timescales, one can compare the number of molecular knots found in various studies. We discovered more knots along the molecular outflows of W43-MM1 protostars than those reported for low-mass protostars, even though the latter have already been studied for three decades. If ejection episodicity reflects variability in the accretion process, protostellar accretion could be more variable or episodicity is simply easier to detect in W43-MM1 than in nearby clouds.

To conclude, we argue that studying protostellar clusters is an efficient method to statistically constrain the episodicity of the ejection process, and possibly the accretion process. It makes it possible to circumvent the difficulty of measuring inclination angles in very embedded protostars where proper motions are difficult to study. Interpreting the outflow time variability, measured to have a $\overline{\Delta {\text{t}}_{\text{corr}}}~500{\pm }_{100}^{300}$ yr typical timescale, requires dedicated models of the accretion episodicity and its link to outflow variability. We emphasize that complete physical ingredients should be implemented into these models to fully represent the complex environment of massive disk surrounding the intermediate- to high-mass protostars of W43-MM1.

5 Conclusion

We used ALMA to investigate, at high-spatial, 2600 AU, resolution and with CO(2–1) and SiO(5–4) molecular line tracers, the molecular outflows developing in W43-MM1. Our main conclusions can be summarized as follows:

• 1.

  We discovered a rich cluster of 46 outflow lobes driven by 27 protostars (see Fig. 1 and Figs. 3, 4). They are typically ~0.1 pc long, strongly collimated jet-like structures with high-velocity components reaching typical velocities of ΔV ~ 50 km s⁻¹ (see Table 2).

• 2.

  The protostellar cores driving these molecular outflows span a large mass range, from 1 to 100 M_⊙, yet no significant variation of the maximal length and maximal velocity of outflows with the cores’ mass is observed. 18 out of the 28 molecular outflows are bipolar and only one core is associated with two outflows.

• 3.

  We showed a clear anti-correlation between the outflow maximal length and the visual extinction of the crossed background cloud (see Fig. 8b). We propose that the high-density, dynamical cloud of the central region (up to A_V = 3000 mag and ~10 km s⁻¹) limits the propagation of protostellar outflows.

• 4.

  The SiO(5–4) emission of protostellar outflows presents a remarkable morphological and kinematical coincidence with the CO(2–1) emission (see Figs. 9 and A.2). This excellent correspondence, which remains unusual, is explained by the fact that the SiO(5–4) critical density is easily reached in the background crossed by the W43-MM1 outflows (see Sect. 4.1.2).

• 5.

  The jet-like component of molecular outflows consists of a continuous emission plus a series of knots (see Figs. 3, 4). Observed in the CO(2–1) and SiO(5–4) integrated maps as local maxima, knots are typically separated by 1.5′′, or ~8000 AU. These knots appear in PV diagrams at constant distances from the core (see Figs. 5–7) and are interpreted as internal bow shocks (see Sect. 4.2.2). They could be created by the temporal variations of the velocity and/or the mass flux of the protostellar ejection. The large number of knots detected, 86 along 38 outflow lobes, suggests that W43-MM1 protostars undergo accretion bursts.

• 6.

  The detailed study carried out for numerous PV diagrams revealed complex kinematic structures, including fingers with velocities increasing approximately linearly with the distance from the core (see, e.g., Figs. 5a and c). These structures, sometimes referred to as Hubble laws, could be due to CO gas of the envelope dragged and accelerated by the fast atomic and/or molecular jet.

• 7.

  We estimated the timescales associated to each knot and discovered, thanks to the large statistics achieved in W43-MM1, a typical timescale between ejecta. Once corrected with a uniform inclination angle of 57.3°, the ejection episodicity timescale in W43-MM1 is ~500${\pm }_{100}^{300}$ yr. This timescale is consistent with the few variability timescales observed for low-mass protostars. The physical interpretation of this typical timescale, in terms of instabilities in the close environment of the protostars, remains unclear.

Acknowledgements

We are grateful to Paolo Cortes, Catherine Dougados, and Jonathan Ferreira for useful discussions. This paper makes use of the following ALMA data: #2013.1.01365.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This work was supported by the Programme National de Physique Stellaire and Physique et Chimie du Milieu Interstellaire (PNPS and PCMI) of CNRS/INSU (with INC/INP/IN2P3) co-funded by CEA and CNES. This project has received funding from the European Union’s Horizon 2020 research and innovation programme StarFormMapper under grant agreement No 687528.

Appendix A Additional figures

Fig. A.1 Moment 1 – velocity field – map of 13CS(5–4) line in the central region (a) and the south-western region (b). Hatched ellipses locate the cores with detected molecular outflows, black contours at 3 mJy beam−1 (a), and 7 mJy beam−1 (b) outline the continuum 1.3-mm emission of the 12-m array.

Fig. A.2 Molecular outflows observed in SiO(5–4) compared with their CO(2–1) emission in the central region (a and b, Z2 and Z3 fields of Fig. 3a) and in the south-western region (c). The SiO line is integrated over 43− 93 km s−1 (continuous blue contours) and 103−153 km s−1 (continuous red contours), with contours of 6, 11, 18 to 198 by steps of 15 (a and b) and 6, 11, 18, 28 (c), in units of σSiO = 23 mJy beam−1 km s−1. The CO line is integrated over 21−64 km s−1 (dotted blue contours), with contours 7, 15, 30 to 270 by steps of 40 (a and b) and 5, 12, 20, 30 (c), in units of σCO,B= 32 mJy beam−1 km s−1. It is integrated over 128−158 km s−1 (dotted red contours), with contours of 7, 15, 30 to 280 by steps of 50 (a and b) and 5, 15, 30, 45 (c), in units of σCO,R= 22 mJy beam−1 km s−1. Contours are overlaid on the 1.3 mm continuum emission of the 12-m array (gray scale). Green ellipses locate the W43-MM1 cores, arrows (sometimes broken) indicate the direction of their outflows, dashed lines outline the best-developed outflow cavities.

Fig. A.3 Jet maximal velocity, ΔVmax vs. mass of the launching core, Mcore, (a) and the visual extinction of the cloud background crossed by the outflow, AV, (b). Filled and empty squares pinpoint outflow lobes developing in the central and south-western part of W43-MM1, respectively.Blue and red colors show the blue- and red-shifted lobes, respectively. No clear correlation is revealed.

References

All Tables

All Figures

Fig. 1 Cluster of molecular outflows in the W43-MM1 massive protocluster. Contours of the integrated blue- and red-shifted emission of the CO(2–1) (a) and SiO(5–4) (b) line wings are overlaid on the 1.3 mm continuum emission of the 7 m array (gray scale, with a 5.3′′ beam and 10 mJy beam−1 rms). (a) The CO line is integrated over 42−85 km s−1 for blue lobes and 111−154 km s−1 for red lobes. Contours are 15, 35 to 235 by steps of 50 in units of σblue = 33 mJy beam−1 km s−1 and σred = 42 mJy beam−1 km s−1. (b) The SiO line is integrated over 43−93 km s−1 for blue lobes and 103−153 km s−1 for red lobes.Contours are 8, 18, 28 to 203 in units of σ = 23 mJy beam−1 km s−1. Ellipses in the lower-left corner represent the angular resolutions of each image, and a scale bar is shown on the top panel to indicate the size in physical units. Boxes outline the zooms-ins of Figs. 3, 4.

In the text

Fig. 2 Typical CO(2–1) (a and b) and SiO(5–4) (c and d) lines toward molecular outflows, here that of core #8. (a, c): lines averaged over the complete extent of the outflow, thus including both the blue and red lobes. (b, d): lines averaged over one high-velocity knot, here knot R2 of core #8. Line absorption around the velocity at rest, indicated by a dashed line, is due to missing short spacings and optical thickness. The transition between the low- and high-velocity components of the molecular outflows is indicated in panels b and d with a blue dotted line.

In the text

Fig. 3 Molecular outflows driven by cores in the central region of the W43-MM1 ridge. Outflows detected at high velocity (a) are shown in detail with their low- and high-velocity components in zooms (b), (c), and (d). Contours of the integrated blue- and red-shifted emission of the CO(2–1) line wings are overlaid on the 1.3-mm continuum emission of the 12-m array (gray scale). For blue lobes, emission is integrated over 42− 64 km s−1 (HV, blue contours) and 82−88 km s−1 (LV, cyan contours). For red lobes, emission is integrated over 108−119 km s−1 (LV, orange contours) and 128−158 km s−1 (HV, red contours). For blue lobes, contours are 7, 15, 30 to 230 by steps of 40 (HV), and 7, 15 to 120 by steps of 15 (LV), in units of σ = 20 mJy beam−1 km s−1. For red lobes, contours are 10, 20 to 160 by steps of 20, in units of σLV, R = 37 mJy beam−1 km s−1 (LV) and 7, 15, 30 to 280 by steps of 50, in units of σ (HV). Green ellipses locate the W43-MM1 cores, arrows (sometimes broken) indicate the direction of their outflows, dashed line outlines the best developed outflow cavities. Knots of the blue lobe of core #8 are labelled B1 to B7 in zoom (b) (the first knot B0 is not visible on this integrated map). Ellipses in the lower corners represent the angular resolutions of the images, and a scale bar is shown in (a) and (b) to indicate the size in physical units (same physical units in c, d and Fig. 4 than in b). Some isolated features, like those outlined by the blue contours below the Z2 box in a, correspond to diffuse gas not associated with the W43-MM1 outflows.

In the text

Fig. 4 Molecular outflows driven by cores of south-western part of W43-MM1. Contours correspond to the integrated blue- and red-shifted emission of the CO(2–1) line shown with the same color code as in Fig. 3. They are overlaid on the 1.3-mm continuum emission of the 12-m array (gray scale). Forblue lobes, contours are 5, 15, 30 by steps of 15 (HV) and 10, 20, 40, 60 (LV), in units of σ = 20 mJy beam−1 km s−1. For red lobes, contours are 5, 15, 25, 35, in units of σLV, R = 37 mJy beam−1 km s−1 (LV) and 5, 15 to 75 by steps of 15, in units of σ (HV). Same physical units and same convention as in Fig. 3 for lines, arrows, and ellipses. Isolated features, like those outlined by the cyan contours at the top of panel a, correspond to diffuse gas not associated with the W43-MM1 outflows.

In the text

Fig. 5 PV diagrams along molecular outflows of cores #8 (a), #3 (b), and #7 (c). Position offsets, r on the X-axis, are given with respect to the core center or the connecting point between the red and blue lobes. Negative and positive position offsets correspond to the red and blue lobes, respectively. Each figure has the same physical scales on the X-axis and on the Y -axis. Contours are 5, 10, 22 to 67 by 15 steps in units of σCO = 2.5 mJy beam−1. Additional 16σCO contours are drawn in (a) and (c). Vertical and horizontal segments in red locate the detected knots (at rknot in Table 3) and outline their associated velocity range from which Vknot was measured. Segments with arrows in (a) and (b) indicate the knots whose velocity range is limited by absorption features. Fingers, with Hubble-law gas distributions connecting the knots with lower velocity structures, are indicated with continuous black lines, gas layers possibly associated with lateral forward shocks are indicated with dotted red lines. The PV area confused by another outflow is hatched in (b).

In the text

Fig. 6 PV diagrams along the molecular outflows of cores #9 (a), #49 (b), #2 (c), #39 (d), #44 (e), and #67 (f). Same convention of position offsets, segments, and lines as in Fig. 5. Contours are 5, 10, 22 to 67 by steps of 15 in (a), 4, 8, 15 to 65 by steps of 10 in (b) and (d–f), 5, 10, 22 to 122 by steps of 25 in (c), in units of σCO =2.5 mJy beam−1. Additional 2σCO contours aredrawn in dotted line in (b) and (f) to highlight weaker emission. PV areas confused by another outflow in (a) and hot core lines in (c) are hatched.

In the text

	Fig. 7 PV diagrams along the outflows of cores #15 (a) and #26 (b). Same convention of position offsets, segments, and lines as in Fig. 5. Contours are 4, 8, 15 to 65 by steps of 10 in units of σCO =2.5 mJy beam−1. PV areas with confusion between the blue lobe of core #15 and the red lobe of core #26 are hatched.
In the text

Fig. 8 Maximal outflow length observed projected on the plane of the sky, Lmax, vs. the mass of the launching core, Mcore (a) and the visual extinction of the cloud background crossed by the outflow, AV (b). For the longest outflows, measurements at their tip (in parenthesis in Table 2) are used instead of that averaged over their complete extent. Filled and empty squares pinpoint outflow lobes developing in the central and south-western parts of W43-MM1, respectively. Blue and red colors show theblue- and red-shifted lobes, respectively. The anti-correlation found in (b) (best fit line $L_{\text{max}}\propto A_{\text{V}}^{\alpha }$ with α = −0.5 ± 0.1) suggests that molecular outflows are shorter in denser environments.

In the text

Fig. 9 Molecular outflows observed in SiO(5–4) compared with their CO(2–1) emission in the central region (a, Z1 field ofFig. 3a) and in the surrounding of core #67 (b). The SiO line is integrated over 43− 93 km s−1 (continuous blue contours) and 103−153 km s−1 (continuous red contours), with contours of 6, 11, 18 to 198 by steps of 15 (a) and 5, 10, 20 (in b), in units of σSiO = 23 mJy beam−1 km s−1. The CO line is integrated over 21−64 km s−1 (dotted blue contours), with contours 7, 15, 30 to 270 by steps of 40 (a) and 5, 12, 19 (b), in units of σCO,B= 32 mJy beam−1 km s−1; and 128−158 km s−1 (dotted red contours), with contours of 7, 15, 30 to 280 by steps of 50 (a) and 5, 15, 30, 45 (b), in units of σCO,R= 22 mJy beam−1 km s−1. Contours are overlaid on the 1.3-mm continuum emission of the 12-m array (gray scale). Green ellipses locate the W43-MM1 cores, arrows (sometimes broken) indicate the direction of their outflows. Knots of the molecular outflow of core #8 arelabelled B1 to B7.

In the text

Fig. 10 Outflow dynamical timescale corrected for a mean projection angle of 57.3°, tdyn, corr, vs. mass of the launching core, Mcore (a) and the visual extinction of the cloud background crossed by the outflow, AV (b). The dynamical timescale of a given core is evaluated in two ways: as the oldest ejection event (triangles) and as the Lmax ∕ΔVmax ratio (squares). Both dynamical timescales are plotted in (a), showing that they are generally consistent within errorbars. A single estimate is used in (b), which corresponds to the oldest ejection event (triangles) when it can be measured and the Lmax ∕ΔVmax ratio otherwise (squares). No significant correlation is observed in (a). The anti-correlation between tdyn and AV in (b), $t_{\text{dyn}}\propto A_{\text{V}}^{\alpha }$ with α = −0.7 ± 0.1 (magenta line) and a Pearson correlation coefficient of −0.67, is reminiscent of that between Lmax and AV in Fig. 8b.

In the text

Fig. 11 Distribution of timescale differences between successive ejection events. The top and bottom axes show the timescales before and after the correction by a homogeneous inclination angle of 57.3°, Δtcorr, and Δt, respectively. The 75% completeness limit, ~250 yr after deprojection, is represented by a dashed line. Outflows display a variability with a typical timescale, corrected for projection effects, of $\overline{\Delta t_{\text{corr}}}~500$ yr.

In the text

Fig. 12 Mean projected distance between knots, ${\overline{\Delta r}}^{\text{lobe}}$, vs. the mean velocity of knots, ΔVlobe, for 21 outflow lobes with at least 2 knots. Error bars represent the ± 1σ dispersion. Measurements and core numbers are indicated in red for outflow lobes discussed below. (a): blue lines locate the measurements for the mean inclination angle of a random distribution of outflows (57.3°) and for the inclination angles of the first and last quartiles (41°, 76°), assuming a constant variability timescale of Δtcorr = 500 yr. (b): blue lines locate the measurements for the median timescale difference between knots, Δtcorr = 500 yr, and the first and last quartiles (Δtcorr = 400 and 800 yr), assuming a constant inclination angle of 57.3°.

In the text

	Fig. A.1 Moment 1 – velocity field – map of 13CS(5–4) line in the central region (a) and the south-western region (b). Hatched ellipses locate the cores with detected molecular outflows, black contours at 3 mJy beam−1 (a), and 7 mJy beam−1 (b) outline the continuum 1.3-mm emission of the 12-m array.
In the text

Fig. A.2 Molecular outflows observed in SiO(5–4) compared with their CO(2–1) emission in the central region (a and b, Z2 and Z3 fields of Fig. 3a) and in the south-western region (c). The SiO line is integrated over 43− 93 km s−1 (continuous blue contours) and 103−153 km s−1 (continuous red contours), with contours of 6, 11, 18 to 198 by steps of 15 (a and b) and 6, 11, 18, 28 (c), in units of σSiO = 23 mJy beam−1 km s−1. The CO line is integrated over 21−64 km s−1 (dotted blue contours), with contours 7, 15, 30 to 270 by steps of 40 (a and b) and 5, 12, 20, 30 (c), in units of σCO,B= 32 mJy beam−1 km s−1. It is integrated over 128−158 km s−1 (dotted red contours), with contours of 7, 15, 30 to 280 by steps of 50 (a and b) and 5, 15, 30, 45 (c), in units of σCO,R= 22 mJy beam−1 km s−1. Contours are overlaid on the 1.3 mm continuum emission of the 12-m array (gray scale). Green ellipses locate the W43-MM1 cores, arrows (sometimes broken) indicate the direction of their outflows, dashed lines outline the best-developed outflow cavities.

In the text

	Fig. A.3 Jet maximal velocity, ΔVmax vs. mass of the launching core, Mcore, (a) and the visual extinction of the cloud background crossed by the outflow, AV, (b). Filled and empty squares pinpoint outflow lobes developing in the central and south-western part of W43-MM1, respectively.Blue and red colors show the blue- and red-shifted lobes, respectively. No clear correlation is revealed.
In the text