© ESO 2021

1 Introduction

In the past few decades, several studies have been performed to investigate the properties of fragmentation in molecular clouds. These works have concentrated on the study of the physical processes that control fragmentation such as turbulence (Zhang et al. 2009, 2015; Pillai et al. 2011; Wang et al. 2011, 2014; Lu et al. 2015), disk fragmentation (e.g., Adams et al. 1989; Kratter et al. 2010; Vorobyov & Basu 2010), rotational fragmentation (e.g., Lim et al. 2016; Arreaga-García 2017), stellar feedback (e.g., Krumholz et al. 2007; Peters et al. 2010; Myers et al. 2013), and magnetic fields (e.g., Boss 2004; Vázquez-Semadeni et al. 2005; Commerçon et al. 2011a; Hennebelle et al. 2011). On the other hand, a number of studies have found that the fragmentation of dusty cores seems to be related to the density of the cores (e.g., Gutermuth et al. 2011; Lombardi et al. 2013; Ragan et al. 2013; Imara 2015; Lee et al. 2015; Mairs et al. 2016; Nguyen-Luong et al. 2016; Liu et al. 2016; Pokhrel et al. 2016, 2018; Hacar et al. 2017; Mercimek et al. 2017; Alfaro & Román-Zúñiga 2018; Li et al. 2019; Orkisz et al. 2019; Sanhueza et al. 2019; Sokol et al. 2019; Svoboda et al. 2019; Zhang et al. 2019), indicating that the fragmentation of cores is consistent with a thermal Jeans process in these cases.

Among the aforementioned processes potentially controlling fragmentation, both the magnetic field and stellar feedback have been poorly studied so far from an observational point of view. While for the magnetic field there are very recent works aimed at studying its effects in samples of massive dense cores where both fragmentation and magnetic field properties have been observed (e.g., Añez-López et al. 2020; Palau et al. 2021), the stellar feedback effects have not been broadly studied yet.

In addition, the concept of hierarchical fragmentation (Hoyle 1953) has been explored in several recent observational studies toward low- and high-mass star-forming regions (Teixeira et al. 2006, 2007, 2016; Zhang et al. 2009; Wang et al. 2011, 2014; Hacar et al. 2013; Pineda & Teixeira 2013; Kainulainen et al. 2017; Pokhrel et al. 2018), suggesting that hierarchical fragmentation is a common process (Wang et al. 2011; Zhang et al. 2015; Teixeira et al. 2016; Kainulainen et al. 2017). Therefore, it is very important to study the processes that control fragmentation in more star-forming regions, paying special attention to the massive ones, which are typically associated with complex structures and feedback processes. However, given the large distances of massive star-forming regions, the spatial scales of 200 AU (~0.001 pc) have been poorly explored to date.

The W51 IRS2 region (Ginsburg et al. 2017) harbors a well-known H_II region with a luminosity of ~3 × 10⁶ L_⊙ (Erickson & Tokunaga 1980), and it is located at 5128(±1867) pc (Reid et al. 2019). In the W51 IRS2 region, there are two highly obscured massive sources called W51 North and W51 d2. W51 d2 was reported as an ultracompact H_II (UCH_II) region (Gaume et al. 1993; Lacy et al. 2007), and W51 North, located to the east of W51 d2, is an active massive star-forming region (Zhang et al. 1998; Sollins et al. 2004; Zapata et al. 2008, 2009, 2010). W51 North shows strong emission at millimeter and submillimeter wavelengths, very faint centimeter free-free emission (Gaume et al. 1993; Zhang et al. 1998; Eisner et al. 2002; Zapata et al. 2009), and no mid-infrared counterparts (Kraemer et al. 2001; Okamoto et al. 2001; Lacy et al. 2007). In addition, very strong OH, SiO, and H₂ O maser emissions have been detected (Schneps et al. 1981; Gaume & Mutel 1987; Morita et al. 1992). This emission is associated with the formation of high-mass stars. Fortunately, W51 IRS2 was already been observed by the Atacama Large Millimeter and Submiilimeter Array (ALMA) at different wavelengths with several angular resolutions (the best angular resolution of 0.028′′ corresponds to 144 AU). The high angular resolution capability of ALMA gives us a chance to reveal the small-scale properties of this molecular cloud.

The paper has been organized as follows. Section 2 introduces the archival data and the reduction of all observations. Section 3 presents our results directly derived from the observations. Sections 4 and 5 provide detailed analyses and a discussion on the fragmentation in W51 IRS2. The main findings are summarized in Sect. 6.

2 Observational data

In this work, we used three ALMA datasets of W51 IRS2 from the following projects: 2017.1.00293.S, 2015.1.01596.S, and 2013.1.00994.S. Using the pipeline in CASA (McMullin et al. 2007) (Version 4.7), we obtained the visibilities to perform self-calibration and generate the final continuum images. The continuum images of calibrated visibilities were produced with the task TCLEAN of the CASA package. During the imaging process, for all datasets, a “briggs” weighting ROBUST parameter of 0.5 was set to balance the sensitivity and angular resolution.

2.1 ALMA 3 mm data

Band 3 (3 mm) observations were carried out with ALMA from 2017 October 6 to 15 as part of Cycle 5 program 2017.1.00293.S. The project was carried out in five executions using between 48 and 51 antennas with a diameter of 12 m. The baselines range from 33 to 15983 m (9.389–5.33 × 10³ kλ), which corresponds to the largest recoverable angular scale of 9.7′′ (following Eq. (A.5) of Palau et al. 2010). The observing field of W51 IRS2 was centered at the following sky position: RA(J2000) = 19^h 23^m40^s, Dec (J2000) = +14°31′05′′. Strong maser lines were included in the 3 mm data (Ginsburg & Goddi 2019). Additionally, W51 North is a well-known chemically rich hot core. In order to avoid line contamination, the 3 mm continuum image was only generated using the line-free channels. The channel selection was performed in the uv domain, and the line-free channels were discerned in the four spectral windows, which have a bandwidth of 1.875 GHz, and they are centered at rest frequencies of 86.919 (spw0), 87.828 (spw1), 98.229 (spw2), and 99.978 GHz (spw3). The quasar J2148+0657 was adopted as a flux and bandpass calibrator, and the phase gain calibrator was J1922+1530.

We performed six rounds of phase self-calibration and one round of amplitude self-calibration to visibilities. The phase-only self-calibrations decreased the root mean square (rms) noise of continuum image from 0.5 to 0.15 mJy beam⁻¹. The amplitude self-calibration did not introduce any systematic error, yielding a final continuum image with a rms noise of 0.08 mJy beam⁻¹ and a synthesized beam of 0.078′′ × 0.051′′ (position angle: –45°).

2.2 ALMA 1.3 mm data

2.2.1 High angular resolution

Band 6 (1.3 mm) high angular resolution observations were carried out with ALMA in 2015 October 30 and 31 as part of Cycle 3 program 2015.1.01596.S. The project was performed in two executions using 42 antennas with a diameter of 12 m. The baselines range from 40 to 16161 m (51.563–1.27 × 10⁴ kλ), corresponding to the largest recoverable angular scale of 1.6′′. The phase center of the observations is at the sky position RA(J2000) = 19^h23^m40^s.05, Dec (J2000) = +14°31′05′′.50. The continuum image is obtained by averaging the line-free channels from ten spectral windows. There are two spectral windows with a 1.875 GHz bandwidth, and the remaining eight spectral windows have a 234 MHz bandwidth. The rest frequencies are 217.840 (spw0), 221.310 (spw1), 220.650 (spw2), 219.560 (spw3), 220.400 (spw4), 233.296 (spw5), 232.164 (spw6), 232.861 (spw7), 231.901 (spw8), and 235 GHz (spw9). The flux, bandpass, and phase calibrators were J2148+0657, J1751+0939, and J1920+1540, respectively.

