The ALMA Survey of 70 μm Dark High-mass Clumps in Early Stages (ASHES). XI. Statistical Study of Early Fragmentation

Fragmentation during the early stages of high-mass star formation is crucial for understanding the formation of high-mass clusters. We investigated fragmentation within 39 high-mass star-forming clumps as part of the Atacama Large Millimeter/submillimeter Array Survey of 70 μm Dark High-mass Clumps in Early Stages (ASHES) survey. Considering projection effects, we have estimated core separations for 839 cores identified from the continuum emission and found mean values between 0.08 and 0.32 pc within each clump. We find compatibility of the observed core separations and masses with the thermal Jeans length and mass, respectively. We also present subclump structures revealed by the 7 m array continuum emission. Comparison of the Jeans parameters using clump and subclump densities with the separation and masses of gravitationally bound cores suggests that they can be explained by clump fragmentation, implying the simultaneous formation of subclumps and cores within rather than a step-by-step hierarchical fragmentation. The number of cores in each clump positively correlates with the clump surface density and the number expected from the thermal Jeans fragmentation. We also find that the higher the fraction of protostellar cores, the larger the dynamic range of the core mass, implying that the cores are growing in mass as the clump evolves. The ASHES sample exhibits various fragmentation patterns: aligned, scattered, clustered, and subclustered. Using the Q -parameter, which can help distinguish between centrally condensed and subclustered spatial core distributions, we finally find that in the early evolutionary stages of high-mass star formation, cores tend to follow a subclustered distribution.


INTRODUCTION
High-mass stars (> 8 M ⊙ ) are mostly formed in clusters (Lada & Lada 2003).Investigating the fragmentation process of molecular clouds provides insights into the formation process of these stellar cluster members through molecular cloud contraction and has implications to the stellar initial mass function (IMF;Motte et al. 2022).In particular, Infrared dark clouds (IRDCs) are thought to be the best place to study the initial phase of high-mass stellar clusters prior to being significantly disturbed by feedback processes (e.g., outflows, stellar winds, and radiation; Rathborne et al. 2006;Chambers et al. 2009;Sanhueza et al. 2010Sanhueza et al. , 2012;;Tan et al. 2013;Rosen et al. 2020).Among IRDCs, those that are 70 µm dark are likely the most pristine, with no evidence of star formation in the IR (e.g., Sanhueza et al. 2013;Tan et al. 2013;Guzmán et al. 2015;Sanhueza et al. 2017;Contreras et al. 2018).Stars form in gas condensations where gravity dominates over any supporting mechanism such as turbulence, pressure, and magnetic fields.To understand what dominates core formation, separations and masses of cores have been compared with those expected from Jeans instabilities (Jeans 1902).Considering a uniform, infinite, isothermal medium at rest with a density ρ and pressure P 0 = ρc 2 s , the mean separation of fragments expected from the gravitational collapse, called the thermal Jeans length, is defined as where c s = (k B T /µm H ) 1/2 is the isothermal sound speed, G is the gravitational constant, k B is the Boltzmann constant, and m H is the mass of the hydrogen atom.H 2 and He govern the thermal velocity dispersion, and the mean molecular weight per free particle µ can be set to be 2.37 (Kauffmann et al. 2008), which is calculated from the cosmic abundance ratios.The mass of a sphere associated with the Jeans length is called Jeans mass and is expressed as If we include the impact of non-thermal motions as well as the thermal velocity dispersion, the Jeans length and mass are called turbulent Jeans length (λ tu J ) and mass (M tu J ), respectively.
The recently observed core separations in IRDCs or high-mass star-forming regions are mostly comparable to the thermal Jeans length of clumps rather than to the turbulent Jeans length (e.g., Beuther et al. 2015Beuther et al. , 2018;;Palau et al. 2018;Liu et al. 2019;Sanhueza et al. 2019;Lu et al. 2020;Beuther et al. 2021, and Ishihara et al. submitted).For core masses, observations of IRDCs using ALMA revealed that a large portion of cores have low-to intermediate-mass (≲ 30M ⊙ ), which also preferentially exhibit thermal Jeans fragmentation (e.g., Sanhueza et al. 2019), although less sensitive observations at lower resolution with the ALMA 7 m-array or the Submillimeter Array (SMA) have found most massive cores comparable to the turbulent Jeans mass (> a few M ⊙ ; Zhang et al. 2009Zhang et al. , 2015)).Only a few ALMA studies of more evolved high-mass star-forming regions affected by feedback from massive stars show less effective fragmentation, favoring turbulent fragmentation (Rebolledo et al. 2020;Jiao et al. 2023).Some recent studies, on the other hand, propose hierarchical fragmentation (i.e., from clump to sub-clump and sub-clump to cores; Palau et al. 2018;Pokhrel et al. 2018;Svoboda et al. 2019;Rosen et al. 2020;Zhang et al. 2021) based on two different spatial resolution observations or double peak of core separation distribution.It should be noted that most of these studies are case studies (limited sample) or singlepointing observations (limited spatial area).The general trend of core separations and core masses remains unclear, and a statistical study offers the opportunity to conclusively answer the following questions: Are average core properties explained by thermal Jeans fragmentation of clumps?Is clump fragmentation hierarchical?
We use 839 cores identified from continuum emission of thirty-nine 70 µm-dark IRDC clumps using the dendrogram algorithm in the ASHES survey (Sanhueza et al. 2019;Morii et al. 2023).This is the largest sample of cores embedded in IRDCs and is thought to be the best sample for a statistical study of cores in the very early phase of evolution to date.Morii et al. (2023) revealed that the majority of the ASHES clumps only host low-to intermediate-mass cores, implying the need for core growth, and that core mass segregation does not clearly appear at such very early stages, although cores are likely segregated in terms of their density.This paper is constructed as follows: Section 2 describes the observation setups and data reduction process.Section 3 summarizes the analysis to measure core separations and estimate core masses.We compare observed core properties (separations and masses) with thermal and turbulent Jeans parameters and with Jeans parameters estimated from clump and sub-clump densities, and discuss the hierarchical fragmentation, and diversity found within our sample in Section 4. Our conclusions are listed in Section 5.

