Mass Assembly in Massive Star Formation: A Fragmentation Study of ATLASGAL Clumps

The mass assembly in star-forming regions arises from the hierarchical structure in molecular clouds in tandem with fragmentation at different scales. In this paper, we present a study of the fragmentation of massive clumps covering a range of evolutionary states, selected from the ATLASGAL survey, using the compact configuration of the Submillimeter Array. The observations reveal a wide diversity in the fragmentation properties with about 60% of the sources showing limited to no fragmentation at the 2″ scale, or a physical scale of 0.015–0.09 pc. We also find several examples where the cores detected with the Submillimeter Array are significantly offset from the clump potential, suggesting that initial fragmentation does not result in the formation of a large number of Jeans mass fragments. The fraction of the clump mass that is in compact structures is seen to increase with source evolution. We also see a significant correlation between the maximum mass of a fragment and the bolometric luminosity of the parent clump. These suggest that massive star formation proceeds through clump fed core accretion, with the initial fragmentation being dependent on the density structure of the clumps and/or magnetic fields.


INTRODUCTION
Although massive stars play a dominant role in the physical and chemical evolution of the Universe, several aspects of their formation, especially in the very early stages, are still poorly understood.Recent high resolution observations have shown the presence of gravitational unstable disks around massive stars suggesting that massive stars grow in mass by disk accretion similar to the formation of low-mass stars (e.g.Motogi et al. 2019;Beuther et al. 2017).However, the initial conditions at larger spatial scales are still the subject of considerable debate.
Broadly, there are two models for the formation of massive stars: core accretion and competitive accretion (Tan et al. 2014).In core accretion, the initial fragmentation of a clump creates gravitationally bound cores with a range of masses, with massive stars forming from massive cores with surface densities larger than ∼ 1 g cm −2 (McKee & Tan 2003).There is a direct correlation between the mass of an individual core and the mass of stars that form from the core with multiplicity depending on whether or not there is fragmentation in the gravitationally unstable disk (Krumholz et al. 2009).Thus, the matter that constitutes a massive star is derived primarily from that of the parent core.It also follows that the the stellar initial mass function is related to the core mass function (Alves et al. 2007).
In competitive accretion, a massive clump fragments into a large number of Jeans mass fragments followed by competitive accretion of matter from the clump reservoir, with the process naturally leading to the stellar initial mass function (Bonnell et al. 2003(Bonnell et al. , 2004)).Thus, most of the matter that constitutes the massive star is derived from the entire clump rather than an individual core/fragment (Smith et al. 2009).The process of competitive accretion also leads to the massive star being located at the center of the gravitational potential of the clump (Bonnell et al. 2001) and naturally results in the formation of massive stars in clusters.
One of the observational techniques to discriminate between the theories is to trace the distribution of mass from large clump scales ( 0.5 pc) to that of cores ( 0.1 pc).Historically, observational studies targeting the early stages of massive star formation have faced challenges due to their rarity, rapid evolutionary time-scale, and deeply embedded clustered environment.However, the completion of blind surveys of the Galactic plane from millimeter to near-infrared wavelengths has led to a significant advancement in the ability to identify massive star forming regions at different evolutionary stages.
There have been numerous studies targeting the fragmentation of massive clumps at various scales, with the number of studies being too numerous to list comprehensively.These studies reveal a wide diversity in the fragmentation properties with some clumps undergoing limited fragmentation (e.g.Zhang et al. 2009;Bontemps et al. 2010;Palau et al. 2013;Zhang et al. 2015;Csengeri et al. 2017;Beuther et al. 2018;Figueira et al. 2018), while other sources manifest significant fragmentation (e.g.Bontemps et al. 2010;Palau et al. 2013;Cyganowski et al. 2017;Beuther et al. 2018;Sanhueza et al. 2019;Svoboda et al. 2019).Since this diversity provides evidence both in favor of and against both theories, there is still significant debate regarding the mechanism leading to the formation of massive stars.
In this paper, we report on the results of a study targeting the fragmentation of 16 massive clumps spanning a range of evolutionary stages using the Submillimeter Array (SMA) in its compact configuration.We describe the sample selection process and sample characterization in Section 2. The SMA observations and data reduction are described in Section 3, while Section 4 outlines the results, methodology for identifying individual fragments and measuring their properties along with notes on selected sources.A discussion of the fragmentation observed along with the implications for the formation process is presented in Section 5. Finally, a summary of the results obtained are listed in Section 6.