We performed five rounds of phase self-calibration and one round of amplitude self-calibration to visibilities. The phase-only self-calibrations decreased the rms noise of the 1.3 mm high-resolution continuum image from 0.2 to 0.1 mJy beam⁻¹, while an additional 0.05 mJy beam⁻¹ systematic error from amplitude self-calibration was introduced to the 1.3 mm high-resolution continuum data, yielding a final continuum image with a rms noise of 0.15 mJy beam⁻¹ and a synthesized beam of 0.031′′ × 0.025′′ (position angle: –45°).

2.2.2 Low angular resolution

Band 6 (1.3 mm) low angular resolution observations were carried out with ALMA in 2013 July 18 as part of Cycle 2 program 2013.1.00994.S. The project was carried out using 38 antennas with a diameter of 12 m. The baselines range from 124 to 1573 m (9.164–1573.35 kλ), corresponding to the largest recoverable angular scale of 9.9′′. The phase center of the observations is at the sky position RA(J2000) = 19^h23^m39^s.95, Dec (J2000) = +14°31′05′′.50. The continuum image was obtained by averaging all channels from four spectral windows, which have a bandwidth of 1.875 GHz and they are centered at rest frequencies of 224.984 (spw0), 226.984 (spw1), 239.015 (spw2), and 241.015 GHz (spw3). The quasars J1924−2194 and J1751+0939 were used as bandpass calibrators. The flux and phase calibrators are Titan and J1922+1530, respectively. The aim of this project is to search for polarization in W51 IRS2, consequently the polarization calibrator was J1924-2914, but we only took data from the Stokes ‘I’ parameter to generate the continuum image.

We performed seven rounds of phase self-calibration and one round of amplitude self-calibration to the visibilities. The phase-only self-calibrations decreased the rms noise of the 1.3 mm low-resolution continuum image from 1.5 to 0.7 mJy beam⁻¹. The amplitude self-calibration did not introduce any systematic errors, yielding a final continuum image with a rms noise of 0.5 mJy beam⁻¹ and a synthesized beam of 0.237′′ × 0.217′′ (position angle: 28°).

Fig. 1 Upper panel: color image and white contours represent the 1.3 mm high-resolution continuum map of the W51 IRS2 region. Contour levels are [−5, −3, 3, 5, 10, 20, 40, 80, 120, 160] × σ (σ = 0.15 mJy beam−1). Red contours correspond to the 3 mm continuum map, and they are [−3, 3, 5, 10, 20, 40, 80, 120, 160, 200, 300, 400, 500] × σ (σ = 0.08 mJy −1). The red and white ellipses in the bottom-right corner indicate the synthesized beams of the 3 and 1.3 mm high-resolution continuum maps, respectively. Bottom panel: color image represents the spectral index map that was derived from the 1.3 and 3 mm high-resolution continuum maps. The black contours are the 1.3 mm high-resolution continuum, and the contour levels are the same as in the upper panel.

3 Results

3.1 Dust emission and spectral index

In the upper panel of Fig. 1, we present the 3 and 1.3 mm ALMA high-resolution continuum images obtained in this study. The color scale and white contours correspond to the intensity distribution of the 1.3 mm high-resolution continuum, and the red contours correspond to the intensity distribution of the 3 mm continuum, respectively. The sources corresponding to W51 North and W51 d2 are labeled in this figure. W51 IRS2 was previously observed by several projects, where four clumps (SMA1–4) were reported in this region following a filamentary structure along the east-west direction (Zapata et al. 2009; Tang et al. 2013; Koch et al. 2018). The upper panel of Fig. 1 shows that most of the continuum objects are well resolved, and the majority of these sources have 3 mm counterparts. Most of the detected objects are part of elongated structures, consistent with previous observations. Our data confirm that W51 IRS2 is a complex and structured region.

The bottom panel of Fig. 1 presents a spectral index map that was derived from the 3 and 1.3 mm ALMA high angular resolution continuum images. The 1.3 mm high-resolution continuum data are overlaid in the spectral index map as black contours to indicate the structure of the continuum emission. The spectral index map was derived from the “spix” mode in the “immath” task of the CASA package (McMullin et al. 2007), which computes the spectral index using two continuum images with different frequencies based on the relationship between flux density and spectral index of S_ν ∝ν^α. The spectral index is defined as α = ln(I_3mm∕I_1.3mm)∕ln(ν_3mm∕ν_1.3mm), where I is the intensity value of each pixel and ν is the frequency. Since these two continuum images have different synthesized beams, before generating the spectral index map, we tapered them to the same resolution of 0.079′′ × 0.053′′ using a pixelsize of 0.01′′. In order to produce this spectral index map, we also included a cut-off with rms noises higher than 3σ. This helped us discriminate the fake sources in the field. Based on the spectral index map, most sources have spectral indices of around 2, suggesting that the dust emission may be partly optically thick since optically thin dust emission is expected to have a substantially steeper spectral index in the range of 3–4. The averaged spectral index over the W51 North region is ~2.5, which is comparable to the previously reported spectral index of 2.8 in W51 North using Submillimeter Array observations (Zapata et al. 2009). W51 d2 and the region that we call CC1, which are two small regions located to the south of the surveyed area, show spectral indices around 0.2, which is related to optically thin free-free emission from a H_II region. This is also consistent with previous works finding that W51 d2 is a very bright centimeter continuum source (Zapata et al. 2009; Ginsburg et al. 2016, 2017).

3.2 Fragment identification

To study the fragmentation properties, it is necessary to identify fragments in this region. The 1.3 mm high-resolution image is the best option to do this job as it has the highest angular resolution of 0.028′′ (corresponding to a spatial scale of 144 AU at 5128 pc distance). As we mention in Sect. 2.2.1, the largest recoverable angular scale for this dataset is 1.6′′. This means that the flux from the 1.3 mm high-resolution continuum may not only come from compact structures, but it also may come from their extended surrounding environment. As shown in the upper panel of Fig. 1, the majority of the 1.3 mm continuum objects from the continuum image seem to have angular sizes smaller than 0.2′′ (~1000 AU). Therefore, to identify the real compact and small-scale fragments, we constrained the 1.3 mm high-resolution continuum image by using only visibilities with a uv range larger than 910 kλ (corresponding to the largest recoverable angular scale of 0.1′′). However, cutting off short baselines may lead to missing flux problems and create artifacts in the continuum image. In order to test whetherthe short baselines removed here can significantly affect our results, we simulated the ALMA observations using the “SIMOBSERVE” task in the CASA package. The detailed description and results of the simulations are presented in Appendix A. Based on the results of our simulations, it seems that cutting off baselines shorter than 910 kλ does not create significant artifacts, but it allowed us to recover the real compact sources.

Figure 2 shows the 1.3 mm continuum image obtained from the data with a uv range larger than 910 kλ. This restriction on the uv range results in a rms noise of σ = 0.11 mJy beam⁻¹. In this figure, the white contours correspond to the uv-range constrained 1.3 mm continuum. To indicate the original 1.3 mm high-resolution continuum distribution, we overlaid the 6σ contour of the full uv-range 1.3 mm high-resolution continuum image in Fig. 2 (red contour). To give close-inviews of the continuum objects, we divided the whole region into seven subclumps (clump A to G are named CA, CB, CC, CD, CE, CF, and CG), and the subplots are also presented in Fig. 2.