OBSERVATIONS AND DATA REDUCTION
We have used observations from the ASHES survey, which was carried out with ALMA in Band 6 (∼224 GHz; ∼1.34 mm) through three cycles: Cycle 3 (2015.1.01539.S, PI: P. Sanhueza), Cycle 5 (2017.1.00716.S, PI: P. Sanhueza), and Cycle 6 (2018.1.00192.S, PI: P. Sanhueza).The data was taken with the main 12 m array and the Atacama Compact Array (ACA), including both the 7 m array and total power (TP).Targets are thirty-nine 70 µm-dark IRDC clumps with the potential for high-mass star formation (see Table 1 in Morii et al. 2023).The whole IRDC clumps were covered by Nyquist-sampled ten-pointing and three-pointing mosaics with the 12 m array and the 7 m array, respectively.The mosaicked area corresponds to 0.97 arcmin 2 within 20% power point, equivalent to the effective field of view (FOV) of ∼1 ′ per target.The detailed observation setups such as on-source time and maximum recoverable scale for all sources are summarized in Table 2 of Morii et al. (2023).
Data reduction was carried out using the CASA software package versions 4. 5.3, 4.6, 4.7, and 5.4.0 for calibration and 5.4.0 and 5.6.0 for imaging (CASA Team et al. 2022).Continuum images were produced by averaging line-free channels.The effective bandwidth for continuum emission was ∼3.7 GHz.After subtracting continuum emission, we combined the 12 m array data with the 7 m array data using the CASA task concat, and then they were cleaned together.Additionally, we also produced continuum images only from the 7 m array data following the same procedure.In this work, we only used TP data of C 18 O (J = 2-1) line to estimate the velocity dispersion of the clumps because TP antennas do not provide continuum emission.Some of C 18 O (J = 2-1) data has already been presented in the ASHES pilot survey (Sabatini et al. 2022).We used TCLEAN with Brigg's robust weighting of 0.5 to the visibilities and an imaging option of MULTISCALE with scales of 0, 5, 15, and 25 times the pixel size, considering the extended components of IRDCs.The average 1σ root-mean-square (rms) noise level of the combined image is ∼0.094 mJy beam −1 , with a beam size of ∼1.′′ 2 (Morii et al. 2023).The root-mean-square (rms) noise level of the 7 m-array image is summarized in the second column of Table 4.
For molecular lines, we used the automatic cleaning algorithm for imaging data cubes, YCLEAN (Contreras 2018; Contreras et al. 2018) to CLEAN each spectral window with custom-made masks.We adopted a Briggs's robust weighting of 2.0 (natural weighting) to improve the signal-to-noise (S/N) ratio.The average synthesized beam size is ∼ 1. ′′ 4. The average rms noise levels are ∼0.03K for cubes with the velocity resolution of 1.3 km s −1 and ∼0.09K for cubes with the velocity resolution of 0.17 km s −1 .All images have 512 × 512 pixels with a pixel size of 0. ′′ 2, and all images shown in this paper are the ALMA 12 and 7 m combined, before the primary beam correction, while all measured fluxes are derived from the combined data corrected for the primary beam attenuation.

Core Sample
Using the dendrogram technique (Rosolowsky et al. 2008), Morii et al. (2023) identified cores using the 1.3 mm continuum images.They set a minimum value, F min , as 2.5σ, a minimum significance to separate them, δ, as 1.0σ, and the minimum number of pixels to be contained in the smallest individual structure, S min , as the half-pixel numbers of the beam area.Here, σ is a root-mean-square (rms) noise level of the continuum image.Since these parameters are optimistic values, they applied the additional constraint to the flux density to exclude suspicious structures.They have excluded cores with a flux density smaller than 3.5σ.Additionally, they eliminated cores at the edge of FOV.The total number of identified leaf structures is 839 from 39 clumps (see an example in Figure 1).For the analysis carried out in this work, we have adopted the core masses from Morii et al. (2023).

