Episodic accretion constrained by a rich cluster of outflows

Context. The accretion history of protostars remains widely mysterious, even though it represents one of the best ways to understand the protostellar collapse that leads to the formation of stars. Aims. Molecular outﬂows, which are easier to detect than the direct accretion onto the prostellar embryo, are here used to characterize the protostellar accretion phase in W43-MM1. Methods. The W43-MM1 protocluster hosts a sufﬁcient number of protostars to statistically investigate molecular outﬂows in a single, homogeneous region. We used the CO(2–1) and SiO(5–4) line datacubes, taken as part of an ALMA mosaic with a 2000 AU resolution, to search for protostellar outﬂows, evaluate the inﬂuence that the environment has on these outﬂows’ characteristics and put constraints on outﬂow variability in W43-MM1. Results. We discovered a rich cluster of 46 outﬂow lobes, driven by 27 protostars with masses of 1 − 100 M (cid:12) . The complex environment inside which these outﬂow lobes develop has a deﬁnite inﬂuence on their length, limiting the validity of using outﬂows’ dynamical timescale as a proxy of the ejection timescale in clouds with high dynamics and varying conditions. We performed a detailed study of Position–Velocity diagrams of outﬂows that revealed clear events of episodic ejection. The time variability of W43-MM1 outﬂows is a general trend and is more generally observed than in nearby, low- to intermediate-mass star-forming regions. The typical timescale found between two ejecta, ∼ 500 yr, is consistent with that found in nearby protostars. Conclusions. If ejection episodicity reﬂects variability in the accretion process, either protostellar accretion is more variable, or episodicity is easier to detect in high-mass star-forming regions than in nearby clouds. The timescale found between accretion events could result from instabilities associated with bursts of inﬂowing gas arising from the close dynamical environment of high-mass star-forming cores.


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 The CO(2-1) data cube is only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http:// cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/636/A38 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 highmass (>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 ambient gas. 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 −1 (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 by episodic 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 , 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 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 Notes. (a) 1σ rms in units of mJy beam −1 . 1 mJy beam −1 corresponds to 0.12 K at 230 GHz.
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 2 cloud ridge has a 2 × 10 4 M mass and qualifies as a "mini-starburst," because its star formation activity is reminiscent of that of starburst galaxies (SFR ∼ 6000 M Myr −1 , 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.

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-monly 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 −1 (0.015 and 0.007 K).

Results and analysis
Using the CO(2-1) and SiO(5-4) lines 1 , 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). wings of the CO(2-1) and SiO(5-4) lines centered at the cloud velocity at rest, V LSR ∼ 98 km s −1 , 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 −1 and 111−154 km s −1 for CO(2-1); 43−93 km s −1 and 103−153 km s −1 for SiO(5-4) (see Fig. 2). These velocity ranges are different because CO lines present effects of selfabsorption 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 A&A 636, A38 (2020) 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.

Discovery of a rich cluster of outflows
features have a negligible impact on our wide-band integrations, except for the absorption observed at ∼80 km s −1 , 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 −1 from the cloud velocity at rest, V LSR ∼ 98 km s −1 ) 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 −1 (resp. 10 ; 21 km s −1 ) and ∆V = −56 ; −34 km s −1 (resp. 30 ; 60 km s −1 ) 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).

Association with dust cores
The 13 CS(5-4) line was chosen to estimate the cores' velocity because it traces high-density gas (critical density of ∼10 6 cm −3 ). The 13 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 13 CS(5-4) line integrated over the core, except for cores #1, #2, and #4 for which we used the O 13 CS(18-17) line. This line traces extremely dense and warm gas, and is thus less confused than 13 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 13 CS nor in O 13 CS(18-17), we assumed that their velocity at rest would be close to that of their neighboring cores and used 95 km s −1 for undetected cores in the south-westernmost part, and 96 km s −1 for core #31. In W43-MM1, the cores' velocities at rest, listed in Table 2, range from 95 to 102 km s −1 with a gradient from north-east to south-west and a median of ∼97.7 km s −1 over the whole cloud, and ∼98.6 km s −1 in the central region (see Fig. A.1). Given that the 13 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 −1 . These V LSR values and the gradient are consistent with the ones measured at large scales (28 or ∼0.7 pc) by Nguyen .
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 −1 , cores #4, #13, #31, blue lobes of cores #11 and #36).

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 −1 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 on the 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    angle of the molecular jet, measured from West to North. A star marker indicates molecular jets that are offset by more than 0.2 from their core center. (e) Visual extinction measured over the complete extent of the outflow lobe, uncertain by 30%. A second measure at the tip of the molecular jet is given in parenthesis for lobes covered by more than 3 pixels in the 7 m array continuum map (i.e., with L max ≥ 9.6 × 10 −2 pc) and when it differs by more than 20%. ( f ) Less reliable outflow lobes are pinpointed by mentioning their confused environments (denoted with "Conf") or their nondetection in SiO(5-4) but in CO(2-1) (denoted with "CO"). "(R)" means that the note applies to the red lobe only.
to 2 km s −1 (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 −1 , with a median value of ∼47 km s −1 , 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 any correction 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 respect to 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 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 −1 ). A&A 636, A38 (2020) . 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.

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 A38, page 6 of 22 T. Nony et al.: Episodic accretion constrained by a rich cluster of outflows  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.