SAMPLE SELECTION
Our sample of massive clumps were selected from the Apex Telescope Large Area Survey of the Galaxy at 870 µm (ATLASGAL; Schuller et al. 2009).We selected a total of 16 sources (see Table 1) spanning a range of evolutionary states, inferred using their properties at mid-infrared wavelengths (Troost, T., Masters thesis).We ensured that the sources could be conveniently grouped on the sky so that three to four sources could be targeted in a single observing track with the SMA.
In order to characterize the large scale properties of the clumps targeted in this study, we constructed the spectral energy distributions (SEDs) of the sources by combining the 870 µm ATLASGAL data with 500, 350, 250 and 160 µm data from the Herschel infrared Galactic plane survey (Hi-GAL; Molinari et al. 2010).We did not use the 70 µm data from Hi-GAL since the emission at this wavelength has a significant contribution from warm dust, especially in the more evolved sources.We carried out photometry using the IDL based hybrid photometry and extraction routine (HYPER; Traficante et al. 2015).The variable background in the Hi-GAL maps is modelled as a 2D polynomial followed by a 2D Gaussian fit to the source.In case of crowded fields, a multi-Gaussian fit is performed for de-blending.The unique feature of multi-wavelength photometry using HYPER is that it uses the same aperture, specified at a reference wavelength, for photometry in all bands.This ensures that the same volume of gas and dust is included for flux determination at all wavelengths.We used a reference wavelength of 250 µm for defining the aperture since the emission morphology and resolution is comparable to that of the ATLASGAL data (18 ′′ at 250 µm versus 18.2 ′′ at 870 µm) while the signal to noise ratio is significantly better (see Appendix A for more details).The HYPER flux densities are generally comparable to the Hi-GAL catalogue fluxes except when the sources are affected by blending (especially at 500 µm).The 870 µm flux densities determined using HYPER are also consistent with the ATLASGAL compact source catalog (Csengeri et al. 2014) determined using the GAUSSCLUMPS algorithm (Stutzki & Guesten 1990).The dust temperature (T d ) was determined by fitting the SED with a single component grey body where Sν is the flux density, Ω is the solid angle of the source, Bν is the black body function, τ0 is the dust optical depth at reference frequency ν0 and β is the dust emissivity which was allowed to vary between 1 and 3. Fig. 1 shows an example of SED fits for the sources in our sample.The source AG36.79-0.20,which is dark at 70 µm, is affected by the presence of a second source that is prominent at wavelengths shorter than 250 µm.Since it was not possible to de-blended the two sources, the SED was fit to the 870 -350 µm data using the constraint that the fit should give a 70 µm flux density lower than the completeness limit of the Hi-GAL survey.The isothermal clump mass was determined from the 870 µm flux using where R is the gas to dust ratio that is assumed to be 100, d is the distance to the source, and κν is the dust opacity which is taken to be 1.85 cm 2 g −1 (Ossenkopf & Henning 1994).We have used the kinematic distance for all sources except for those in the W51 region (AG49.25-0.41 and AG49.60-0.25)for which we have adopted the distance measured using trigonometric parallax towards H2O masers (Sato et al. 2010).The kinematic distances were determined using the Galactic rotation curve prescribed by Reid et al. (2014) adopting the radial velocity of NH3 (1,1) (Wienen et al. 2012).The kinematic distance ambiguity was resolved by looking for HI absorption (if radio continuum is associated with the source or region) or self-absorption (if radio continuum is absent) using data from the VLA Galactic Plane Survey (VGPS; Stil et al. 2006) and the Boston University -Five College Radio Astronomy Observatory Galactic ring survey (GRS; Jackson et al. 2006).The resolution of the kinematic distance ambiguity using HI self-absorption is challenging since the natural line profile of HI can be mistaken for an absorption feature.To avoid incorrect resolution of the distance ambiguity, we compared the morphology of 13 CO emssion in GRS with that of HI and only classified a source at the near distance if the morphology of absorption in HI matched with that of 13 CO emssion.Fig. 2 shows examples of sources determined to be at near and far distances using this technique.In the case of source AG47.03+0.24(Fig. 2), while there appears to be a HI self-absorption feature at the source centre, the morphology of the absorption feature does not match with that of the CO emission showing this source to be at the far distance.
The bolometric luminosity of the clumps can be determined from the total flux by integrating the SED.The SED above from sub-millimetre to far-infrared wavelengths is sufficient to characterize sources that are in very early stages of evolution (class A sources) that are dark at 70 µm.However, sources at later evolutionary stages have significant emission up to mid-infrared or even near-infrared wavelengths.In order to determine the bolometric luminosity of these sources, we extended the SED to shorter wavelengths using data from MIPSGAL, MSX (Price et al. 2001), GLIMPSE (Benjamin et al. 2003) and 2MASS (Skrutskie et al. 2006) surveys.We then fit the SEDs using the young stellar object (YSO) models of Robitaille et al. (2007).The objective of this fitting is only to determine the bolometric luminosity of the source and is hence not expected to be sensitive to the precise details of the models themselves.
The effective radius of the source is computed from the HYPER 2D Gaussian fit using where A is the deconvolved area, θmaj and θmin are the full width at half maximum (FWHM) major and minor axes respectively and θ b is the FWHM beam width of the telescope (18 ′′ for the Hi-GAL 250 µm data).The surface density is determined from the mass and effective radius using The gravitational stability of the clumps against collapse can be inferred from the virial parameter which is defined as where σ is the velocity dispersion.We have calculated the virial parameter using the velocity dispersion measured from the GRS 13 CO data.The 13 CO spectrum is extracted at the ATLASGAL source position (listed in Table 1) and fit with a Gaussian to determine the velocity dispersion.Since 13 CO is a trace molecule, the velocity dispersion that characterizes the total column of gas (Fuller & Myers 1992) is computed using where T is the kinetic temperature of the gas and is assumed to be the dust temperature (T d ) and µ is the mean molecular weight of the medium and is assumed to be 2.35 (derived from the cosmic abundance scale of Przybilla et al. 2008).
Table 1 shows that the targeted sample of ATLASGAL clumps are potential sites for forming massive stars.The clump masses are typically significantly higher than 100 M⊙ and the surface densities are well above 0.05 cm 2 g −1 , which is considered as an empirical lower limit for forming massive stars (Urquhart et al. 2013).The clump masses and effective radius also satisfy the empirical relation of Kauffmann et al. (2010) highlighting their potential to form massive stars.The virial parameter ranges from 0.68 to 1.71 showing that almost all sources are unstable to gravitational collapse.The dust temperature ranges from 10.5 to 40.5 K showing that the sources span a range of evolutionary states.
One can obtain a more quantitative estimate of evolutionary state of the clumps from their mass to light ratio.Using the turbulent core model of McKee & Tan (2003), Molinari et al. (2008) derived evolutionary tracks for massive young stellar objects (MYSOs) in the L bol − Menv diagram.Although the tracks have been derived from the turbulent core model, the overall features are expected to be model independent.This is because accretion with minimal loss of the envelope mass is expected to be the dominant process prior to the central object reaching the zero age main sequence (ZAMS), while gas dispersion or envelope clearing is expected to dominate the later phases.This results in vertical evolutionary tracks in the L − M diagram during the accretion phase with the luminosity dominated by accretion luminosity, while the late phase is characterized by horizontal tracks towards lower masses since the luminosity is expected to be constant post ZAMS.Based on a study of MYSOs in 42 regions of massive star formation, Molinari et al. (2008) found that the massive sources with SEDs consistent with the presence of an embedded ZAMS star (called "IR-P" sources) were roughly at the transition between the accretion and clearing phase (see Fig. 9 of Molinari et al. 2008).We thus fit the bolometric luminosity as a function of envelope mass for the "IR-P" sources of Molinari et al. (2008) and use this fit to define an evolutionary state parameter The value of unity for ǫ defines the fit of luminosity vs. mass for the "IR-P" sources, with ǫ < 1 in the accretion phase and ǫ > 1 in the clearing phase.One thus expects ǫ to increase with time as the source evolves.We also note that while one can use (L/M ) directly as a marker for evolutionary stage, we prefer to use ǫ since it is based on the models developed by Molinari et al. (2008).Fig. 3 shows the clumps targeted in this study in the L − M diagram with the evolutionary state parameter being tabulated in Table 1.The value of ǫ ranges from 0.012 to 4.3.This clearly demonstrates that the clumps span a wide range of evolutionary states with a majority of the sources being in the accretion phase.The source sample is thus well suited to address the question of mass assembly in massive star forming regions.