Classification
We classified cores into three evolutionary categories and another three categories in terms of their gravitational stability.As star-formation signatures, we followed Sanhueza et al. (2019) and adopted the detection of molecular outflows traced by CO J = 2-1 and SiO J = 5-4, and the detection of warm gas tracers (H 2 CO J Ka,Kc = 3 2,2 − 2 2,1 , H 2 CO J Ka,Kc = 3 2,1 − 2 2,0 , CH 3 OH J K = 4 2 − 3 1 , and HC 3 N J = 24 − 23), which all have the upper state energy higher than 45 K.We judged cores associated with outflows as outflow cores, and any detection of warm line emission without outflow detections as warm cores.In turn, cores without any detections of outflow and warm gas tracers are classified as prestellar core candidates.This is consistent with the classification of the pilot ASHES survey (Li et al. 2023).In the end, we have 514 prestellar core candidates, and 325 protostellar cores (222 warm cores and 103 outflow cores).
To determine the gravitational state of cores, we calculated the virial parameter (α), the ratio of the virial mass (M vir ) to core mass (M core ).The virial mass was calculated as follows: where σ 2 tot = σ 2 th + σ 2 nt is the velocity dispersion, R is the core radius, a = (1 − b/3)/(1 − 2b/5) is the correction factor for a power-law density profile (ρ ∝ R −b ), and β = (arcsin e)/e is the geometry factor (see Fall & Frenk 1983;Li et al. 2013, for detailed derivation).Here, we adopt a typical density profile index b = 1.6 (e.g., Beuther et al. 2002;Palau et al. 2014), and β = 1.2 (Li et al. 2023).The thermal velocity dispersion and the non-thermal velocity dispersion are given by σ 2 th = kT µmH and σ 2 nt = σ 2 obs − kT m obs , respectively.We assumed that the non-thermal component is independent of the molecular tracer and that σ obs is the observed velocity dispersion.m obs is the molecular weight of the molecule, here 30m H , corresponding to DCO + and N 2 D + .
The identified cores at early stages of evolution are in dense (>10 5 cm −3 ) and cold (<20 K) environments where CO depletion and a high-level deuteration occur (e.g., Caselli et al. 2002;Sabatini et al. 2020;Redaelli et al. 2021Redaelli et al. , 2022;;Sabatini et al. 2022;Sakai et al. 2022;Sabatini et al. 2023).The pilot survey revealed that N 2 D + J = 3 − 2 and DCO + J = 3 − 2 succeeded in tracing such quiescent dense gas, excellent lines to measure the gas velocity dispersion minimizing the contribution from more diffuse intra-clump gas and gas related to protostellar activity like outflows (Sakai et al. 2022).We fitted the core-averaged spectrum of the N 2 D + and DCO + taken from the 12 and 7 m arrays combined data with a single Gaussian.We judged if the emission is detected with signal-to-noise larger than 3. Following the pilot survey (Li et al. 2022), to increase the S/N of the weak line emission, the core-averaged spectrum is spectrally smoothed over two native channels, prior to Gaussian fittings, if it shows marginal ∼3σ confidence in the native spectral resolution.Table 1 summarize the fitting result such as the peak intensity, the velocity center, and the total velocity dispersion for each line emission.Following the discussion in Li et al. (2022), both lines trace the same physical location (dense gas), and if both lines are detected in a core, we generally use σ obs measured from DCO + , except for a few cases.We adopt the σ obs value if either line is detected.As a result, we obtained σ obs information of 492 cores among 839 cores.Following Li et al. (2023), we consider cores with α < 2 to be gravitationally bound cores, and cores with α > 2 are unbound, which might be transient objects if we ignore additional support of magnetic fields or external pressure.For cores without the detection of DCO + and N 2 D + , we classified them as non-detection cores.The derived α ranges from 0.06 to ∼10 with a median value of 1.2.Among 492 cores, 340 cores are classified as bound cores.The detailed results for each core's virial analysis will be presented in a following paper.

Core Separation
We used the minimum spanning tree (MST) method developed by Barrow et al. (1985) to measure core separation.MST connects structures (cores in this case), minimizing the sum of the length and determining a set of straight lines.The right panel of Figure 1 shows an example of the minimum core separation (edge length determined by MST, hereafter δ sep ) as black segments.This separation is the projected separation, and the real separation could be equal or longer1 .Taking the average of the projection effect, we divide the measured separation by a correction factor of π/4 (see Ishihara et al. submitted for the details).The mean δ sep after correction in each clump ranges from 0.08 pc to 0.32 pc, with a median of 0.15 pc.
In addition to the minimum core separation, which is the distance between a pair identified by the MST method, we also estimate the nearest separation for each core, which is the distance to the nearest fragments.The average of nearest core separation after correction in each clump ranges from 0.07 pc to 0.29 pc, with a median of 0.12 pc.

Sub-clump Identification
Inside clumps, cores are not always uniformly distributed but consist of sub-clusters.To study the fragmentation properties of such sub-clusters, we applied the dendrogram technique for the continuum images produced by data only taken by the 7 m-array.We used the same parameters as those used for the core identification but with the rms noise level measured in the 7 m-array continuum images (Table 4).For example, the middle panel of Figure 1 shows the 7 m-array continuum image for G024.524-00.139, and magenta contours represent the identified leaf structures.We define these leaves as sub-clumps.From all clumps, we identified from one to seven sub-clumps per clump, 135 sub-clumps in total.Table 2 gives peak position, beam-convolved size (major and minor full-width half maximum; FWHM), peak intensity, and flux density of sub-clumps identified by the dendrogram algorithm (the properties for all sub-clumps are summarized in a machine-readable table ).
Masses of sub-clumps are estimated from the flux density of 1.3 mm continuum emission assuming optically  thin conditions by where R = 100 is the gas-to-dust mass ratio, κ ν is the dust absorption coefficient, d is the heliocentric distance associated with each ASHES clump, and B ν is the Planck function for a dust temperature T dust .We adopt κ ν of 0.9 cm 2 g −1 from the dust coagulation model of the MRN (Mathis et al. 1977) distribution with thin ice mantles at a number density of 10 6 cm −3 computed by Ossenkopf & Henning (1994).We used the clump dust temperature listed in Morii et al. (2023).As the ASHES clumps are in their early stages, the free-free contamination especially at 1.3 mm is negligible.Here, we assume that the continuum emission comes only from dust emission.
The surface density, Σ, and the molecular volume density, n(H 2 ), were estimated assuming a uniform spherical density distribution as: half of the geometric mean of the FWHM (Table 2).
In this paper, we defined cores overlapped with the subclump as members of the sub-clump (black crosses in the middle panel of Figure 1).The number of cores in each sub-clump varies from zero to more than ten (the maximum is 23).The estimated physical parameters and the number of gravitationally bound cores are summarized in Table 3.