In this work, we assumed an identification threshold of 12σ¹ in the uv-range constrained continuum image. If the 12σ contour corresponds to an isolated source and presents no substructure, this was counted as one single source. On the contrary, if the 12σ contour contains further substructure, we considered different sources within this 12σ contour if they have at least one different closed contour. In Fig. B.1 we illustrate this identification criterion by indicating the 12σ contour for isolated sources along with each first closed contour (above the threshold) that unambiguously separates each identified source². Based on this criteria, a total of 33 continuum sources were identified within the surveyed region, and the positions of the identified continuum sources are marked as red crosses in Fig. 2. A two-dimensional Gaussian fitting was performed toward each identified source, and the fitting results are listed in Table 1.

Generally, the identified sources are embedded within a more extended environment. However, there are some isolated and point-like sources, which seem to be already dynamically decoupled from their parental filamentary structures, and they are probably objects with a different nature. Therefore, in order to avoid over-identification, we made a classification for all identified continuum sources based on the following two criteria. First, the red contour in Fig. 2 is 6σ of the 1.3 mm continuum with full uv range; if a source is enclosed by a red contour individually, and the continuum source does not have an extendedstructure, it is considered as an isolated source, which means that the source was already decoupled from its parental structure. Second, if a source has a ratio of intensity to flux density higher than 0.5, it is considered as an unresolved point-like source. The isolatedand point-like sources were not identified as fragments and were excluded from further analysis and marked with “†” in Table 1. Finally, 29 out of 33 continuum sources were identified as fragments. It should be noted that the terminology of “fragment” in this work just stands for the smallest compact continuum sources which follow the identification criteria described above.

Assuming that the 1.3 mm continuum emission of the identified sources is optically thin, their masses can be estimated using the flux densities as: $$M_{\text{source}}=\frac{S_{\nu }{D}^2}{B(T_{\text{d}}){\kappa }_{\nu }},$$(1)

where, D = 5128(±1867) pc is the distance to W51 IRS2 (Reid et al. 2019) and B(T_d) is the Planck function at temperature T_d. Ginsburg et al. (2017) derived a gas temperature range of 200 K < T_gas < 390 K for W51 North by using methanol lines. In this case, since our surveyed region is much larger than W51 North, we assume a dust temperature of T_d = T_gas = 200 K and that the dust temperature is uniformly distributed in the surveyed region (this assumed temperature has also been used in Goddi et al. 2020). We note that κ_ν is the dust opacity, and the value of κ_ν from Col. 6 of Table 1 in Ossenkopf & Henning (1994) was adopted to estimate the densities of continuum sources (κ = 0.008991 cm² g⁻¹ at 1.3 mm). Furthermore, S_ν is the total flux density measured in the 1.3 mm uv-range constrained continuum image, and it was derived from the Gaussian fitting. The masses of the continuum sources are listed in Table 1.

As mentioned in Sect. 3.1, the majority of sources detected in the surveyed region have spectral indices of ~2, suggesting that they may be partly optically thick. Ginsburg et al. (2017) also reported that the central region of W51 North seems to be optically thick. Assuming that these sources are slightly optically thick, the column density of these sources can be estimated by using $N_{{\text{H}}_2}$=$\frac{1}{{\kappa }_{\nu }}$ = 2.6 × 10²⁵ cm⁻² at 1.3 mm. The masses of the optically thick sources were estimated by M_t = $\pi R^2{N}_{{\text{H}}_2}{m}_{\text{H}}\mu$, where R is the radius of the sources derived from the Gaussian fitting, μ = 2.8 is the mean molecular weight (Kauffmann et al. 2008), and m_H is the mass of H atom. The masses based on the optically thick assumption are also listed in the last column of Table 1. The typical (median) mass estimated by Eq. (1) is 0.4(±0.3) M_⊙, and the typical mass estimated with the optically thick assumption is 0.5(±0.3) M_⊙. This indicates that the optically thin assumption does not strongly affect our results. Hereafter, our calculations are based on the assumption that the dust emission is optically thin.

Fig. 2 Color image and white contours correspond to the uv-range constrained 1.3 mm high-resolution continuum (with a uv range larger than 910 kλ). The white contour levels are [−12, −6, 6, 12, 18, 24, 30, 36, 40, 60, 80, 100, 120, 140, 160] × σ (σ = 0.11 mJy beam−1), and the beam size is shown in the bottom-right corner of each plot. The red contour in each plot is the full uv-range 1.3 mm high-resolution continuum image, and the contour level is 6σ (σ = 0.15 mJ beam−1). We identified fragments as sources with at least one closed contour above the 12σ (1.32 mJy beam−1) identification threshold. A total of 33 continuum objects have been identified in the region. To illustrate the “closed contour” identification method, we plotted the first closed contour above the threshold that separates each identified source as an orange contour in Fig. B.1. To give a better view of the fragments, we divided the whole region into seven subclumps (clump A to G are named CA, CB, CC, CD, CE, CF, and CG), shown in the individual close-in panels, where contours are shown in black to improve the contrast over the background. The crosses in the subplots indicate the central positions of the identified fragments. The names of the fragments are labeled. A Gaussian fitting was performed to obtain the basic parameters of the identified objects and the fitting results are listed in Table 1.

4 Analysis

4.1 Typical separations between fragments

To study the fragmentation process in W51 IRS2 more quantitatively, we estimated the angular separations of each fragment to every other fragment by using their coordinates listed in Table 1, and we considered a pair of fragments as one companion. Thus, the 29 identified fragments generate 406 companions, and every companion has its own corresponding angular separation. We measured the angular separations (θ) of all companions (fragments pairs) in degrees, then we grouped the angular separations of all companions into different annuli with logarithmic bins of log(θ) = 0.1, and we plotted the histogram of these grouped separations in the upper panel of Fig. 3. In this figure, one can see that there are three peaks located in the plot with steps in log(θ) of –4.3, –3.6, and –2.9, corresponding to 925, 4635, and 23 229 AU, respectively. These three peaks are indicated by gray vertical dashed linesin the upper panel of Fig. 3.

A large-scale clump of about ~0.1 pc with its initial physical conditions seems to fragment into medium-scale clumps of ~0.02 pc, then these medium-scale clumps continue to fragment into small-scale objects of ≲0.004 pc and largedensities. Therefore, the whole large-scale clump presents a hierarchical structure with fragmentation taking place at different scale levels. This is called hierarchical fragmentation, and it has been found in different star-forming regions (e.g. Wang et al. 2011; Takahashi et al. 2013; Zhang et al. 2015; Teixeira et al. 2016; Kainulainen et al. 2017). A peak in the histogram indicates that the number of companions is highest at the corresponding scale, meaning that the clustering phenomenon is present at this scale.

The three peaks of the histogram indicate that there seems to be a multilevel hierarchical structure in our surveyed region. However, the third peak is meaningless for our further analysis because it corresponds to the largest spatial scale of the whole surveyed area. Therefore, we only take the first and the second peaks of the histogram into consideration. To a certain extent, the typical separation corresponding to each peak could also reflect the physical process controlling fragmentation at the corresponding scale. For instance, if fragmentation is consistent with a thermal Jeans process, the typical separation defined by the histogram peak should be comparable with its Jeans length. We compare the typical separations with Jeans lengths in the next section.

The correlation between the mean surface density of companions (MSDC) and angular separations (θ) are generally used to study the source distribution and fragmentation of molecular clouds (Gómez et al. 1993; Larson 1995; Bate et al. 1998; Kraus & Hilenbrand 2008; Tafalla & Hacar 2015; Palau et al. 2018). We calculated the MSDC (Σ(θ)) of the identified 29 fragments by following the description of Simon (1997): $$\Sigma (\theta )=\frac{N_{\text{p}}}{2\pi r_{\theta }\delta (\theta )N_{\text{f}}},$$(2)