OBSERVATIONS
The observations towards the 16 ATLASGAL sources were carried out using the SMA in its compact configuration between 2010 May 5 and 2011 October 18.Two sources, AG36.90-0.41 and AG41.05-0.25 were followed up with the SMA in its extended configuration on 2011 August 16.The details of the sources observed and calibrators used are tabulated in Table 2.
All the observations were carried out using the 345 GHz receiver and the ASIC correlator.The correlator was used in the 1 receiver 4 GHz mode which provides a bandwidth of 4 GHz per sideband.The local oscillator (LO) tuning of the 2010 observations was such that rest frequencies of 333.5 to 337.5 GHz were covered in the lower sideband (LSB) while the upper sideband (USB) covered rest frequencies from 345.5 to 349.5 GHz.The 2011 observations were set up such that rest frequencies of 342 to 346 GHz were covered in LSB, while the USB covered rest frequencies from 354 to 358 GHz.In all observations, the correlator was set up with a uniform channel spacing of 0.8125 MHz (∼ 0.7 km s −1 ) across the entire band.
The data were calibrated using the IDL based MIR package.After initial inspection and flagging, the data were weighted by the system temperatures, which ranged from 200 to 600 K depending on elevation.The data were then calibrated for baseline dependent passbands, antenna and time dependent gain and flux using appropriate calibrators.When Neptune was used for flux calibration, the zero spacing flux was used rather than adopting a black body model due to the planet being affected by CO absorption lines.The flux calibration is estimated to be accurate to about 15 per cent.The continuum was then constructed from the line-free channels independently for each sideband.The data were then imported to CASA (version 4.2) and imaged a Noise obtained by imaging compact and extended configurations together using the task 'clean' using multi-frequency synthesis.The uv−range of the SMA in the compact configuration ranges from ∼ 10 − 70 m, providing a synthesized beam of around 2 ′′ , with the 1σ noise ranging from 1.8 to 6.8 mJy beam −1 (see Table 2).The relatively high noise values in some sources is due to the presence of strong emission where the data were limited by dynamic range.The uncertainty in absolute astrometry is estimated to be 0.5 ′′ by comparing the positions of the gain calibrators with their catalog positions.
The data from the extended array, which have a uv−range of ∼ 25 − 230 m, were combined with that of the compact array for both AG36.90-0.41 and AG41.05-0.25 and then imaged.This preserved the total flux density observed in the compact array, but provided the angular resolution of the extended array, which was around 0.6 ′′ × 0.9 ′′ .The shortest baseline of the compact array observations was typically lower than 15 kλ, implying that the observations are sensitive to a largest angular scale of ∼ 13 ′′ .The beam size of an SMA antenna is 25.5 ′′ at the 20% level and ∼ 30 ′′ at the 10% level, which sets the field of view of the array.