Jeans Length and Jeans Mass
We estimated the thermal Jeans length and mass of clumps and sub-slumps from Equations 1 and 2 with their density and dust temperature.The clump density is estimated from the flux density of the continuum emission obtained by the Atacama Pathfinder Experiment Telescope Large Area Survey of the Galaxy (ATLASGAL; Schuller et al. 2009), and summarized in Table 1 of Morii et al. (2023).One example of the continuum image is shown in the left panel of Figure 1.The derived λ th J,cl ranges from 0.06 pc to 0.22 pc, with a median of 0.12 pc, and M th J,cl ranges from 1.1 M ⊙ to 5.5 M ⊙ , with a median of 2.4 M ⊙ .Using the sub-clump density in Table 3, sub-clumps' thermal Jeans length (λ th J,sub−cl ) and mass (M th J,sub−cl ) can be derived.The density of sub-clumps is about one order of magnitude larger than clumps, and the estimated λ th J,sub−cl is a few times smaller than λ th J,cl , ranging from 0.018 pc to 0.12 pc, with a median of 0.05 pc.M th J,sub−cl ranges from 0.3 M ⊙ to 3.5 M ⊙ , with a median of 1.1 M ⊙ .
We also estimate the turbulent Jeans parameters of clumps.We applied a 1D Gaussian fitting toward the line profile of C 18 O J = 2-1 averaged within the clump, which was observed by TP.The obtained velocity dispersion (σ obs ) is summarized in  4.
Figure 2 shows the histogram of the ratios of the observed separation and mass divided by Jeans length and mass of clumps, respectively, in logarithmic scale.The top two panels show the case for the thermal Jeans fragmentation (δ sep /λ th J,cl and M core /M th J,cl ), while the bottom two panels are the cases of turbulent Jeans fragmentation (λ tu J,cl and M tu J,cl ).The δ sep /λ th J,cl distributions take a peak around 1-2.For mass, the distribution has a large variation across more than two orders of magnitudes and more than half of cores have lower mass than the thermal Jeans mass of clumps.The peak is around 0.2 but there is a secondary peak around unity.For the turbulent Jeans fragmentation, both ratios of δ sep /λ tu J,cl and M core /M th J,cl show their peaks at values much smaller than unity.

DISCUSSION
In this section, we aim to characterize the fragmentation at early stages of high-mass star formation found in the ASHES sample by comparing core separation and mass with thermal and turbulent Jeans parameters, and by comparing clump and sub-clump Jeans parameters to study hierarchical fragmentation.

Thermal versus Turbulent Jeans Fragmentation
What dominates the fragmentation is one key question in star formation that has been investigated in some previous studies (e.g., Zhang et al. 2009;Palau et al. 2013Palau et al. , 2015;;Sanhueza et al. 2017Sanhueza et al. , 2019;;Beuther et al. 2018Beuther et al. , 2024)).We have the largest sample in IRDCs thanks to     b Velocity dispersion (σ) was obtained by the fitting of the line profile of C 18 O (J = 2 − 1) averaged within the clump with a 1D Gaussian, which was observed by Total Power (TP).
c Number of bound cores within r = 0.45 pc (see Section 4.3).
d The fraction of protostellar core to all bound cores in each clump.
e The parameter to describe how centrally concentrated the core spatial distribution is.See 4.4, for more details.
the mosaicked high spatial resolution and high sensitivity ALMA observations.IRDCs are thought to be the best for answering this question in the very early phase of high-mass star and cluster formation.
To visualize that clearly, Figure 3 displays the massnearest separation relation, as in Wang et al. (2014, following this work, for this figure, we use the nearest separation between cores rather than the minimum separation defined in MST).In this figure, blue and green shaded regions show what is expected from thermal Jeans fragmentation and turbulent Jeans fragmentation, respectively.Our identified cores and sub-clumps are denoted as crosses and circles, respectively.The blue line shows thermal Jeans fragmentation with T = 15 K (mean temperature of the ASHES sample) varying the density from 10 2 cm −3 to 10 6 cm −3 .The blue-shaded region shows the same density range but with T from 10 K to 30 K. The green-shaded region shows turbulent Jeans fragmentation to the same density and temperature range but varying the velocity dispersion from 0.8 to 2.7 km s −1 .These ranges cover the ASHES sample properties (see Table 1 of Morii et al. 2023).The majority of cores (99.6 %) and all sub-clumps are plotted below the green area, and on average, both cores and sub-clumps are in the thermal Jeans fragmentation regime, while the variation of core masses is large, about two orders of magnitude.Especially gravitationally bound cores (red crosses) are located around blue-shaded regions, preferring thermal Jeans fragmentation.
This trend can also be seen in Figure 2. The δ sep /λ th J,cl distribution takes a peak around unity, and both the mean and median values of M core /M th J,cl are around unity (1.28 and 0.393, respectively).On the contrary, the Jeans length and mass in the turbulence-dominated case are both much larger than λ th J,cl and M th J,cl , and the ratios become much smaller than unity (see the bottom panels).Thus, based on the estimated masses and separations, we conclude that core formation found in the ASHES sample (and likely in IRDCs in general) is regulated by thermal Jeans fragmentation rather than turbulent fragmentation.
Some other ALMA observations with similar setups like the ASHES survey still suggest turbulencedominated fragmentation (e.g., Rebolledo et al. 2020;Xu et al. 2023b), but most are not in a quiet/early stage.They focus on more evolved, active high-mass star-forming regions where some feedback from the newborn stars can suppress the fragmentation by inducing additional turbulence and warming up the surrounding gas (e.g., Krumholz & McKee 2008).