where N_p is the number of fragment pairs (or the number of companions) that fall in the annulus with a radius r_θ and width δ(θ), and N_f is the total number of fragments. The logarithmic relationship between Σ(θ) and the angular separation (θ) is shown in the bottom panel of Fig. 3. The orange squares denote the averaged Σ(θ) of each annulus, and the error bars correspond to the standard deviation. In order to compare the MSDC of our identified fragments with the MSDC corresponding to uniformly, randomly distributed sources, we ran 1000 rounds of calculations of MSDC for a uniformly random distribution of sources within the same area. We plotted the averaged result and errors over 1000 rounds of random distributions as a black solid line and gray shaded region in the bottom panel of Fig. 3, respectively.

Only the region from –4.0 ≤ log(θ) ≤ –3.4 shows a flat behavior (within errors), consistent with our observed fragments being randomly distributed (blue area in Fig. 3). We define a break as the spatial scale where the MSDC is no longer consistent with the random distribution, and it is related to the spatial scale of the parental structure undergoing fragmentation. From the bottom panel of Fig. 3, two breaks are found at log(θ) = –4.0 (1845 AU) and log(θ) = –3.4 (7346 AU). In the bottom panel of Fig. 3, the orange vertical dashed lines indicate the breaks of Σ(θ).

The values of the MSDC for separations larger than the second break at log(θ) = –3.4 (7346 AU) are below the gray area. This is because the surveyed area is undergoing hierarchical fragmentation, but our random sources were generated within the entire surveyed region. If there was only a one-level fragmentation process, the flat region would extend up to large scales, and only one break would appear in the plot. On the other hand, the values of the MSDC for separations smaller than the first break at log(θ) = –4.0 (1845 AU) are clearly located well above the gray shaded area, revealing an enhancement of mean surface density, andthus a significant clustering phenomenon. In addition, these points at separations smaller than log(θ) = –4.0 (1845 AU) presenta linear behavior, mainly contributed to by the fragments in W51 North. We performed a linear fit toward these points and found a slope of –0.8(±0.1) (the result is plotted as a red dashed line in Fig. 3). This value is very consistent with the values obtained for the linear regime of the MSDC in the Pipe Nebula (–0.68) and the Lupus I cloud (–0.72) (Román-Zúñiga et al. 2019), where this linear regime is found at similar or slightly larger spatial scales of 2000–8000 AU.

Fig. 3 Upper panel: histogram of angular separations of companions. The Y-axis is the number of companions, and the X-axis is angular separations. The gray vertical dashed lines represent the peak number of the companions. Bottom panel: orange squares with error bars represent the mean surface density of companions (MSDC) versus the angular separation of the identified fragments. The orange vertical dashed lines denote the breaks of the MSDC. The black solid line and the gray region represent MSDC and errors that were derived from the 1000 rounds of uniformly, and randomly distributed 29 fragments in the same region. The linear fitting toward the data points at separations smaller than the first break, marked as a red dashed line in the bottom panel, results in a slope of –0.8(±0.1).

4.2 Ensembles defined as groups of fragments

In Sect. 4.1, two MSDC breaks were found corresponding to 1845 and 7346 AU. As we mentioned before, these breaks reflect the beginnings of clustering phenomena, and they are related to the size of parental structures that undergo fragmentation. Thus, we use such spatial scales to define groups of fragments at different spatial scales. In this way, we draw a circle with aradius corresponding to the identified breaks for every fragment. If a circle includes the center of any other fragment, we consider them to be related fragments that belong to the same “parental structure”. For the convenience of the subsequent discussion, we define an “ensemble” as a group of subobjects belonging to the same parental structure (following Zhang et al. 2021 nomenclature).

Based on the above criteria, a total of five ensembles were identified with a spatial scale of 1845 AU. The center of each ensemble was calculated by taking the geometric average of the positions of the fragments that belong to the same ensemble, and a spatial scale of 1845 AU was adopted as its radius. Taking a spatial scale of 7346 AU as a radius, only two ensembles were identified. Since the spatial scale of the first break is smaller than the second one, we call the ensembles identified using the 1845 AU radius “small-scale ensembles” (SEs). Each small ensemble is referred to as an SE plus numbers. The ensembles identified using the 7346 AU radius are called “large-scale ensembles” (LEs). The coordinate of each LE is the geometric average of the positions of the SEs that belong to the same LE. The details regarding the ensembles are listed in Table 2. The identified ensembles are shown as small and large circles in Fig. 4, and the names of SE and LE are labeled in the figure. It should be noted that every fragment is assigned to a SE (except a fragment in “CA4” outside the circle of SE1, which is still assigned to SE1).

From Fig. 4, the hierarchical fragmentation scenario in this region is well illustrated. The entire region fragments into two separated subregions, LE1 and LE2, then each large-scale ensemble fragments into different small-scale ensembles, respectively. Finally, each small-scale ensemble splits up into several fragments that have been identified using the uv-range constrained continuum image.

4.3 Fragmentation

In order to investigate the fragmentation process occurring in our surveyed region, we estimated the fragmentation levels and densities of each ensemble. We defined the fragmentation level of SEs as the number of identified fragments belonging to that SE. For LEs, we defined their fragmentation level as the number of SEs belonging to the corresponding LE. The fragmentation level of each ensemble is listed in Table 2.

Assuming that the identified ensembles are spheres with a uniform density distribution, then their densities can be estimated as follows: $$n_{\text{ens}}=\frac{M_{\text{ens}}}{\frac{4}{3}\pi R^3\mu m_{\text{H}}},$$(3)

where R is the radius of the ensemble, which has been considered to be 1845 AU for SEs and 7346 AU for LEs, μ = 2.8 is the mean molecular weight (Kauffmann et al. 2008), m_H is the mass of a H atom, and M_ens is the mass of an ensemble, which can be estimated from the total flux density within the circle region in Fig. 4 (radii of 1845 AU for SEs and 7346 AU for LEs) following Eq. (1). A dust temperature of T = 200 K was used to estimate the masses of the ensembles, assuming that the dust temperature is constantly distributed in the region. Since the 1.3 mm uv-constrained high-resolution data are employed to identify fragments and as the 3 mm data are affected by the free-free emission, the data set at 1.3 mm with low angular resolution (0.227′′, Sect. 2.2.2) was used to estimate the masses of the ensembles. This data set has the largest recoverable angular scale of ~ 10′′, which covers the surveyed region very well. As we mentioned before, the high-resolution 1.3 mm continuum image with a uv range larger than 910 kλ was adopted to identify fragments. In order to estimate the densities of the ensembles from the 1.3 mm low-resolution data, we generated a 1.3 mm low-resolution continuum image using only the uv range of 9.164 ~ 910 kλ. Since the 1.3 mm low-resolution data has a resolution of 0.227′′ (1163 AU), which is smaller than the diameter of SEs, we directly used the flux density from 1.3 mm low-resolution data to estimate the masses and densities of ensembles. The estimated densities of SEs and LEs are listed in Table 2.

As can be seen from Table 2, SE1 has the highest fragmentation level and is also the densest ensemble. Other SEs have densities smaller than SE1, and also lower fragmentation levels. While this is the first result consistent with thermal Jeans fragmentation, further tests need to be performed to make this result more robust.