RESULTS
We detected continuum emission from the SMA observations towards all sources at greater than 5σ level.Figure 4 shows an example of the 870 µm image from ATLASGAL alongside that from the SMA in its compact configuration.The full set of figures are available in Appendix B. Clearly, the maps at high resolution show a range of morphologies, with some sources fragmenting into many cores, while others are dominated by a single core.
The SMA images reveal a diversity in morphology that is a direct consequence of the hierarchical structure in molecular clouds.The hierarchical structure can be analyzed statistically (e.g. through principal component analysis; Heyer & Brunt 2004) or through segmentation by dividing the data into different structures.Two prominent algorithms used for the latter analysis are CLUMPFIND (Williams et al. 1994) and GAUSSCLUMPS (Stutzki & Guesten 1990;Kramer et al. 1998).An alternate method to describe hierarchical structure is through dendrograms (Houlahan & Scalo 1992;Rosolowsky et al. 2008).The dendrogram technique identifies local maxima as leaves which merge into a larger structures at nodes.The two main parameters that control the identification of these structures are (a) the height of a leaf above a node (i.e. the difference between the peak of a local maximum and the level at which it merges into a larger structure) and (b) the number of image pixels that are larger than its neighbors over a specified box size.The two parameters are used to control spurious local maxima that result from noise in real data.
We have carried out the dendrogram analysis for all the sources in our sample through the following approach.First, we determine local maxima that are significant at the 5σ level.Next, contiguous regions were identified down to ∼ 1σ around each local maximum.Each distinct contiguous region was defined as a large scale "core", and corresponds to a root in the dendrogram.All local maxima in a root that have a peak intensity that is at least 3σ, and are at least 1σ above a node (i.e. the  level at which local maxima merge into a larger structure), were defined as leaves.Figure 5 shows an example of the hierarchical structure in AG49.25−0.41 through a dendrogram.While the dendrogram analysis is useful to examine the hierarchical structure in the clumps, one cannot directly use the analysis to determine the masses of the fragments.This is because the dendrogram only gives the height of a leaf above that of a node (i.e. the height of a peak above a saddle point).This by itself is not an estimate of the total flux density associated with the leaf.One would thus have to carry out an analysis using CLUMPFIND or GAUSSCLUMPS to determine the flux density and the mass associated with each fragment.
Both CLUMPFIND and GAUSSCLUMPS have merits and demerits with regards to determining the properties of the fragments.CLUMPFIND identifies isolated contours as individual clumps and blended contours that surround more than one local maximum are split using a "friends-of-friends" algorithm (Williams et al. 1994).In sources with several peaks and diffuse extended emission, this can give rise to fragments with odd shapes.An additional factor to be kept in mind is that the output from CLUMPFIND is at times a sensitive function of the choice of contour levels that are provided as an input to the algorithm.For our work, we mitigate this issue by using the flux densities of nodes that are identified in the dendrogram analysis as the input contour levels for CLUMPFIND.Since each node is a saddle point at which two local maxima merge together, the flux density of a node is the appropriate threshold contour level that can separate the two local maxima into individual fragments.This also ensures that spurious fragments are not identified due to closed contours that may arise due to noise in the image.Figure 6 shows the output of CLUMPFIND for the source AG49.25−0.41where 6 distinct fragments have been identified, some with odd shapes.On the other hand, GAUSSCLUMPS fits the emission from the source using multiple Gaussian components.While this works well when the emission is comprised of several compact sources (which may be blended), difficulties are encountered when the overall morphology is complex with a mix of compact and extended emission.In such cases, some Gaussian components can be introduced purely for the sake of fitting the overall emission structure even though they may not be bona-fide fragments.This may in addition result in the flux densities of some fragments being under-estimated since part of the emission may be fit by an extended component.Considering the complex emission morphology in some sources, we choose to use CLUMPFIND in tandem with the dendrogram to identify fragments and determine their properties (see Table 3).It should be noted that both the dendrogram and CLUMPFIND analysis were first carried out using the CLEAN images which have uniform noise across the image, after which the source properties were extracted from the images that have been corrected for the primary beam response of the SMA antenna.
The masses in Table 3 have been calculated using the same temperature as for the ATLASGAL clump (Table 1).However, it is likely that the temperature at the scales probed by the SMA is larger than what is observed at the clump scale of ATLASGAL.An extreme example is the case of the hot core G351.77−0.54 which was observed to have brightness temperatures close to 1000 K at a scale of 18 − 40 AU (Beuther et al. 2018).Since our observations are at scales of ∼ 3000 − 20, 000 AU, we do not expect to see such high temperatures.Nonetheless, it is likely that the temperatures are higher than what is seen at 18 ′′ resolution, especially for sources with ǫ 1.To explore this in more detail, we have looked at earlier studies of massive star forming regions at high resolution.Wiseman & Ho (1998) studied the Orion ridge region using observations of ammonia with the Very Large Array (VLA) at a resolution of 8.5 − 9.0 ′′ (which corresponds to physical scales comparable to that probed in our study since our sources are at much larger distances), and found temperatures to range from ∼ 15 − 30 K in filaments although values as high as ∼ 80 K are observed in the central core.Issac et al. (2020) studied the extended green object G19.88−0.53 using ALMA and found temperatures ranging from 47 − 116 K, with high temperatures being characteristic of hot cores which correlates with the presence of an ionized radio jet in the region.On the other hand, Wang et al. (2012) studied the region G28.34+0.06 using ammonia observations with the VLA at ∼ 2 ′′ resolution and found temperatures to range from 8 − 30 K with higher temperatures coinciding with heating induced by outflows impinging on molecular gas.Based on this, we expect high temperatures in the sources AG48.61+0.02 and AG50.28−0.39,which have ǫ > 1 and L > 10 5 L⊙ with accompanying significant internal heating.Although these sources also have the highest dust temperatures at 18 ′′ resolution, the temperatures at high resolution may be significantly higher.An increase in temperature by a factor of 2 will reduce the masses of the cores by a similar factor.Thus, the masses of the massive fragments associated with high luminosity ATLASGAL clumps should be taken as upper limits, with the true masses being a factor of 1.5 to 2 smaller depending on the extent of central heating in the sources.However, considering the results of Wiseman & Ho (1998); Wang et al. (2012), we do not expect the masses of fragments that have no counterparts at mid-infrared wavelengths to be very different from what is reported in Table 3. Note-The columns show (1) the source name (2)−(3) J2000 equatorial coordinates (4) peak flux density (primary beam corrected) (5) integrated flux density (primary beam corrected) (5) mass ( 6) deconvolved solid angle (7) effective radius and (8) mass surface density.Columns ( 9), ( 10) and ( 11) indicate the presence or absence of 70 µm counterpart in the Hi-GAL survey, 24 µm counterpart in the MIPSGAL survey ( † indicates saturation) and an MIR counterpart in the GLIMPSE survey respectively.
a Not a point source; possible Extended Green Object (EGO) 4.1.Notes on selected sources AG36.90−0.41 and AG41.05−0.25 -Both sources are dominated by a single structure at ∼ 2 ′′ resolution, although some extended emission is seen surrounding the central source.Figure 7 shows the sources at 0.6 ′′ resolution by combining the compact and extended configurations of the SMA.Clearly, the emission continues to be dominated by a single core even at high resolution.Both sources are associated with a 24 µm point source with AG36.90−0.41 also hosting a 6.7 GHz methanol maser (Pandian et al. 2011).AG48.61+0.02-There are two fragments in this source, with SMA1 dominating the emission with a mass of 393 M⊙ (although the true mass may be much lower depending on the temperature of the core at this scale).The other fragment, SMA2 has a mass of only 21.8 M⊙.The infrared emission from the clump is also very strong with the entire clump being saturated in the MIPSGAL survey.The morphology of emission in the GLIMPSE survey is suggestive of there being an outflow along with extended green emission (i.e.excess emission at 4.5 µm).
AG53.14+0.07-There are seven fragments detected in this clump with masses ranging from 0.7 M⊙ to 19 M⊙.The clump is also associated with a 6.7 GHz methanol maser whose location is closest to that of SMA2 (Pandian et al. 2011).SMA2 is also associated with a point source in both GLIMPSE and MIPSGAL surveys although no point source is listed in the catalogs due to possible saturation.