Hierarchical fragmentation
We found that the overall properties of identified cores are consistent with thermal Jeans fragmentation as shown in Figure 3.However, as Figure 1 shows, cores are embedded in sub-clump structures inside clumps.We next address if cores form from the fragmentation of clumps or dense sub-structures (here we refer to them as sub-clumps) by using clumps' and sub-clumps' Jeans parameters.In this section, we use only gravitationally bound cores which will likely form stars. Core separations are recalculated only using bound core positions.
Figure 4 compares the fragmentation of clumps and sub-clumps with the observed core properties.The top two panels show the ratios of separation and masses normalized by thermal Jeans parameters estimated by using clump density, and the bottom two panels show the case of using sub-clump density.For the bottom panels, only cores inside sub-clumps are considered.The left panels present the ratios of separation of bound cores normalized by thermal Jeans length.Their distribution peaks are located around 1-2 and 2-3 for the top and bottom, respectively.The right panels display the histogram of the mass ratios and bound cores are highlighted in red.The peak is just around unity for clump fragmentation.Although the bottom panel shows no clear single peak, it has a broad peak around 1-5.These histograms also show that unbound cores and cores without the detection of dense gas tracers generally have masses smaller  than the Jeans mass of clumps.Overall, the observed core properties favor the fragmentation from clumps rather than sub-clumps.
Hierarchical fragmentation is expected in one of the theoretical scenarios for high-mass star formation, called Global Hierarchical Collapse (GHC) scenario (Vázquez-Semadeni et al. 2019).Such hierarchical fragmentation has been reported in some previous studies observing high-mass starless clump candidates (Svoboda et al. 2019;Zhang et al. 2021), IRDCs (Wang et al. 2011(Wang et al. , 2014)), and also OMC-1S (Palau et al. 2018).Our analysis also revealed sub-clumps, intermediate structures connecting clumps and cores, but the core properties cannot be explained well by sub-clump fragmentation.This implies that a step-by-step fragmentation (clump to sub-clump to core) is unlikely.
We also investigated the fragmentation from clump to sub-clump by comparing Jeans parameters with subclump separations and masses.We find no clear peak in the distributions due to the small number of statistics.However, thermal Jeans fragmentation may still be favored over turbulent fragmentation as shown in Figure 3.We note that the sub-clumps locate, in Figure 3, in the area at which the parental clumps should have densities of ∼10 3 cm −3 , inconsistent with observations.This could be explained if sub-clump structures form as a result of clump fragmentation into cores.Subclumps could be the ensemble of cores that evolve as cores evolve, not in a step-by-step hierarchical fragmentation, but rather in a simultaneous formation process.
We point out that this simultaneous formation of subclumps and cores, in which the core properties are not determined by sub-clumps, may only be valid at the very early stages of high-mass star formation traced in the ASHES survey.Later, in more active high-mass starforming regions, once subsequent fragmentation occurs, core properties may be explained differently.

Fragmentation Level and Clump Properties
Jeans fragmentation invokes the idea that the fragmentation level or the number of cores depends on clump density.We investigated the correlation of clump density or Jeans number (the ratio of clump mass to Jeans mass) with the number of bound cores, n(Bound core).To calculate n(Bound core), we count the number of bound cores within the same physical area for all clumps and impose a mass threshold to reduce the effect of having different distances and sensitivities for each clump.We have limited the sample for this discussion, excluding clumps that are located too close (< 3.5 kpc) and too far (> 5.5 kpc), and two more with the worst mass sensitivity (> 0.45M ⊙ ).As a result, the 30 clumps remaining are located between 3.5 and 5.5 kpc and have a mass sensitivity between 0.086 and 0.41 M ⊙ .We count cores within a circle with a radius of 0.45 pc centered on the mean positions of cores.The circle size almost corresponds to the FoV of the closest clump.In addition, we impose for all clumps a mass threshold of 0.41 M ⊙ , which corresponds to the worst mass sensitivity among the 30 clumps.The measured n(Bound core) are listed in Table 4).
The top-left two panels of Figure 5 show moderate to strong correlations of n(Bound core) with surface density Σ cl and mean clump number density n(H 2 ) cl .It indicates denser clumps produce more number of cores (higher fragmentation level).We found stronger correlation between n(Bound core) and Σ cl with a Spearman's rank correlation coefficient of r s = 0.74, while that is r s = 0.69 for n cl .Both p-values are much smaller than 0.01.The less scattered plot of n(Bound core)-Σ cl implies that the clump surface density is the best indicator of the fragmentation level: the higher the clump surface density, the more likely it is to have a larger number of cores.To confirm that these correlations do not result from the co-dependence on the distance, we re-calculated clump mass (M cl,0.45pc ) and surface density (Σ cl,0.45pc ) using flux within the same physical area (r=0.45 pc).The right panel displays the strong correlation between Σ cl,0.45pc and n(Bound core).We also compared n(Bound core) with the Jeans numbers (M cl,0.45 pc /M th J,0.45 pc ), which is proportional to Σ 3 cl,0.45pc /n(H 2 ) 3/2 cl,0.45 pc .It is confirmed that the measured n(Bound core) has a strong correlation with the number of cores expected from thermal Jeans fragmentation ( r s = 0.70 and p-value<< 0.01).
These correlations are still found in the case of subclump to core fragmentation, but in this case, the coefficients are relatively weaker.The bottom three panels also imply the strong or moderate correlation between the number of bound cores in each sub-clump and subclump surface density, volume density, and Jeans number, although relatively weaker than the clump case.It suggests that sub-clumps are also involved with core formation.
To summarize, we revealed that a higher fragmentation level (or higher core number density) can be expected from a region with a higher surface density and it is consistent with Jeans fragmentation.The tight correlation between the fragmentation level and the clump/cloud surface density has been observed in more evolved star-forming regions as well, indicating that this correlation begins early on in IRDCs and prevails during the evolution of high-mass star-forming regions (e.g., Palau et al. 2014;Sokol et al. 2019).