A second test for thermal Jeans fragmentation is the comparison between the typical fragment separations and masses with the expected values from a Jeans analysis. To do this, the densities derived from the 1.3 mm low-resolution data with a uv range of 9–910 kλ were used to estimate the Jeans lengths and Jeans masses for all the identified ensembles. Assuming an idealized sphere with uniform distributions of density and dust temperature, the Jeans length and Jeans mass can be estimated as follows: $${\lambda }_J=c_{\text{s}}{\left(\frac{\pi }{G\rho }\right)}^{\frac{1}{2}},$$(4) $$M_J=\frac{4}{3}\pi \rho {\left(\frac{{\lambda }_J}{2}\right)}^3,$$(5)

where c_s is the sound speed at temperature T = 200 K, G is the gravitational constant, and ρ is the mass density. The estimated Jeans lengths and Jeans masses for the ensembles are listed in Table 2. We also estimated the typical (median) separation between every two subobjects within each SE and LE, along with the typical (median) masses of the subobjects, and listed the results in Table 2.

From Table 2, SE1 results in a λ_J = 817(±112) AU and M_J = 0.9(±0.1) M_⊙. The Jeans length obtained for SE1 is consistent (within errors) with the typical separations of fragments within SE1. Additionally, as we mentioned before, the first peak shown in the upper panel of Fig. 3 is mainly contributed to by SE1 (W51 North), and the scale identified by this first peak (925 AU) is also consistent with its Jeans length. Similarly, a comparison of the Jeans lengths of all the other identified SEs and LEs to the typical observed separations indicate that they are fully consistent within errors.

Regarding the typical masses of the fragments in SEs, these are generally similar to or even smaller than the Jeans mass. In no case is the difference (after taking the uncertainties reported in the table into account) larger than a factor of 3, and this could be explained by a combination of both dust temperature and dust opacity uncertainties. Thus, for SEs it is plausible to think that the typical measured masses are also consistent with the Jeans mass. As for LEs, the masses of the subobjects (SEs) are in general larger by up to a factor of 5 (LE1) or 2 (LE2) than the Jeans mass. Although the typical mass of subobjects for LE1 is clearly higher than the corresponding Jeans mass, by only considering a temperature of 300 K instead of 200 K in the Jeans analysis, the difference would only be a factor of 3 (instead of 5). These results are further discussed in Sect. 5.1.

Fig. 4 Color image: 1.3 mm low-resolution continuum image with the uv range 9–910 kλ, white contours indicate the 1.3 mm continuum with a uv range larger than 910 kλ. The white contour levels are [−12, −6, 6, 12, 12, 18, 24, 30, 36, 40, 60, 80, 100, 120,140,160] × σ (σ = 0.11 mJy beam−1), and the beam size is shown in the bottom left corner. The red contour in the plot isthe 1.3 mm continuum with a full uv range, and the contour level is 6σ (σ = 0.15 mJ beam−1). The small and larger circles denote the “small-scale ensembles” (SEs) and “large-scale ensembles” (LEs), respectively. These SEs and LEs are defined by the spatial scales corresponding to the first and second breaks from the MSDC analysis (see Sect. 4.2 for details). Overall, five SEs and two LEs are defined in the region, and their names are labeled in the figure. The information about these ensembles is listed in Table 2.

5 Discussion

5.1 Thermal Jeans fragmentation could be at work in W51 IRS2

In Sect. 4.3 two tests for the thermal Jeans fragmentation scenario were performed. First, it was studied if the fragmentation level of the ensembles was related to the density of these ensembles. From this study, it was found that SE1, the ensemble associated with W51 North, presented both the highest fragmentation level along with the highest density. In addition, the Jeans length and mass estimated for SE1 are very comparable to the measured median separations and masses of the fragments in SE1. Thus, from our previous analysis it seems that the observed fragmentation in SE1 is consistent with thermal Jeans fragmentation.

While this scenario seems quite clear for SE1, the low fragmentation levels in the other SEs is more intriguing. In fact, one can estimate the number of fragments expected for each SE by calculating the Jeans number: $$N_{\text{Jeans}}=\frac{M_{\text{ens}}\text{CFE}}{M_{\text{Jeans}}},$$(6)

where CFE stands for core formation efficiency and is a factor accounting for the fact that, given a density structure, it is not expected that it is entirely fragmented at the same time, but that the highest density substructures fragment first and thus there must always be a fraction of the parental structure which has not fragmented yet (e.g., Palau et al. 2015; Vázquez-Semadeni et al. 2019). If we calculate the Jeans number for SE1 assuming that CFE = 100%, then N_Jeans = 83, while the number of observed fragments is 20, indicating that assuming CFE = 100% is probably not correct. In order to have N_Jeans = 20 in SE1, the CFE must be ~24%, which is consistent with the values found by Palau et al. (2015) for a sample of massive dense cores. Thus, one can calculate N_Jeans for the other SEs, adopting the same CFE obtained for SE1. The resulting N_Jeans are 4, 5, 4, and 2 for SE2, SE3, SE4, and SE5, respectively. This expected number of fragments for the SEs matches the observed number of fragments well, especially for the case of SE2 and SE5, while the case deviating more strongly from this prediction is SE3 where only one fragment is identified while five fragments were expected.

This suggests that either other mechanisms in addition to thermal pressure are at work in SE3, or that the Jeans mass has been underestimated for SE3. In the first case, one could think of both turbulence or magnetic fields as possible additional supporting mechanisms against gravity. However, the role of turbulence is greatly debated in the community because in order to act as an effective form of support against gravity, it should be isotropic, which is highly questioned because turbulence is characterized by an energy spectrum containing most of the energy at large scales (see Sect. 2.2 of Vázquez-Semadeni et al. 2019, see also Dobbs et al. 2005; Burkhart & Mocz 2019). In addition, this has not been supported by observational tests (e.g., Palau et al. 2015). Regarding the support offered by the magnetic field, there are a number of theoretical and numerical works suggesting that the magnetic field should be efficient at suppressing fragmentation (e.g., Boss 2004; Vázquez-Semadeni et al. 2005; Commerçon et al. 2011a; Hennebelle et al. 2011), and also very recent observational works seem to be consistent with this prediction (e.g., Añez-López et al. 2020; Palau et al. 2021). However, for the case of W51 IRS2, there are magnetic field measurements by Koch et al. (2018) and the magnetic field does not seem to be particularly strong in SE3 compared to SE1.

In conclusion, it seems that the most plausible explanation for the low fragmentation found in SE3 is that it might be warmer than the other SEs naturally inhibiting fragmentation. If a temperature of 400 K is assumed in the Jeans calculation, the Jeans mass increases to ~ 4 M_⊙, and the Jeans number decreases to ~2, consistent with our observations within the Poisson uncertainty. In order to test if SE3 is actually hotter than the other ensembles of W51 IRS2, we checked the 1.3 mm molecular line data of W51 IRS2 (high-resolution dataset), and vibrational transitions of CH₃CN (v8 = 1) have been found in absorption in W51 North, W51 d2, and SE3 (Tang et al., in prep.). This indicates that these three regions should have gas temperatures higher than the other identified ensembles, and this is consistent with the assumption that SE3 has a temperature higher than 200 K. We note that the high temperatures required in the Jeans analysis to explain our observations indicate that W51 IRS2 is experiencing a reduced or suppressed fragmentation event, that is if the ambient or parental cloud had not been heated by nearby massive stars, W51 IRS2 would most likely undergo a more powerful fragmentation event in terms of the number of fragments.

In summary, if temperatures for the Jeans analysis are adopted to be in the range 200–400 K, then our observations are very consistent with the expected numbers for thermal Jeans fragmentation. Our observations suggest that the fragmentation in W51 IRS2 is thermally inhibited, as already suggested by Ginsburg et al. (2017), and that a thermal Jeans fragmentation process is probably at work with an unusually high temperature due to the stellar feedback from the nearby massive stars.