-beam, by:
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 temperature T dust . Because W43-MM1 is a warm and dense cloud, we assumed T dust = 23 K and κ 1.3mm = 0.01 cm g −1 (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 .
Equation (1) then leads to the following equation for the visual extinction, A V : 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.

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 4a,b and Figs. 9a,b,. 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. , 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 r knot in Table 3) and outline their associated velocity range from which V knot 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).
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 the outflow 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 A38, page 9 of 22 A&A 636, A38 (2020) 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 driven by 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 −1 , 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 −1 and ∼0.2 (∼0.005 pc), respectively. Table 3 also lists the mean velocity offset of knots detected in each lobe, calculated as ∆V lobe = ∆V knot 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 −1 . While the dispersion is relatively small, typically ∼4−5 km s −1 (see, e.g, Figs. 5b or 6a), four lobes have a dispersion larger than 10 km s −1 , 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

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 (

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(Motte et al. , 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).

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.    Table 2) are used instead of that averaged over their complete extent. Filled and empty squares pinpoint outflow lobes developing in the central and southwestern parts of W43-MM1, respectively. Blue and red colors show the blue-and red-shifted lobes, respectively. The anti-correlation found in (b) (best fit line L max ∝ A α V with α = −0.5 ± 0.1) suggests that molecular outflows are shorter in denser environments.
correlation coefficient of −0.70 and a p-value of 7 × 10 −8 . 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 and 0.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 velocity characteristics 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 4 −10 7 cm −3 (see Sect. 3.4, Nguyen Luong et al. 2013;Louvet et al. 2014Louvet et al. , 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, ∆V knot blue 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 −1 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 gets more intense at low velocity, reaching I CO = 0.5 Jy beam −1 vs. I CO = 0.3 Jy beam −1 . 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 −1 gradients Louvet et al. 2016) and inflowing gas at several tens of km s −1 (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 −1 . 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 A38, page 12 of 22 T. Nony et al.: Episodic accretion constrained by a rich cluster of outflows 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.

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,. Such a high detection rate of the molecular outflows in SiO is frequent in high-mass star-forming regions (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 6 cm −3 . 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 2 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 −1 , 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 −1 ). 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 highvelocity 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 A&A 636, A38 (2020) of W43-MM1, and not only where high-velocity (>25 km s −1 ) 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.

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(Motte et al. , 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 midlength, 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 −1 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 −1 ), 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 −1 ) 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.

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

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 Fig. 10. Outflow dynamical timescale corrected for a mean projection angle of 57.3 • , t dyn, corr , vs. mass of the launching core, M core (a) and the visual extinction of the cloud background crossed by the outflow, A V (b). The dynamical timescale of a given core is evaluated in two ways: as the oldest ejection event (triangles) and as the L max /∆V max 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 L max /∆V max ratio otherwise (squares). No significant correlation is observed in (a). The anti-correlation between t dyn and A V in (b), t dyn ∝ A α V with α = −0.7 ± 0.1 (magenta line) and a Pearson correlation coefficient of −0.67, is reminiscent of that between L max and A V in Fig. 8b. 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 4 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 a random distribution of inclination angles, the mean value is given by i = π/2 0 i sin (i) di = 1 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 that cloud 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.

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

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 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 • , ∆t corr, 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 ∆t corr ∼ 500 yr.
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 −1 (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, ∆r 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 A38, page 16 of 22 T. Nony et al.: Episodic accretion constrained by a rich cluster of outflows median distance between knots of ∆r lobe = 1.5 (or 8300 AU) and a median velocity of ∆V lobe = 45 km s −1 . 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 ∆r 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.

Episodic ejection constraining variable accretion
Our study revealed the existence of a typical timescale for the molecular outflow episodicity in W43-MM1, ∆t corr ∼ 500± 300 100 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 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 3 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 ∆t corr ∼ 500± 300 100 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.

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 −1 (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 −1 ) 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± 300 100 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.