Nature of fragments
We have detected a total of 43 fragments in 16 clumps.We have searched for counterparts at mid-infrared wavelengths at wavelengths shorter than 70 µm using the Hi-GAL, MIPSGAL and GLIMPSE catalogs using a search radius of 3 ′′ (see Table 3).We find that 8 fragments have point source counterparts in the GLIMPSE catalog.An additional three fragments have emission from a point source in the GLIMPSE maps although they are not listed in the catalog.The infrared colors in the Spitzer-IRAC bands are consistent with 5 out of 8 fragments being Class I young stellar objects (YSOs) and one source (AG410.7−0.12SMA3) being a Class II YSO (Gutermuth et al. 2008).The nature of the two other GLIMPSE sources is not clear since they are detected in only two IRAC bands.
Seven fragments have a counterpart in the MIPSGAL catalog, with an additional two sources having a point source that is saturated in the central region.We also find that the 24 µm emission in MIPSGAL is fully saturated around AG48.61+0.02 and AG50.25−0.39.Eleven fragments have associated 70 µm counterparts in the Hi-GAL survey.Thus, around 28 fragments have no counterpart at mid-infrared wavelengths and may be starless, although further studies using deeper observations are required to establish the same.Of these, AG36.79 SMA1, AG36.79 SMA2, AG47.03SMA1, AG47.05 SMA1, AG52.82 SMA2 and AG52.92SMA2 have masses exceeding 20 M⊙, and may be suitable candidates for being massive starless cores.

Fragmentation of the ATLASGAL clumps
The fragmentation properties and the underlying processes that lead to fragmentation is one of the outstanding questions in the study of massive star formation.According to the Jeans instability criterion, fragmentation occurs when the mass exceeds the Jeans mass (MJ ) over a Jeans length (λJ ).The Jeans length for thermal fragmentation is given by where c 2 s = kT /µmH is the sound speed and ρ is the mean density in the clump.The Jeans mass is given by The Jeans length and Jeans mass for fragmentation that is governed by turbulence rather than thermal motion can be derived by replacing the isothermal sound speed with the turbulent velocity dispersion Table 4 shows the thermal and turbulent Jeans length and Jeans mass for our sample along with the range of nearest neighbor separations for the fragments detected with the SMA.It is to be noted that these separations are observed in projection and are thus lower limits of the true inter-fragment separation.We have used the dust temperature (see Table 1) to estimate the thermal sound speed, while the mean density has been computed from the clump mass and effective radius using Figure 8 shows a histogram of the ratio of the nearest neighbor separations to the thermal and turbulent Jeans length.
Table 4 and Figure 8 show that the separation and masses of the fragments detected in our study are more consistent with thermal rather than turbulent fragmentation.The masses of all fragments are smaller than the turbulent Jeans mass, and the inter-fragment separations are much smaller than the turbulent Jeans length except for one source.Figure 8 shows that there are cases where the nearest neighbor separation is also smaller than the thermal Jeans length, similar to what was observed by Beuther et al. (2018).However, all these cases pertain to fragments that are part of a larger structure, but are identified as a distinct leaf in the dendrogram.The sub-Jeans scale separations in these cases are likely to be a consequence of the local clump densities being higher than the mean density used to compute the thermal Jeans length.Though less likely, it is also possible that some fragments are extended structures pertaining to the main fragment rather than being independent fragments themselves.The results obtained here are in agreement with several statistical and detailed studies to date (e.g.Palau et al. 2015;Beuther et al. 2018;Palau et al. 2018;Beuther et al. 2019;Sanhueza et al. 2019;Svoboda et al. 2019;Saha et al. 2022) where the fragmentation has been found to be governed by thermal rather than turbulent motion.
One can also observe significant diversity in the fragmentation properties of the sample.Five clumps are dominated by a single massive core with masses ranging from 38 to 276 M⊙ with no low-mass cores to the sensitivity limit.Of these, the observations are sensitive to the thermal Jeans mass at the 5σ limit for three sources suggesting a real paucity of low-mass cores.Further, two sources continue to be dominated by a single fragment down to 0.6 ′′ resolution (∼ 0.025 pc scale).Another four clumps have the lowest fragment mass to be greater than 5 M J,th , with the number increasing to nine clumps at the limit of 3 M J,th , although one has to bear in mind the caveat that the observations are not sensitive to the thermal Jeans mass limit in some cases.This is similar to the observations of Zhang et al. (2009); Bontemps et al. (2010); Palau et al. (2013); Zhang et al. (2015); Csengeri et al. (2017); Figueira et al. (2018) where limited fragmentation was seen with the most massive fragment having masses several times that of the thermal Jeans mass.On the other hand, two clumps display significant fragmentation with the lowest fragment mass comparable to the thermal Jeans mass.This is similar to what was observed by Cyganowski et al. (2017); Henshaw et al. (2017); Beuther et al. (2018); Zhang et al. (2021) where a significant population of low-mass fragments were found to form coevally with massive fragments.We also do not see any correlation between the degree of fragmentation, quantified by the number of clumps, and the mean density.Since the Jeans mass decreases with increasing density, one expects to observe a higher degree of fragmentation with increasing density.Thus, although the fragment separations are closer to that expected from thermal fragmentation, there are other factors that control the degree of fragmentation.As pointed out by Palau et al. (2014); Beuther et al. (2018), the diversity in fragmentation may be due to the density profile of the parent clump, with steeper density profiles seen to limit fragmentation.We plan to carry out a study of the density profiles of our ATLASGAL sample to verify whether this can explain our observed fragmentation diversity.Alternately, magnetic fields may also play an important role in the degree of fragmentation (Commerçon et al. 2011;Fontani et al. 2018).For example, radiation magnetohydrodynamic simulations of Commerçon et al. (2011) are able to inhibit fragmentation in highly magnetized cores, with the level of fragmentation dependent on the strength of the magnetic field (Fontani et al. 2018).