5.2 A cluster-like spatial distribution of fragments in W51 North (SE1)

In Sect. 4.1, a linear regime in the MSDC plot of W51 IRS2 was reported for separations < 1845 AU (Fig. 3-bottom), and the fragments corresponding to these separations are essentially the fragments associated with W51 North. The value of the slope of the linear regime is related to the physical process determining the separation between sources. As already explained in Larson (1995) and Simon (1997), a particular scale is identified in several star-forming clouds above which a slope of − 0.6 was measured in the MSDC plot, associated with the “clustering regime”. On the other hand, below this particular scale, a slope of − 2 was measured for the same clouds, associated with the so-called “binary regime”. The clustering regime is associated with the fragmentation of the cloud, while the binary regime is affected by additional processes such as the dynamical interactions between thedifferent members of a multiple system. The value of the slope found here for W51 North (SE1), –0.8 ± 0.1 (Sect. 4.1), is fully consistent with the expected value for the clustering regime in other clouds of the Milky Way such as Taurus, Ophiucus, the Pipe Nebula, Lupus I, or OMC-1S (Gómez et al. 1993; Simon 1997; Palau et al. 2018; Román-Zúñiga et al. 2019)and clearly shallower than the typical value of the binary regime. Although W51 North is very young and part of the fragments identified here could merge with time, the fact that the slope of the MSDC plot for W51 North precisely corresponds to the expected value for the clustering regime suggests that a protocluster is in the making in W51 North. However, one should keep in mind that if a fraction of the fragments merge with time, the slope of the linear regime could change. Further kinematic studies should be carried out to explore the possible evolution of the fragments.

It is important to compare the slope and scale of the linear regime found in Fig. 3 for W51 North to other star-formingregions also presenting the clustering regime. The scale at which the MSDC changes from a clustering to binary regime (what we call here the “clustering-to-binary scale”) is found to vary from region to region and, interestingly, Simon (1997) reported an anticorrelation between the clustering-to-binary scale and the stellar density of each cloud. In order to test if both W51 IRS2 and OMC-1S fit the anticorrelation found by Simon (1997), we made a rough estimate for the stellar density in W51 IRS2 and OMC-1S by considering that the fragments identified here are true (precursors of) protostellar objects. Although most of the fragments identified in this work have no reported signs of stellar activity such as outflows, such an estimate for the stellar density is probably accurate as a first order of magnitude because of several reasons. First, the identified fragments have masses that are generally consistent with the Jeans mass (within uncertainties) and therefore could be gravitationally unstable. Second, thanks to the uv restriction used to identify fragments (Sect. 3.2), only the fragments that are already very compact and centrally peaked have been considered, and we have left out the larger noncentrally peaked structures which could be transient structures. Using the data reported in Table 1, we calculated the densities of the fragments and their median value is 8 × 10⁹ cm⁻³. Such a high density is precisely the density at which the first hydrostatic core, the direct precursor of the protostellar object, formed in the model of Larson (1969), and in subsequent more recent works (e.g., Commerçon et al. 2011b; Bhandare et al. 2018). Third, although there are no infrared or X-ray counterparts associated with the identified fragments due to a lack of observations with a similar resolution as our 1.3 mm continuum data, Palau et al. (2018) identified fragments in OMC-1S adopting the same method as described in Sect. 3.2, and they found that 50% of them are associated with star-forming signposts. Therefore, based on the above reasons, it is feasible that an important fraction of the fragments identified here are already protostellar or will become protostellar in the near future. Another fraction of the fragments detected here might not collapse either because of thermal support, mergers of fragments, or tidal forces stretching the fragments (e.g., Lee & Hennebelle 2018; Hennebelle et al. 2019).

Thus, the stellar density for both W51 IRS2 and OMC-1S was estimated by dividing the number of identified fragments in each region by the volume corresponding to the area where these fragments have been identified (assuming they are spherically distributed). The resulting stellar densities for both W51 IRS2 and OMC-1S, along with the previously reported stellar densities for Taurus, Ophiucus, and the Trapezium, are listed in Table 3. Given the reasons listed in the previous paragraph, a 50% uncertainty is assumed for the stellar density estimation.

Regarding the clustering-to-binary scale, for both W51 IRS2 and OMC-1S, the MSDC clearly reveals a linear regime with a slope in the range of − 0.75 to −1, which is well below the slope for the binary regime, but the transition to the binary regime itself is not detected. Thus, we consider that the lower limit of our range of separations in the MSDC plot is the upper limit of the clustering-to-binary scale for W51 IRS2 and OMC-1S.

The anticorrelation discovered by Simon (1997) is presented in Fig. 5, where the data for both W51 IRS2 and OMC-1S have already been included. The figure shows that the two regions added here follow the anticorrelation reasonably well between the clustering-to-binary scale and the stellar density, but now at higher stellar densities. This finding, along with the slope value for the linear regime associated with W51 North, which is very similar to the values reported for the clustering regime in other clouds of the Galaxy, indicates that W51 North seems to be forming a protocluster with spatial properties very similar to the properties of other clustered regions in the Galaxy.

In Sect. 5.1, it is shown that fragmentation in W51 IRS2 is consistent with a thermal Jeans process, but with the particular case of high temperatures (200–400 K), probably produced by the stellar feedback from nearby massive stars. Therefore, fragmentation in W51 IRS2 has probably been thermally inhibited (as already suggested by Ginsburg et al. 2017), and much more fragments might have formed in the region if it had not been previously heated to such high temperatures. However, the inhibited fragmentation does not seem to prevent the formation of a protoclusterassociated with W51North as the spatial distribution of the associated fragments (most likely being or becomingprotostellar) are fully consistent with the spatial distribution of typical cluster-like regions in the Galaxy.

Finally, we would like to mention that not only does W51 IRS2 fit the anticorrelation, but OMC-1S also corresponds to stellar densities that are even higher. The fact that the slopes of the clustering regime in the MSDC plot are very similar in all the regions compiled here (between − 0.6 and − 1), which in turn cover a very broad range of stellar densities, is consistent with the fractal dimension (ratio of the change in detail to the change in scale) of the molecular cloud structure being similar for a wide range of gas densities and spatial scales, as already suggested by Simon (1997).

Fig. 5 Anticorrelation between the clustering-to-binary scale in the MSDC plot and stellar density for different clouds. The data point of W51 IRS2 is taken from this work; the data point of OMC-1S is taken from Palau et al. (2018); and the data points of Taurus, Ophiucus, and Trapezium are taken from Gómez et al. (1993) and Simon (1997).

6 Conclusions

To study the fragmentation properties of the W51 IRS2 region, ALMA data from three projects have been used. Based on our analyses, the main findings of this work are summarized as follows:

1. A spectral index map was derived from 3 to 1.3 mm high-resolution continuum images. Based on the spectral index map, we found that most areas in the surveyed region have spectral indices around 2.0, with only two regions, W51 d2 and CC1, associated with spectral indices around 0.2 indicative of free-free emission, indicating that the majority of the detected sources might be partly optically thick.

2. A 1.3 mm uv-range constrained high-resolution continuum image (uv range >910 kλ) was generated to identify fragments, which is sensitive only to structures with angular scales smaller than 0.1′′ (~500 AU). Based on a 12σ identification threshold, a total of 33 spatially separated continuum objects have been identified and 29 of them are defined as fragments, with masses ranging from 0.2 to 2 M_⊙.