Fragmentation Diversity
The ASHES survey also reveals a diversity in both the range of core masses per clump (mass dynamic range) and fragmentation patterns.Morii et al. (2023) reported that most clumps host low-to intermediate-mass cores.However, the dynamic range in core masses varies from clump to clump.For example, the left panel in Figure 6 shows some (relatively) massive cores surrounded by some low-mass cores, and the right panel shows a cluster of low-mass cores with a small dynamic range in mass.We find that such a mass dynamic range correlates with the fraction of protostellar cores (cores with outflow or warm line detection).Figure 7 shows the fraction of bound protostellar cores to the total number of bound cores as a function of the maximum and minimum mass difference (M max −M min ), the ratio (M max /M min ), the standard deviation of core masses (σM ), and interquartile range.The fraction of protostellar cores is summarized in Table 4.The x-axis of the top two panels is derived from the maximum and minimum core masses and can be significantly affected by a single peculiar object if exists.The interquartile range, the difference of 25 percentile and 75 percentile, is on the contrary, less affected by the most massive object.The Spearman's rank correlation coefficients are r s ∼ 0.5 and p-values are less than 0.01 except for the case of M max /M min .One outlier in these plots is G340.232-00.146,shown as an unfilled marker.As discussed in Sanhueza et al. (2019), in this peculiar clump, the most massive core is rather large with a radius of ∼10 4 au and more fragmented structures are expected in higher angular resolution observations from visual in-  spection of the continuum image.The correlation coefficients displayed in Figure 7 become higher if we exclude it as shown in the values in parenthesis (bottom, right of each panel).These plots imply that if we take the fraction of protostellar cores as an indicator of cluster evolution, the dynamic range of mass increases with cluster evolution.This correlation can be interpreted as clump-fed accretion onto cores.Clumps would initially fragment into cores of similar mass (small dynamic mass range), and as time goes on, some cores grow in mass more than others, resulting in a larger mass range difference (large dynamic mass range).
Among 39 ASHES targets, some clumps show aligned fragmentation as the left panel of Figure 6, some show concentrated, and some show spread core distributions with several sub-clumps (e.g., right panel of Figure 6).We calculated the Q-parameter to investigate the cluster members' distribution.The Q-parameter was defined by Cartwright & Whitworth (2004) as The term m is the normalized mean edge length of MST and is defined as where N c is the number of cores in the region, Σ Nc i=1 L i is the total length of all the lines MST connected, hereafter 'edges', and A is the area of the cluster and estimated by A = πR 2 cluster .Here, the radius of the cluster R cluster is defined as the distance from the mean position of cores to the farthest core position.The second term in the first equation of Equation 6is the factor to normalize the mean edge length of cores (the first term is defined as l MST ) having different areas (A) and/or different numbers of cores (N c ).The term s is the ratio of the mean core separation to the cluster radius (R cluster ).Here, the core separation is different from the minimum core separation (δ sep ), and it is the core separation within the region, not only considering the minimum separation.Now both m and s are independent of the number of cores in the cluster-forming clump.
For clusters with a smooth radial density gradient (n ∝ r −α ), Q increases from ∼0.8 to 1.5 as the degree of concentration increases from α = 0 to 2.9, and for sub-clustering clusters Q becomes smaller than 0.8 (Cartwright & Whitworth 2004).The parameters estimated for the ASHES sample range from 0.6 to 0.9, and their average is 0.76.Some clumps indicate a uniform distribution of cores (Q ∼0.8), but most (70%) prefer sub-clustering (Q < 0.8).Clumps with aligned fragmentation have a Q-parameter of 0.7-0.8.However, we note that, as can be inferred from the definition of the Q-parameter, this parameter is unable to identify aligned fragmentation.Clumps with several sub-clumps or showing spread fragmentation have Q < 0.8.
The origin of such variation in the fragmentation pattern is not yet clear from the current data (e.g., clump mass, density, virial parameter, protostellar, or core fraction), and further information on the magnetic field or clump-scale properties such as large-scale gas dynamics seems to be necessary.For example, Tang et al. (2019) suggests that the balance among the magnetic field, turbulence, and gravity determines the core fragmentation pattern such as no fragmentation, aligned fragmentation, and clustered fragmentation.It should be noted that we find no clear correlation between the dynamic range in mass and the fragmentation pattern; both aligned fragmentation and spread fragmentation show large and small mass dispersion.