The initial fragmentation of massive clumps
One of the primary differences between the theories of core accretion and competitive accretion is regarding the initial fragmentation of massive clumps.While the theory of core accretion predicts initial fragmentation into a population of cores whose mass function mimics that of the stellar initial mass function, the competitive accretion scenario has the clump fragmenting into a large number of fragments with mass comparable to the thermal Jeans mass.Observational studies obtain a snapshot of the state of fragmentation well after the initial stages, with it being difficult to establish the period of time that has elapsed since the beginning of fragmentation to the time of observation.It is hence not straightforward to make interpretations about the initial conditions of fragmentation based on observational results.
However, the observation of sources that are dominated by a single core (e.g.CygX-N63 of Bontemps et al. 2010, NGC 7129-FIRS2 andI20126+4104 of Palau et al. 2013, several ATLASGAL clumps in our work, and other works cited in the references above), and the paucity of low-mass fragments in some clumps (e.g.Zhang et al. 2015) suggests that the parent clumps in these cases did not fragment into a large number of Jeans mass fragments.This does not mean that massive stars form in isolation in these sources since the core mass is typically a small fraction of the clump mass, and low-mass stars may form in the clump at later times (Zhang et al. 2015).An interesting observation in the context of initial fragmentation is the offset between the locations of the SMA fragments and that of the ATLASGAL clump.We see four sources (AG36.83−0.04,AG41.07−0.12,AG47.05+0.25 and AG52.92+0.41)where the offset exceeds 5 ′′ .Figure 9 shows the offsets for AG41.07−0.12 and AG47.05+0.25 respectively, with all three fragments in the former offset from the ATLASGAL source in excess of 13 ′′ , whereas the latter has a single fragment that is offset from the mean clump position by over 6 ′′ .A caveat in the context of this discussion is that the position of the ATLASGAL clump depends on the algorithm used to define the source.For example, the positions of ATLASGAL sources differ by several arcseconds between the compact source catalog (Contreras et al. 2013) and the catalog generated by GAUSSCLUMPS (Csengeri et al. 2014).However, based on visual inspection, the positions of the SMA fragments are significantly offset from the mean clump position in all the cases above.Thus, the offset between the locations of the fragments and the large scale clump is real, subject to the astrometric accuracy of the two data sets.
The astrometry of the ATLASGAL data have been found to be accurate to that of the pointing of the APEX telescope which is ∼ 2 − 3 ′′ .The astrometric accuracy of the SMA data is determined by the uncertainty in the positions of the gain calibrators, separation between the calibrators and target sources, uncertainty in the antenna positions and the signal to noise ratio of the target source itself.The positions of the gain calibrators are from VLBA calibrator surveys which have position acccuracies better than 1 milli-arcsecond (mas), while antenna positions are known to 0.1−0.2λ(where λ is the wavelength).As shown in Table 2, the SMA observations included two gain calibrators on all days of observations.For example, the 2010 October 4 data included 1751+096 and 2025+337 as gain calibrators, which are located about 26 • from G47.05+0.25 respectively.The positions of 1751+096 and 2025+337 are known to better than 0.2 mas (Johnston et al. 1995) and 1.1 mas (Beasley et al. 2002) respectively.The uncertainty in astrometry from the phase transfer from gain calibrators to target sources can be estimated by examining the phase solutions over a time duration corresponding to the separation between the calibrators and targets (1.3 hours for 20 • separation).An examination of the gain solutions shows that the astrometric uncertainty from this process is smaller than the size of the synthesized beam (∼ 2 ′′ ).To further verify the accuracy of astrometry, we calibrated the data using only one gain calibrator, 1751+096, and imaged the second gain calibrator, 2025+337, which is almost 43 • away.The position of 2025+337 agreed with the value given in the VLBA calibrator survey to within 0.5 ′′ .Since our target sources are much closer to the gain calibrator compared to the separation between the two gain calibrators above, we infer that the astrometric accuracy of the SMA fragments are better than 0.5 ′′ .Thus, the detected offsets between the positions of the ATLASGAL and SMA sources is larger than the astrometric uncertainty of the source positions and is clearly significant.
Finally, one has to consider whether the offsets are an outcome of differences in resolution.However, a difference in resolution by itself cannot introduce an offset between a compact structure (that is unresolved at poorer angular resolution) and the overall large scale structure.Alternately, a smoothing operation of an SMA dataset to the lower resolution of ATLASGAL cannot by itself change the location of the peak emission unless there are multiple compact sources of comparable strength present, in which case the smoothed data would show a peak at the mean position of the multiple sources.In cases where there is a single compact source, or where one source is significantly stronger than other sources in the vicinity, the smoothed data must have a peak emission that is identical to or close to the location of the strong compact source.Thus, the detected offsets between the ATLASGAL emission and our SMA maps are unlikely to be due to the effect of resolution alone.Considering that observations using the compact configuration of the SMA filter out emission on scales larger than ∼ 13 ′′ (see Section 3), the observed offsets are most likely depicting the actual differences between the location of smoother diffuse clump emission that is measured by ATLASGAL and the compact cores that are discernible upon filtering of this extended emission by the interferometer (for example, in the case of AG47.05+0.25, over 90% of the emission is filtered out by the SMA).
The angular offsets that exceed 5 ′′ correspond to a physical projected separations that are a significant fraction of the effective radius of the clump (for example, 0.57 pc in the case of AG41.07−0.12 and 0.22 pc in the case of AG47.05+0.25).If the process of forming a massive core involves competitive accretion of the gas reservoir onto a large number of Jeans mass fragments, the massive core must be located near the center of the gravitational potential (Bonnell et al. 2001).Since fragments compete among themselves to accrete the gas that is distributed within the clump, fragments that are far from the center of the potential do not accrete significantly and retain low masses (Bonnell & Bate 2006).The observation of sub-millimeter fragments that are significantly offset from the center of the gravitational potential of the clump in at least four sources strongly suggests that the cores did not form from competitive accretion of the clump material onto a large number of fragments throughout the clump.We stress that these results primarily suggest that a clump does not undergo initial fragmentation into a large number of Jeans mass fragments as invoked by models of competitive accretion.It does not rule out the phenomenon of competitive accretion to transfer mass from clumps to multiple cores within that may have formed through limited fragmentation.