3. The mean surface density of companions (MSDC) analysis reveals two spatial scales at 1845 and 7346 AU associated with overdensities, suggesting a two-level hierarchical fragmentation. We used these two spatial scales to define ensembles as groups of fragments. Five small-scale ensembles (SEs, corresponding to a radius of 1845 AU) and two large-scale ensembles (LEs, corresponding to a radius of 7346 AU) are defined in our surveyed region.

4. The fragmentation level and densities of both SEs and LEs have been estimated. The ensemble with the highest density, associated with W51 North, corresponds to the ensemble with the highest fragmentation level. Combining this with a Jeans analysis, it seems that fragmentation in SEs and LEs is consistent with a thermal Jeans process if high temperatures in the range from 200–400 K are considered, naturally explaining the small number of fragments along with their masses and separations.

5. A linear regime is found in the MSDC plot for separations <1845 AU, with a slope of –0.8 ± 0.1. Such a slope is very similar to the slope observed in previous works for the clustering regime in the MSDC, and the anticorrelation found by Simon (1997) between the transition scale from the clustering to the binary regime and the stellar density seems to hold for the higher stellar densities studied here.

Overall, W51 IRS2 seems to be a case of thermal Jeans fragmentation operating at high temperatures due to stellar feedback from nearby massive stars, implying a thermal inhibition of fragmentation in the region. However, such an inhibition does not seem to be powerful enough to prevent the formation of a protocluster in W51 North, for which the spatial distribution of the future protostars seems to be fully consistent with the spatial distribution of well-known clusters in the Galaxy.

Acknowledgements

The authors deeply appreciate the crucial comments from the referee, which have undoubtedly improved the quality of this paper. A.P. is grateful to Enrique Vázquez-Semadeni for very insightful discussions. This paper makes use of the following ALMA data: project IDs: 2013.1.00994.S, 2015.1.01596.S and 2017.1.00293.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC 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. L.A.Z. acknowledges financial support from CONACyT-280775, and UNAM-PAPIIT IN110618 grants, México. A.P. acknowledges financial support from the UNAM-PAPIIT IN111421 grant, the Sistema Nacional de Investigadores of CONACyT, and from the CONACyT project number 86372 of the ‘Ciencia de Frontera 2019’ program, entitled ‘Citlalcóatl: a multiscale study at the new frontier of the formation and early evolution of stars and planetary systems’, México. S.-L.Q. is supported by the key research project of the National Natural Science Foundation of China (project number: 12033005). M.T. is supported by China Scholarship Council (CSC), and PhD research startup foundation of Chuxiong Normal University.

Appendix A ALMA simulations

In Section 3.2, we used the 1.3 mm high-resolution uv-range constrained continuum image to identify the continuum objects. However, cutting off short baselines affects the uv coverage and may create artifacts in the final continuum image. In order to study whether cutting-off baselines shorter than 910 kλ can significantly affect our results, we simulated ALMA observations using the “SIMOBSERVE” task in the CASA package.

However, CASA does not provide the specific antenna list used for the 1.3 mm high-resolution observations. Thus, we manually generated an antenna list completely consistent with the 1.3 mm high-resolution observations to perform these simulations. Fig A.1 presents the results of two different simulations with full uv coverage and constrained uv coverage. Panel (a) shows the 1.3 mm high-resolution continuum image with a full uv range. Firstly, we performed a Gaussian fitting to W51 North, then we took the ellipse obtained by Gaussian fitting as the source size to simulate an elongated structure with a comparable spatial scale to W51 North. Secondly, we ran two different simulations. In the first run, the model did not include any embedded compact objects. In the second run, two point sources were included, and the typical intensity of the identified continuum sources in Table 1 were considered to be the peak intensity of the point sources. Finally, the TCLEAN task in CASA was employed to generate continuum images for the two simulations. The realistic residual map from TCLEAN of the 1.3 mm continuum observations was added to the simulated continuum images as the noises. Panels (b) and (d) show continuum images with a full uv range generated from the two different simulations, respectively. According to panels (b) and (d), there is almost no difference between the two images. The continuum images of two different simulations with only a uv range > 910 kλ are presented in panels (c) and (e), respectively. One can see that two compact sources appear in panel (e), but no compact source is found in panel (c), indicating that cutting off baselines shorter than 910 kλ does not create strong artifacts in the 1.3 mm high-resolution image, but this allows us to identify the embedded sources.

Short baselines are sensitive to large-scale emission, and cutting off short baselines may lead to severe missing flux problems. To test how our continuum images change with the uv-range cutoff, we performed several simulations toward different elongated structures with a variety of uv-range cutoffs, as shown in Fig A.2. In panels (a1) - (a5) of Fig. A.2, the source size has been set as 0.1′′ ×0.05′′ (ratio of the major axis to minor axis of 2:1), which is similar to the typical size (0.06′′×0.03′′) of continuum objects identified in Section 3.2. This simulation can be used to investigate how the morphology of a typical compact source changes with a uv-range cutoff. In panels (b1) - (b5), the source size has been set to 0.25′′×0.05′′ (ratio of the major axis to minor axis of 5:1), which reveals how an elongated structure changes with a uv-range cutoff. In panels (c1) - (c5), the source size has been set to 0.5′′ ×0.05′′ (ratio of themajor axis to minor axis of 10:1), which presents how an extremely elongated structure changes with uv-range cutoffs. The first to fifth columns of Fig. A.2 present the continuum images with a full uv range, 910 kλ, 1500 kλ, 3000 kλ, and 5000 kλ uv-range cutoff, respectively. From panels (a1) and (a2) of Fig. A.2, one can see that theintensity and morphology are almost the same in both a full uvrange and 910 kλ cutoff, indicating that a 910 kλ uv-range cutoff does not induce severe missing flux problems for the typical source sizes of our observations. When comparing panels (b1), (b2), (c1), and (c2), we found that none of the elongated structures are split in the 910 kλ uv-range cutoff, which indicates that cutting off baselines shorter than 910 kλ does not create significant artifacts. Raising up the uv-range cutoff threshold induces interference fringes, and some artifacts indeed appear when cutting off a uv range shorter than5000kλ as shown in panel (c5). However, such a uv-range cutoff leads to severe missing flux problems and the peak intensities of artifacts are just ~5×10⁻⁴ Jy beam⁻¹, which is certainly below our 12σ (1.32×10⁻³ Jy beam⁻¹) identification threshold. Therefore, we are confident that cutting off the baselines shorter than 910 kλ does not create artifacts in our 1.3 mm high-resolution image, and the continuum objects identified in Section 3.2 are most likely real compact sources.

Fig. A.1 Panel (a): 1.3 mm continuum image of the W51 North region with a full uv range. The red ellipse is given by Gaussian fitting toward W51 North. This ellipse has been adopted as a model to perform ALMA simulations. Panel (b): Simulated full uv-range continuum image without any embedded source. Panel (c): Simulated uv-range constrained (> 910 kλ) continuum image without any embedded sources. Panel (d): Simulated full uv-range continuum image with two point sources. Panel (e): Simulated uv-range constrained (> 910 kλ) continuum image with two point sources. The contour levels in all the panels start from 12σ (1.32 mJy beam−1). The residual map of the 1.3 mm continuum image generated by the TCLEAN task of CASA has been added to the simulated images to realistically simulate the noise.

Fig. A.2 Simulated continuum images with different source sizes and uv ranges. No simulated sources have embedded point sources. The source sizes and uv-range cutoff thresholds are labeled in each panel. The beam size is shown in the bottom left corner of each panel.

Appendix B Fragments identified using the dendrogram algorithm