Early Fragmentation picture
We have revealed that the observed mean core separation and masses are comparable to thermal Jeans lengths and masses, respectively, and much smaller than turbulent Jeans parameters.It implies that turbulence is not a dominant source characterizing core formation.This is consistent with the study by Traficante et al. (2020), indicating that gravity dominates over turbulence once the regions become dense (e.g., Σ > 0.1 g cm −2 ).ASHES targets are all dense with a surface density larger than 0.1 g cm −2 .Compared with other ALMA studies in high-mass star-forming regions, in which core properties are better explained by turbulent Jeans fragmentation (e.g., Rebolledo et al. 2020;Xu et al. 2023b), the ASHES sample contains 70 µm-dark, cold regions not affected by feedback mechanism from massive stars.Thus, our finding implies that the initial fragmentation in massive clumps, prior to the changes due to gravitational accretion and some feedback effects, is described by thermal Jeans fragmentation.
However, our sample still contains super-Jeans cores with a mass more than 10 times M th J,cl .Those cores may have grown in mass by acquiring gas from the surrounding environment.Infall rates in the range 10 −4 -10 −3 have been measured in two ASHES targets (Contreras et al. 2018;Redaelli et al. 2022), allowing the cores to quickly grow in mass in a free-fall time.Alternatively, the magnetic field may play a role in suppressing fragmentation as suggested by theoretical studies (Hennebelle & Teyssier 2008;Commerçon et al. 2011).The theoretical prediction that the magnetic field can inhibit fragmentation in high-mass star-forming regions has garnered support from observational studies in more evolved high-mass star-forming regions.Dust polarization emission from infrared-bright sources has highlighted the significant role magnetic fields play in exerting pressure on the medium, from clump to core scale, effectively curbing fragmentation during gravitational collapse (Zhang et al. 2014;Hull & Zhang 2019).Supporting evidence for fragmentation suppression due to magnetic fields has also been presented by Frau et al. (2014).In addition, Das et al. (2021) showcased the magnetic field's impact on reducing fragment numbers, while Palau et al. (2021) reported a tentative correlation between fragment quantity and the mass-to-flux ratio among massive dense cores, as suggested by theoretical and numerical works.Future observations of dust polarisation toward such sub-clumps or massive cores would verify this effect on suppressing fragmentation.Studying how significantly the magnetic field contributes to the fragmentation is also important to the understanding of the diversity of fragmentation patterns seen in Figure 6.Although at relatively smaller scales, these points are one of the goals of the Mag-MaR (Magnetic Fields in Massive Star-forming Regions; Fernández-López et al. 2021;Cortés et al. 2021;Sanhueza et al. 2021) survey once the whole survey sample is analyzed.
We found sub-structures inside clumps using the 7 marray data, which is consistent with the measured Q of 0.6-0.8,implying that the initial core distribution is not yet so concentrated but rather sub-clustered.These sub-structures are located in the thermal Jeans fragmentation-dominated regime in Figure 3, as well as cores.We therefore investigated whether cores are directly formed from such sub-clumps rather than from clumps.We find no evidence that sub-clump fragmentation is the preferred mechanism to explain the observed core properties.The comparison of the number of cores or the degree of fragmentation with the properties of clumps or sub-clumps also suggests a stronger link between clumps and cores.A possible picture is that sub-clumps and cores are simultaneously formed from clumps, and core properties are determined from clump properties.
It should be noted that such sub-clumps are likely to contribute to density segregation since denser cores are generally embedded in sub-clumps.If they can grow by more effective gas feeding, this may later lead to mass segregation as discussed in Morii et al. (2023).Gas feeding or gravitational collapse of clumps and subclumps would increase the mass dynamic range as seen in Figure 7. Xu et al. (2023a) suggests that the gravitational concentration or gas accretion towards the center of mass would result in the appearance of mass segregation and in the increase of the Q-parameter.Comparing our results with more evolved clusters would confirm this hypothesis.

CONCLUSIONS
We have studied the fragmentation properties in 39 clumps as a part of the ALMA Survey of 70 µm dark High-mass clumps in Early Stages (ASHES), which aims to characterize the very early phase of high-mass star formation.Using the 839 cores identified in the continuum images, we compared their masses and separations with Jeans parameters.We have obtained the following conclusions: 1.The mean core separation measured by the MST method ranges from 0.08 pc to 0.32 pc in each region, and core masses range from 0.05 M ⊙ to 81 M ⊙ .The core mass and core separation are explained by thermal Jeans fragmentation ruling out turbulent Jeans fragmentation at the very early stages of high-mass star formation.
2. Comparing the Jeans parameters of clumps and sub-clumps with the observed core properties, core properties, especially for bound cores, are likely determined from clumps.We interpret this as a simultaneous formation of sub-clumps and cores within clumps.
3. The fragmentation level or the number of cores within each clump shows a strong correlation with the Jeans number, the ratio of clump mass to Jeans mass, implying that early core formation can be described with thermal Jeans fragmentation.It also has a strong correlation with clump surface density.
4. Furthermore, our sample shows the diversity of fragmentation in terms of mass dynamic range and spatial distribution.The correlation between the protostellar core fraction and the mass dynamic range is likely a sign of the clump-fed accretion scenarios.We have revealed aligned, spread, clustered, and sub-clustered fragmentation patterns, and the measured Q-parameter also implies that the early fragmentation seen in ASHES fields is not centrally concentrated.Facility: ALMA 26 Software: CASA (v4.5.3, 4.6, 4.7, 5.4, 5.6;CASA Team et al. 2022)