Mass assembly from clumps to cores
The "core formation efficiency" (CFE) is defined by Bontemps et al. (2010); Palau et al. (2013); Csengeri et al. (2017) as the fraction of the clump mass that is in compact structures.We have estimated the CFE for each source by taking the ratio of the total flux density of all fragments detected in the SMA map to that in the ATLASGAL map.The latter is computed by taking the total flux density that is within the 10% level of the SMA primary beam centered at the SMA pointing location.It is important to note that estimating the CFE as the ratio of flux densities is based on the assumption that the temperatures at high and low resolution are identical.As noted in Section 4, the actual temperature at high resolution may be significantly higher, especially in sources with ǫ > 1.This results in the CFE to be overestimated in such sources.To account for this effect, we incorporate an additional uncertainty in temperature by a factor of 2 for sources with ǫ > 1 and a factor of 1.5 for sources with ǫ > 0.1 for estimating the CFE.
Figure 10 shows the CFE as a function of ǫ (eq.1), the evolutionary state parameter (left panel) and the mean number density of H2 (right panel) in the ATLASGAL clump.We observe a positive correlation between the CFE and ǫ, with a Spearman correlation coefficient of 0.69 and a p−value of 0.003.This shows that the two quantities are correlated with only a 0.3% probability of the correlation occurring by chance.It is to be noted that the correlation coefficient does not take into account the uncertainties associated with the two quantities.However, as seen from Figure 10, accounting for the uncertainty in CFE only reduces the correlation coefficient and doesn't remove it altogether.In fact, considering the CFE to be at the lower limits implied by the uncertainty only reduces the correlation coefficient to 0.31.Thus, we find that compact structures encompass a larger fraction of the clump mass as the source evolves.This in turn shows that mass is channeled from clumps to cores over time and strongly suggests that the evolution of a core does not happen in isolation as in models of core accretion (e.g.Krumholz et al. 2009).We also observe a weak correlation between the CFE and the density of the ATLASGAL clump, with a Spearman correlation coefficient of 0.49 and a p−value of 0.06.This is similar to the observations of Bontemps et al. (2010); Palau et al. (2013); Csengeri et al. (2017), although our data does not adequately sample H2 number densities greater than 10 6 cm −3 .
We also find a good correlation between the maximum mass of a fragment (Mmax) in a clump and the bolometric clump luminosity (left panel of Figure 11, with the uncertainty in maximum mass incorporating the heating effect as outlined in the previous paragraph).A Spearman correlation test gives a correlation coefficient of 0.81 with a p−value of 1.4 × 10 −4 .Since both the mass and luminosity are dependent on the distance, we also carried out a partial Spearman correlation test to remove the effect of distance.This gives a correlation coefficient of 0.6 which indicates that the two quantities are well correlated.The degree of correlation seen in our work is much stronger than that of Palau et al. (2013) who observed only a weak correlation in their sample.As indicated by Palau et al. (2013), the presence of this correlation argues against a formation scenario purely by competitive accretion since such a mechanism precludes any correlation between the mass of small scale condensations and the final stellar mass (Bonnell et al. 2004).We also do not see any positive correlation between the maximum mass of a fragment and the number of fragments (right panel of Fig 11).The correlation is in fact negative with a correlation coefficient of −0.65 with a p−value of 0.006, whereas Bonnell et al. (2004) predict an increase in the number of stars with increasing maximum mass.Since competitive accretion models predict the number of stars to increase with the maximum mass in the cluster, the lack of a positive correlation between the number of fragments and the maximum mass of a fragment also suggests that massive stars do not form purely by competitive accretion.The transfer of mass from clump to core scales with evolution coupled with evidence against pure competitive accretion suggests that massive stars grow through clump fed core accretion (e.g.Wang et al. 2010).This is also corroborated by observations such as that of Zhang & Wang (2011); Yuan et al. (2018);Peretto et al. (2020) and references therein.Competitive accretion is likely to play a role in the process of mass transfer from clump to core scales when multiple cores are present in the clump.