The continuum sources listed in Table 1 were identified through an intensity threshold and a "closed contour" criterion applied to a uv-restricted image. In each plot of Fig. B.1, the first closed contour above the threshold that separates each identified source is also presented as an orange contour. To explore whether this identification method can affect the results of our work, we adopted the python-based continuum objects identification algorithm of “dendrogram³” to identify fragments in the surveyed region. The identification threshold was set to 12σ (1.32 mJy beam⁻¹), which is the same as the threshold used in Section 3.2. The differencing threshold (min_delta parameter) of the dendrogram was set to 0.1σ (0.011 mJy beam⁻¹). We considered the “leaves” given by the dendrogram algorithm as the continuum sources. A total of 32 continuum sources were identified bythe dendrogram algorithm, and we plotted them in Fig. B.1. The sources identified in Section 3.2 are also shown in Fig. B.1. We found that the majority of continuum sources identified by the dendrogram algorithm are indeed consistent with the sources identified in Section 3.2, and only one source (CB10) in subregionCB has no counterpart from the dendrogram algorithm. This indicates that our identification method used in Section 3.2 yields results comparable to other identification methods. There is a little offset between the positions of the continuum objects defined by the Gaussian fitting and the dendrogram method. This offset may be caused by the fact that the dendrogram algorithm adopts the continuum peak as the central position of the identified sources, while in Section 3, we used a Gaussian fitting to obtain the central positions of the continuum objects. Taking the beam size into consideration, these offsets do not affect our results or discussion in Sections 3, 4, & 5.

Fig. B.1 Panels (a) to (g) correspond to subregions CA to CG as in Fig 2. The continuum objects identified in this work (Section 3.2) are marked with red crosses, and the continuum objects identified by the dendrogram algorithm are marked with black stars in each panel. The color scale is the 1.3 mm uv-range constrained image, and the white contour corresponds to the 12σ (1.32 mJy beam−1) identification threshold. The beam size is shown in the bottom right corner of each panel. A total of 32 continuum objects have been identified by the dendrogram algorithm in the surveyed region. The first closed contour above the threshold that separates each identified source is marked as an orange contour to illustrate the “closed contour” identification method described in Section 3.2.

References

All Tables

All Figures

Fig. 1 Upper panel: color image and white contours represent the 1.3 mm high-resolution continuum map of the W51 IRS2 region. Contour levels are [−5, −3, 3, 5, 10, 20, 40, 80, 120, 160] × σ (σ = 0.15 mJy beam−1). Red contours correspond to the 3 mm continuum map, and they are [−3, 3, 5, 10, 20, 40, 80, 120, 160, 200, 300, 400, 500] × σ (σ = 0.08 mJy −1). The red and white ellipses in the bottom-right corner indicate the synthesized beams of the 3 and 1.3 mm high-resolution continuum maps, respectively. Bottom panel: color image represents the spectral index map that was derived from the 1.3 and 3 mm high-resolution continuum maps. The black contours are the 1.3 mm high-resolution continuum, and the contour levels are the same as in the upper panel.

In the text

Fig. 2 Color image and white contours correspond to the uv-range constrained 1.3 mm high-resolution continuum (with a uv range larger than 910 kλ). The white contour levels are [−12, −6, 6, 12, 18, 24, 30, 36, 40, 60, 80, 100, 120, 140, 160] × σ (σ = 0.11 mJy beam−1), and the beam size is shown in the bottom-right corner of each plot. The red contour in each plot is the full uv-range 1.3 mm high-resolution continuum image, and the contour level is 6σ (σ = 0.15 mJ beam−1). We identified fragments as sources with at least one closed contour above the 12σ (1.32 mJy beam−1) identification threshold. A total of 33 continuum objects have been identified in the region. To illustrate the “closed contour” identification method, we plotted the first closed contour above the threshold that separates each identified source as an orange contour in Fig. B.1. To give a better view of the fragments, we divided the whole region into seven subclumps (clump A to G are named CA, CB, CC, CD, CE, CF, and CG), shown in the individual close-in panels, where contours are shown in black to improve the contrast over the background. The crosses in the subplots indicate the central positions of the identified fragments. The names of the fragments are labeled. A Gaussian fitting was performed to obtain the basic parameters of the identified objects and the fitting results are listed in Table 1.

In the text

Fig. 3 Upper panel: histogram of angular separations of companions. The Y-axis is the number of companions, and the X-axis is angular separations. The gray vertical dashed lines represent the peak number of the companions. Bottom panel: orange squares with error bars represent the mean surface density of companions (MSDC) versus the angular separation of the identified fragments. The orange vertical dashed lines denote the breaks of the MSDC. The black solid line and the gray region represent MSDC and errors that were derived from the 1000 rounds of uniformly, and randomly distributed 29 fragments in the same region. The linear fitting toward the data points at separations smaller than the first break, marked as a red dashed line in the bottom panel, results in a slope of –0.8(±0.1).

In the text

Fig. 4 Color image: 1.3 mm low-resolution continuum image with the uv range 9–910 kλ, white contours indicate the 1.3 mm continuum with a uv range larger than 910 kλ. The white contour levels are [−12, −6, 6, 12, 12, 18, 24, 30, 36, 40, 60, 80, 100, 120,140,160] × σ (σ = 0.11 mJy beam−1), and the beam size is shown in the bottom left corner. The red contour in the plot isthe 1.3 mm continuum with a full uv range, and the contour level is 6σ (σ = 0.15 mJ beam−1). The small and larger circles denote the “small-scale ensembles” (SEs) and “large-scale ensembles” (LEs), respectively. These SEs and LEs are defined by the spatial scales corresponding to the first and second breaks from the MSDC analysis (see Sect. 4.2 for details). Overall, five SEs and two LEs are defined in the region, and their names are labeled in the figure. The information about these ensembles is listed in Table 2.

In the text

	Fig. 5 Anticorrelation between the clustering-to-binary scale in the MSDC plot and stellar density for different clouds. The data point of W51 IRS2 is taken from this work; the data point of OMC-1S is taken from Palau et al. (2018); and the data points of Taurus, Ophiucus, and Trapezium are taken from Gómez et al. (1993) and Simon (1997).
In the text

Fig. A.1 Panel (a): 1.3 mm continuum image of the W51 North region with a full uv range. The red ellipse is given by Gaussian fitting toward W51 North. This ellipse has been adopted as a model to perform ALMA simulations. Panel (b): Simulated full uv-range continuum image without any embedded source. Panel (c): Simulated uv-range constrained (> 910 kλ) continuum image without any embedded sources. Panel (d): Simulated full uv-range continuum image with two point sources. Panel (e): Simulated uv-range constrained (> 910 kλ) continuum image with two point sources. The contour levels in all the panels start from 12σ (1.32 mJy beam−1). The residual map of the 1.3 mm continuum image generated by the TCLEAN task of CASA has been added to the simulated images to realistically simulate the noise.

In the text

	Fig. A.2 Simulated continuum images with different source sizes and uv ranges. No simulated sources have embedded point sources. The source sizes and uv-range cutoff thresholds are labeled in each panel. The beam size is shown in the bottom left corner of each panel.
In the text

Fig. B.1 Panels (a) to (g) correspond to subregions CA to CG as in Fig 2. The continuum objects identified in this work (Section 3.2) are marked with red crosses, and the continuum objects identified by the dendrogram algorithm are marked with black stars in each panel. The color scale is the 1.3 mm uv-range constrained image, and the white contour corresponds to the 12σ (1.32 mJy beam−1) identification threshold. The beam size is shown in the bottom right corner of each panel. A total of 32 continuum objects have been identified by the dendrogram algorithm in the surveyed region. The first closed contour above the threshold that separates each identified source is marked as an orange contour to illustrate the “closed contour” identification method described in Section 3.2.

In the text