Figure 1 .
Figure 1.Continuum images for G024.524-00.139obtained by (left) a single-dish telescope (Schuller et al. 2009), (middle) ALMA 7 m array, and (right) ALMA 12 m array and 7 m array.The beam sizes of each image are shown as a black circle or ellipse at the bottom left.(left) Gray-solid contours represent 3 × 2 n σ (n=1, 2, 3, ...), where σ = 71 mJy beam −1 is the rms noise level.The black contour represents the FoV of ALMA observations.(middle) Gray-solid contours represent 3 × 2 n σ (n=1, 2, 3, ...) with σ = 0.53 mJy beam −1 .Magenta thick contour represents leaf structures identified by the dendrogram algorithm.Black and red crosses show core peak positions included in sub-clumps and not-included ones, respectively.The black ellipse in the bottom left corner represents the synthesized beam size.The black line indicates the spatial scale in the bottom right corner.(right) Black segments show the outcome from the minimum spanning tree, which corresponds to the set of straight lines that connect cores in a way that minimizes the sum of the lengths.Gray contour levels are the same with the middle panel but with σ = 0.095 mJy beam −1 .Magenta thick contour represents leaf structures identified by the dendrogram algorithm.

Figure 2 .Figure 3 .
Figure 2. Core separation (δsep) and core masses normalized with Jeans length and Jeans mass of clumps, respectively.Blue and green bins represent the case of thermal or turbulent Jeans fragmentation, respectively.The thick solid lines represent kernel density distribution.The vertical lines correspond to the ratio of unity.

Figure 4 .
Figure 4. Core separation (δsep) and core masses normalized with thermal Jeans length and masses, respectively.The top three panels show ratios normalized by thermal Jeans lengths and masses of clumps, and the bottom three show the case of sub-clump fragmentation.Only cores inside sub-clumps are used for this analysis.The thick solid lines represent kernel density distribution.The vertical lines correspond to the ratio of unity.(Left) Core separations of bound cores normalized by Jeans length.(Right) Core masses normalized by Jeans mass.Gravitationally bound and unbound are colored red and blue, respectively (cores without detections of dense gas tracers are in gray).

Figure 5 .
Figure 5. (Top) The number of bound cores above 0.41 M⊙ located around the clump peak (r < 0.45 pc), n(Bound core), as a function of clump surface density Σ cl , clump mean number density n(H2) cl , Jeans number, the ratio of M cl,0.45pc to thermal Jeans mass (M th J,0.45pc ), and clump surface density with a limited area Σ cl,0.45pc for selected 30 clumps.Spearman's rank correlation coefficients are denoted inside each panel.All p-values are much smaller than 0.01.(Bottom) The number of bound cores above 0.41 M⊙ in sub-clumps, N (Bound core), as a function of sub-clump surface density Σ sub−cl , mean number density n(H2) sub−cl , and Jeans number (M sub−cl /M th J,sub−cl ).

Figure 6 .Figure 7 .
Figure 6.ALMA 1.3 mm continuum image of (left) G024.010+00.489 and (right) G028.273-00.167.The circle size represents the core mass, and the position is centered at the continuum peak of each core.The three different colors have meanings the same as Figure 3.

Table 2 .
Sub-clump Properties Obtained from Dendrogram

Table 4 .
The non-thermal velocity dispersion is given by σ 2 nt = σ 2 obs − kT m obs with m obs = 28m H .By replacing c s with c 2 s + σ 2 nt in Equa-tion 1, the turbulent Jeans length (λ tu J,cl ) can be estimated.It is about five times larger than the thermal Jeans length ranging from 0.31 pc to 2.35 pc, and the turbulent Jeans masses (M tu J,cl ) are estimated in a range of 10 2 to 7×10 3 M ⊙ .All estimated values for each clump are summarized in Table

Table 3 .
Sub-clump Physical Parameters Note-The radius is calculated from the geometric mean of the FWHM divided by 2. (This table is available in its entirety in machine-readable form.)

Table 4 .
Clump Information and Jeans Parameters K.M. is financially supported by Grants-in-Aid for the Japan Society for the Promotion of Science (JSPS) Fellows (KAKENHI Number 22J21529) and supported by FoPM, WINGS Program, the University of Tokyo.K.M. is also supported by JSPS Overseas Challenge Program for Young Researchers (202280210).PS was partially supported by a Grant-in-Aid for Scientific Research (KAKENHI Number JP22H01271 and JP23H01221) of JSPS.GS acknowledges the projects PRIN-MUR 2020 MUR BEYOND-2p ("Astrochemistry beyond the second period elements", Prot.2020AFB3FX) and INAF-Minigrant 2023 TRIESTE ("TRacing the chemIcal hEritage of our originS: from proTostars to planEts"; PI: G. Sabatini).Data analysis was in part carried out on the Multi-wavelength Data Analysis System operated by the Astronomy Data Center (ADC), National Astronomical Observatory of Japan.This paper uses the following ALMA data: ADS/JAO.ALMA#2015.1.01539.S, ADS/JAO.ALMA#2017.1.00716.S, and ADS/JAO.ALMA#2018.1.00192.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), M OST 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.