CONCLUSIONS
We have carried out a study of fragmentation in massive clumps using observations of ATLASGAL sources with the SMA in its compact configuration at an angular resolution of 2 ′′ .Our main conclusions are 1.We observe a wide diversity of fragmentation in the sample.Some clumps are dominated by a single core or have limited fragmentation with the least massive fragment having a mass several times the Jeans mass, while other clumps show significant fragmentation with masses down to the Jeans mass.Two sources continue to be dominated by a single core even when observed with the extended configuration of the SMA at a resolution of 0.6 ′′ , or a physical scale of ∼ 0.025 pc.
2. The separation and masses of the fragments are smaller than the Jeans length and Jeans mass inferred from the turbulent velocity dispersion.Thus, the data are more consistent with thermal rather than turbulent fragmentation.
3. Four sources show a significant projected separation between the clump and the cores within, which taken together with the observation of limited fragmentation in many sources suggests that clumps do not undergo initial fragmentation into a large number of thermal Jeans mass fragments.
4. We find evidence for a larger fraction of clump mass to be incorporated into cores with evolution.
5. We also find a good correlation between the mass of the most massive fragment in the clump and the bolometric luminosity.We do not see the number of fragments to increase with increasing mass of the most massive fragment.These results suggest that massive stars do not form purely by competitive accretion.
Our results support a broad picture of massive star formation wherein the initial fragmentation of massive clumps is governed by thermal fragmentation in combination with the underlying density distribution and/or magnetic fields.Subsequently, cores continue to gain mass from the clump, a process that may incorporate aspects of competitive accretion depending on the number of fragments within the clump.Individual massive stars or binary systems are likely form within these cores, as seen in high resolution observations of Beuther et al. (2017); Motogi et al. (2019).

Figure 2 .
Figure 2. The panels show HI emission from VGPS at the radial velocity of the source with 13 CO contours from GRS overlaid in blue.The left and right panels show examples of sources with and without HI self-absorption respectively.

Figure 3 .
Figure 3.The L − M diagram of the sample.The sources targeted in this work are shown in triangles while the 'IR-P' sources of Molinari et al. (2008) are shown in squares.The dashed line represents the fit to 'IR-P' sources with sources below the dashed line having ǫ > 1.

Figure 4 .
Figure 4.The left panel shows the 870 µm image from ATLASGAL while the right panel shows the 850 µm image using the SMA in its compact configuration.The grey contours show emission at the −3σ level, while the black contours begin at 3σ and increase in factors of √ 2. The synthesized beam is shown in the bottom left corner of the SMA image.The field of view of the SMA map in the right panel is shown by the box in the left panel.The names of the fragments detected (see Table 3) are indicated in the SMA map.The plus signs indicate the locations of 24 µm point sources in the MIPSGAL survey.The images for the other sources are available in Appendix B.
Figure 4.The left panel shows the 870 µm image from ATLASGAL while the right panel shows the 850 µm image using the SMA in its compact configuration.The grey contours show emission at the −3σ level, while the black contours begin at 3σ and increase in factors of √ 2. The synthesized beam is shown in the bottom left corner of the SMA image.The field of view of the SMA map in the right panel is shown by the box in the left panel.The names of the fragments detected (see Table 3) are indicated in the SMA map.The plus signs indicate the locations of 24 µm point sources in the MIPSGAL survey.The images for the other sources are available in Appendix B.

Figure 5 .
Figure 5.The left panel shows the dendrogram while the right panel shows the hierarchical structures in the source AG49.25−0.41.The dashed line in the left panel shows the 1σ noise.The contour levels in the right panel begin at the 1σ level that defines the outer boundary of each root, with subsequent contours indicating the levels of nodes at which the root branches into two sub-structures.The synthesized beam is shown in the bottom left corner of the right panel.

Figure 7 .
Figure 7.The left and right panels show AG36.90−0.41 and AG41.05−0.25 respectively as seen by the SMA at 0.6 ′′ resolution.The images have been made by combining the data from the compact and extended configurations of the SMA.The grey contours show emission at the −3σ level.The thin black contour is at 1σ, while the thick black contours start at 3σ and increase thereafter by factors of √ 2. The synthesized beam is shown in the bottom left.

Figure 8 .
Figure 8.The ratio of the nearest neighbor separations between fragments and the thermal Jeans length (left panel) and turbulent Jeans length (right panel).

Figure 9 .
Figure 9.The left and right panels show the continuum measured by ATLASGAL towards AG41.07−0.12 and AG47.05+0.25 with white contours depicting emission detected by the SMA at 2 ′′ resolution.The position of the ATLASGAL clump is marked by the cyan cross.

Figure 10 .
Figure 10.The core formation efficiency (CFE) as a function of (a) the evolutionary state parameter (left panel) and (b) mean H2 number density (right panel).

Figure 11 .
Figure 11.The left panel shows the maximum mass of a fragment as a function of bolometric luminosity of the clump.The right panel shows the number of fragments as a function of the maximum fragment mass in a clump.

BFigure 12 .
Figure 12.The left panel shows the 870 µm image from ATLASGAL while the right panel shows the 850 µm image using the SMA in its compact configuration.The grey contours show emission at the −3σ level (see Table 2), while the black contours begin at 3σ and increase in factors of √ 2. The synthesized beam is shown in the bottom left corner of each SMA image.The field of view of the SMA map in the right panel is shown by the box in the left panel.The names of the fragments detected in each source (see Table 3) are indicated in the SMA maps.The plus signs indicate the locations of 24 µm point sources in the MIPSGAL survey.
Figure 12.The left panel shows the 870 µm image from ATLASGAL while the right panel shows the 850 µm image using the SMA in its compact configuration.The grey contours show emission at the −3σ level (see Table 2), while the black contours begin at 3σ and increase in factors of √ 2. The synthesized beam is shown in the bottom left corner of each SMA image.The field of view of the SMA map in the right panel is shown by the box in the left panel.The names of the fragments detected in each source (see Table 3) are indicated in the SMA maps.The plus signs indicate the locations of 24 µm point sources in the MIPSGAL survey.

Table 2 .
Details of SMA observations

Table 3 .
Fragments detected at 850 µm using the compact configuration of SMA

Table 4 .
Fragmentation properties of the ATLASGAL clumps