Fragmentation of massive dense cores down to ~1000 AU: Relation between fragmentation and density structure

In order to shed light on the main physical processes controlling fragmentation of massive dense cores, we present a uniform study of the density structure of 19 massive dense cores, selected to be at similar evolutionary stages, for which their relative fragmentation level was assessed in a previous work. We inferred the density structure of the 19 cores through a simultaneous fit of the radial intensity profiles at 450 and 850 micron (or 1.2 mm in two cases) and the Spectral Energy Distribution, assuming spherical symmetry and that the density and temperature of the cores decrease with radius following power-laws. We find a weak (inverse) trend of fragmentation level and density power-law index, with steeper density profiles tending to show lower fragmentation, and vice versa. In addition, we find a trend of fragmentation increasing with density within a given radius, which arises from a combination of flat density profile and high central density and is consistent with Jeans fragmentation. We considered the effects of rotational-to-gravitational energy ratio, non-thermal velocity dispersion, and turbulence mode on the density structure of the cores, and found that compressive turbulence seems to yield higher central densities. Finally, a possible explanation for the origin of cores with concentrated density profiles, which are the cores showing no fragmentation, could be related with a strong magnetic field, consistent with the outcome of radiation magnetohydrodynamic simulations.


INTRODUCTION
One unavoidable ingredient in the process of formation of intermediate/high-mass stars is that they are usually found associated with clusters of lower-mass stars (e. g., Lada & Lada 2003;Rivilla et al. 2013). However, the process yielding the formation of these clusters is not clear, as it remains ambiguous how a cloud core fragments to finally form a cluster. The observed star formation efficiency is by far much smaller than the predicted if gravity dominated the dynamics of the collapse (see e. g., Kruijssen 2013 for a review), and the number of subcondensations typically embedded in massive dense cores is again smaller by one order of magnitude than the predicted from a pure Jeans fragmentation analysis (e. g., Hennemann et al. 2009;Zhang et al. 2009Zhang et al. , 2011Bontemps et al. 2010;Longmore et al. 2011;Pillai et al. 2011;Wang et al. 2011;Lee et al. 2013). Thus, a broad observational evidence indicates that additional physical ingredients competing with gravity must play an important role in the fragmentation process, such as angular momentum, stellar feedback, magnetic fields and turbulence.
Numerical simulations are becoming crucial to establish the relative importance of each of these ingredients in the fragmentation process. Simulations of fragmentation of rotating cores with different equations of state show that cores with some angular momentum easily fragment into multiple objects (e. g., Boss & Bodenheimer 1979; and thus need of additional mechanisms to decrease the fragmentation down to the observed levels (see Goodwin et al. 2007 and references therein). A crucial ingredient in this sense is the presence of a magnetic field. Studies on the effects of magnetic fields in the evolution of a collapsing and rotating core indicate that a strong magnetic field can suppress fragmentation as it constitutes an additional form of support and also via the so-called 'magnetic braking' (e. g., Hosking & Whitworth 2004;Mellon & Li 2008;Hennebelle & Teyssier 2008), which is very effective in removing a significant part of the initial angular momentum through torsional Alfvén waves. Similarly, radiation magnetohydrodynamic simulations of turbulent massive dense cores show that the combined effects of radiative feedback (e. g., Krumholz et al. 2007; Peters et al. 2011) and magnetic field can sum to efficiently suppress fragmentation (e. g., Hennebelle et al. 2011Myers et al. 2013), yielding a better match with the observations (e. g., Palau et al. 2013). In addition, submillimeter observations show that the magnetic field is important in the dynamical evolution and fragmentation of massive dense cores (e. g., Girart et al. 2009Girart et al. , 2013Hezareh et al. 2013). Another ingredient which is probably playing a role in the fragmentation process is supersonic turbulence. Theoretical studies show that turbulence generates density structures due to supersonic shocks that compress and fragment the gas efficiently, especially for large Mach numbers (e. g., Padoan et al. 2001, Schmeja & Klessen 2004Vázquez-Semadeni et al. 2007). More recently it has been shown that the turbulence mode, solenoidal or compressive (e. g., Federrath et al. 2008), affects differently the cloud structure and fragmentation. Solenoidal turbulence corresponds to a null divergence of the driving force (divergence-free) and arises from galactic rotation and magneto-rotational instabilities, and seems to be at work in the Rosette and Polaris Flare clouds (e. g., Hily-Blant et al. 2008;Federrath et al. 2010). Compressive turbulence corresponds to a null curl of the driving force (curl-free) and can be produced by expanding supernova shells or HII regions (e. g., Federrath et al. 2010), as seen in the Pipe Nebula and the Orion B cloud (Peretto et al. 2012;Schneider et al. 2013). Numerical simulations comparing both types of turbulence indicate that compressive turbulence promotes fragmentation and generates high density contrast structures, while solenoidal turbulence slows fragmentation down and keeps a smooth density structure in the cloud (e. g., Federrath et al. 2010;Girichidis et al. 2011).
In addition to these physical agents opposing to gravity, the outcome of numerical simulations also depends on the geometry and physical properties of the collapsing core, in particular the density profile. Hydrodynamical simulations show that highly concentrated cores (e. g., n H 2 ∝ r −2 ) lead to the formation of one single massive star, while flat density profiles (n H 2 ∝ r −1 or flatter) produce an important number of low-mass stars (e. g., Myhill & Kaula 1992;Burkert et al. 1997;Girichidis et al. 2011). This trend appears to be seen in high resolution observations of IRDC G28.34+0.06, where three massive dense cores with different density profiles show varying degrees of fragmentation . Finally, the substructure of molecular gas within a dense core is decreasing with time as the core turns gas into stars, as shown for example by Smith et al. (2009).
From an observational point of view, most of the studies focusing on fragmentation of massive dense cores deal with single regions and/or do not reach high enough spatial scales to study individual fragments in the core (∼ 1000-2000 AU; e. g., Longmore et al. 2011;Pillai et al. 2011;Miettinen 2012;Beuther et al. 2013;Kainulainen et al. 2013;Takahashi et al. 2013). Recently, Bontemps et al. (2010) study a sample of six massive dense cores in Cygnus X and conducted a consistent analysis for all the regions, finding cores with a variety of fragmentation levels. A subsequent study by Palau et al. (2013) focused on a list of 18 massive dense cores (including Bontemps et al. 2010 sources), and observed with similar spatial resolution (∼ 1000 AU), sensitivity (down to ∼ 0.3 M ⊙ ), and field of view (∼ 0.1 pc of diameter), again reveals different fragmentation levels, with about ∼ 30% of the regions showing no fragmentation at all. Palau et al. (2013) study possible correlations between the fragmentation level and properties of the cores such as mass, average density, evolutionary indicators, rotational-to-gravitational energy and turbulent ve-locity dispersion (although these last two properties should be studied in a larger sample), and find no clear trend of fragmentation level with any of these core properties. In that work, two main properties of the cores remained to be studied in relation to fragmentation, mainly the magnetic field and the density structure. While the magnetic field cannot currently be assessed in this sample in a uniform way, the density structure has been studied in roughly half of the sample (e. g., Beuther et al. 2002a;Mueller et al. 2002;Crimier et al. 2010;van der Tak et al. 2013), but the methodology used is different for each work, casting some doubt on the validity of any comparison. Here we present a uniform study of the density structure of each core in the sample of Palau et al. (2013), aiming at accurately assessing its impact on the fragmentation level.
In Section 2 we present the sample of cores and the data used to extract the submillimeter intensity profiles and the Spectral Energy Distributions (SEDs), and we also present the envelope model used to fit these data simultaneously for each core; in Section 3, we present the results of our fits; in Section 4, we discuss the possible trend between fragmentation level and density structure. Finally, in Section 5 we summarize our main conclusions.

Data
We aim at studying the density structure of the 18 massive dense cores presented in Palau et al. (2013) plus DR21-OH, a massive dense core for which fragmentation has been recently studied at similar mass sensitivities and spatial resolution by Girart et al. (2013). Among the total number of 19 massive dense cores, 16 of them are identified as submillimeter sources in the catalog of Di Francesco et al. (2008) Bontemps et al. 2010), we used the submillimeter data of Jenness et al. (1995), who observed also the 450 and 850 µm emission with the JCMT, and the IRAM 30m 1.2 mm data from Motte et al. (2007). The beam (main + error) of the JCMT data was taken from Di Francesco et al. (2008), and the beam (main + error) of the IRAM 30m telescope at 1.2 mm was taken from Greve et al. (1998). Three regions for which the 450 µm data from Di Francesco et al. (2008) were not usable were downloaded from the JCMT archive (IC1396N: m98bu24; AFGL5142: m96bc57; W3IRS5: m02ac32). The single-dish images at 450 and 850 µm (1.2 mm) for the 19 cores, together with the interferometric images showing the fragmentation in each core, are shown in the Appendix .
In order to characterize the spatial structure of the massive dense cores, we computed the circularly averaged radial intensity profile, in rings of 4 ′′ width, as a function of the projected distance from the core center at both 850 and 450 µm when possible, and applied the model described in Fig. 1 and Section 2.2. The intensity profiles at 850 and 450 µm are shown in Figs. 2-5. As can be seen in the figures, the emission is partially resolved for all the cores, revealing the exis- tence of an envelope of decreasing intensity with radius, even at small distances ( 10 ′′ ) from the core center. In some cases, the radial profile was calculated in particular angular 'sectors' centered on the core peak, to avoid strong negative emission in the maps (NGC7129-FIRS2 at 450 µm), or contamination from a nearby strong core (W3IRS5) or from the DR21 ridge (CygX-N48, CygX-N53, DR21-OH). For these cases we did not use the values of the flux density at 850 and 450 µm listed in Di Francesco's catalog, but remeasured the flux densities to assure that we were not including emission from the nearby cores/ridge. For the case of W3IRS5, the sector was taken in the north-south direction also to avoid possible contamination of free-free emission from strong centimeter sources 15 . For the case of CygX-N48, we extracted the intensity profile avoiding the contribution from the north, where the strong DR21-OH core lies; for DR21-OH and CygX-N53 we extracted the profile in a sector with no restrictions in the east-west direction, but allowing it to extend ∼ ±25 ′′ in the north-south direction, which corresponds to the approximate width of the DR21 filament in the east-west direction.
The SED of each massive dense core was built from centimeter wavelengths up to 50-60 µm, as this is approximately the wavelength where envelope and disk contributions start to be comparable (e. g., Crimier et al. 2010;Commerçon et al. 2012), and here we only model the envelope contribution. To build the SEDs we used the data compiled in Palau et al. (2013), which included flux densities at different frequencies from IRAS, JCMT/SCUBA, AKARI and IRAM 30m data when available. In addition, we measured the flux densities at 70 µm of cores with available Spitzer/MIPS2 data. Flux densities obtained in heavily saturated Spitzer/MIPS2 data were considered as lower limits. We also considered as lower limits the millimeter/submillimeter flux densities measured with an interferometer, as part of the core emission is filtered out and we did not take this into account in our model. For OMC-1S we recompiled the SED because the SED in Palau et al. (2013) was built for the source 136−355 (about ∼ 12 ′′ north of the submillimeter peak), while now we have considered the fluxes associated with the submillimeter peak, as we are modeling the submillimeter profiles. In addition, we used the FIR/submillimeter/centimeter fluxes from the literature 16 . 15 Although Richardson et al. (1989a) estimated the contamination from free-free emission to be < 25% at 1.1mm, and only < 3% at 800 µm, the centimeter sources A and B reported by Tieftrunk et al. (1997)  The apertures where fluxes were assessed were figured out from the literature or adopted as reported in the used catalogs and references (see above), and for peak intensities we took the beam as aperture radius.

Model
We modeled the thermal emission of the cores at millimeter and submillimeter wavelengths with a spherically symmetric envelope of gas and dust, with density and temperature decreasing with radius. The density and temperature were assumed to be power-laws of radius, with indices p and q, ρ = ρ 0 (r/r 0 ) −p and T = T 0 (r/r 0 ) −q . For the dust opacity law we assumed that the dust opacity is a power-law of frequency with index β, κ ν = κ 0 (ν/ν 0 ) β , where ν 0 is an arbitrary reference frequency. The values used for the dust opacity were those of Ossenkopf & Henning (1994), for thin ice mantles at a gas density of 10 6 cm −3 , i.e. κ 0 = 0.008991 cm 2 g −1 for ν 0 = 230 GHz, being β a free parameter of the model.
In the case of optically-thin dust emission, in the Rayleigh-Jeans approximation, and for an infinite radius envelope, the intensity as a function of the projected distance to the source center, I ν (b), can be derived analytically (see, for instance, Beltrán et al. 2002a), and is a power-law of the projected distance to the envelope center b, In our model, however, none of the above approximations (optically thin, Rayleigh-Jeans) were made, and there is no analytical expression for the intensity profile. In general, the contribution to the radiated intensity from an element of the envelope of geometrical depth dz, at a radius r from the center (see Fig. 1) is given by with r = √ b 2 + z 2 , is the optical depth from the observer to the element of the envelope at a projected distance b and position z along the line of sight (see Fig. 1), and J ν (T ) is the Planck function in units of temperature, Then, the resulting intensity as a function of the projected distance to the source, I ν (b), has to be evaluated as where z max = b 2 + r max 2 and r max is the maximum radius where the model is calculated. The computational burden is  Table 1 for the exact fitted parameters). Each row corresponds to one core, the left (middle) panel shows the radial intensity profile at 850 (450) µm, with the empty blue circles corresponding to the data, the black solid line corresponding to the model, and the dashed red line showing the beam profile; panels on the right show the SED, with blue empty circles indicating the observed fluxes, the blue full circles indicating observed lower limits, the black solid line showing the model for a fixed aperture, and the red squares corresponding to the model for the same aperture where each flux was measured. low, since both integrals, those of Eqs. 2 and 4, can be evaluated simultaneously in a single run.
Regarding the temperature distribution of the envelope, Chandler et al. (1998) derive from the expressions of Scoville et al. (1976) that, for the radiative heating of a dust cloud by a central luminous source, the temperature power-law index q is a function of the dust emissivity index β, q = 2/(4 + β). However, for the outer part of the envelope the external heating can be important, and a constant temperature of 10 K was adopted for radii larger than the radius where the temperature distribution of index −q drops to 10 K. Regarding the density distribution, we adopted, as maximum radius of the envelope r max , the radius for which the envelope density dropped to a value comparable to that of the ambient gas of the intercore medium, taken to be 5000 cm −3 .
In summary, the four free parameters of the model are the dust emissivity index, β; the envelope temperature at the reference radius r 0 , T 0 ; the envelope density at the reference radius  (Table 1). r 0 , ρ 0 ; and the density power-law index, p. The reference radius r 0 is arbitrary, and has been taken as r 0 = 1000 AU. With this model we fitted the observed radial intensity profiles at 450 and 850 µm or 1.2 mm, and the observed SED from centimeter wavelengths up to 60 µm simultaneously. The uncertainty adopted for the intensity profiles was 20% at 850 µm and 1.2 mm, and 50% at 450 µm (Di Francesco et al. 2008). The uncertainty of the flux densities used to build the SED was adopted to be 50%.
In order to compare the model with the observed intensity profile at each wavelength, we computed the intensity map from the model, and we convolved it with the beam, given as the sum of two circular Gaussians (the main and error beams). Finally, the circularly averaged profile was computed from the convolved map. We also simulated the effect of the finite chop throw used in the observations to cancel the sky emission. Thus, for each map position we estimated the off-source intensity as the average intensity at a chop-throw distance of 120 ′′ from the map (on-source) position, and we computed the difference between the on-and off-source intensities. From this, we computed the radial profile corrected for chopping.
In order to compare the model with the observed SED, for each data point with wavelength higher than 60 µm we computed the model flux density by integration of the model in-  (Table 1).
tensity profile inside the reported aperture used for each measurement. In addition, for points with frequency below 1000 GHz, we added the free-free emission of an unresolved source with a flux density extrapolated from the data points at centimeter wavelengths.
The fitting procedure was the sampling of the fourdimensional parameter space, using the same procedure as that described in Estalella et al. (2012) and Sánchez-Monge et al. (2013). For every set of parameters, we computed the residual χ 2 , where the sum encompasses the points of the intensity profiles at 850 and 450 µm, and the points of the SED fitted by the model. The parameter space was searched for the minimum value of χ 2 . For each source, the initial search ranges for the four parameters was β = 1.5 ± 1.5, T 0 = 300 ± 300 K, ρ 0 = (1.0 ± 1.0) × 10 −16 g cm −3 , and p = 1.5 ± 1.0. The search range was reduced a factor of 0.8 around the best fit value of  Table 1. For CygX-N12 and CygX-N63 the profiles were extracted from the 1.2 mm images of Motte et al. (2007). the parameters found for each loop, consisting in 2000 samples of the parameter space. The final best fit values were taken after 10 loops. Once the best fit parameters were found, their uncertainties were estimated through the increase in χ 2 (see Estalella et al. 2012. We assessed the quality of the fits through the value of the reduced chi-square, In general for all the sources we found fits with χ r ≃ 1, or even less, indicating that the errors of the data were slightly overestimated. In Fig. 14 we show as example, the parameter space around the best fit solution for the source IC 1396N.

RESULTS
Figures 2-4 present the best simultaneous fit to the 850 and 450 µm (or 1.2 mm) radial intensity profiles and the SED for the 19 massive dense cores. A summary of the fitted parameters is presented in Table 1. The four fitted parameters span a wide range of values. The dust emissivity index ranges from 1.0 to 1.9, with a mean value of 1.5. This is comparable to dust emissivity indices measured in high-mass envelopes (e. g., Mueller et al. 2002;Hill et al. 2006 for temperatures between 30 and 50 K), and corresponds to temperature powerlaw indices ranging from 0.34 to 0.40 (mean value: 0.37, consistent with Emerson 1988). The fitted values for the temperature at 1000 AU range from 40 to 260 K, with a mean value of 83 K, and correlate with the bolometric luminosity 17 . The range of densities at 1000 AU is (3-30)×10 −17 g cm −3 , with a mean value of ∼ 8 × 10 −17 g cm −3 . Finally, the index of the density power law ranges from 1.5 to 2.4, with a mean value of 1.85, similar to the results found in previous works focused on massive dense cores (e. g., Pirogov et al. 2009). The cores with steeper density profiles are I22198, NGC 7129, and I20126, while the cores presenting the flattest profiles are OMC-1S, I05358NE, W3IRS5 and CygX-N12.   Palau et al. (2013). N mm is the fragmentation level, taken from Palau et al. (2013) and assessed by estimating the number of millimeter sources within a field of view of 0.1 pc of diameter as detected with a millimeter interferometer with a mass sensitivity of ∼ 0.3 M ⊙ and a spatial resolution of ∼ 1000 AU. For the four Cygnus X cores, N mm has been revised with respect to the values given in Palau et al. (2013). b Free parameter fitted by the model: β is the dust emissivity index; T 0 and ρ 0 are the temperature and density at the reference radius, 1000 AU; p is the density power law index; χ r is the reduced χ as defined in equation (6). Once the model was fitted to the data, we inferred different quantities, which are given in Table 2. The radius where the density reaches ∼ 5000 cm −3 , considered as the ambient density, ranges from 0.1 to 1.0 pc. The regions with smallest radii were I22198, NGC7129-1 and CygX-N63, in part due to the concentrated density profile inferred for these regions, of p 2. On the other hand, the cores with largest radii are I05358NE and CygX-N48. We estimated the 'observed mass' of the core by integrating the density profile up to the observed radius of the core, which ranged from 5 to 1000 M ⊙ . In addition, we estimated the surface density (Σ 0.05 pc ) and the density (n 0.05 pc ) of each core within a radius of 0.05 pc, as this is the radius where we assessed the fragmentation level. We found that Σ 0.05 pc ranged from ∼ 0.1 up to 1.8 g cm −2 , with a mean value of 0.6 g cm −2 . Finally, the radii where the temperature falls down to 10 K ranged from 0.2 to 18 pc, being the cores with higher temperature (at reference radius) those where this radius is larger.

Comparison of density structure with previous works
We compared our inferred density power-law indices to the results of previous works where a similar methodology was applied for some of the cores of our sample. In the two last columns of Table 1 we list p lit , the density power-law index from the literature, together with the corresponding references. Crimier et al. (2010), van der Tak et al. (2000van der Tak et al. ( , 2013 and Mueller et al. (2002) fit one single power-law to the radial intensity profiles, being comparable to the work presented here, while Beuther et al. (2002) used two power-laws with a break at 32 ′′ . Within the regions fitted by the same authors the tendency of flat and concentrated profiles is consistent with the results found in the present work (for example, IC1396N is flatter than NGC7129-FIRS2 and CB3-mm as in Crimier et al. 2010; or I20126 is more concentrated than I22134 and I05358NE as in Beuther et al. 2002). However, the absolute value of the density power-law indices cannot be directly compared between different works because of the different methods applied, which yield differences of 0.2-0.4 in the density power-law index. The most important difference between previous works and this work is the way of determining the temperature profile. While Mueller et al. (2002), 1.2 c 8 a Parameters inferred (not fitted) from the modeling. q is the temperature power-law index, and r 10 K is the radius of the core where the temperature has dropped down to ∼ 10 K; r max is the radius at the assumed 'ambient' density of 5000 cm −3 ; M obs is the mass computed analytically from the model integrating until the radius where the density profile could be measured for each source. Σ 0.05 pc and n 0.05 pc are the surface density and density inside a region of 0.05 pc of radius. b β rot is the rotational-to-gravitational energy ratio as in Palau et al. (2013) and adding the new measured values from NH 3 (1,1) for OMC-1S (Wiseman & Ho 1998) and DR21-OH (Mangum et al. 1992). σ noth−init is the non-thermal velocity dispersion for a quiescent nearby core (thus a first approach to the initial velocity dispersion), as in Palau et al. (2013) and adding a new value for CygX-N48 and DR21-OH (Mangum et al. 1992). σ noth−init was estimated assuming a temperature of 20 K. σ noth−feedb is the non-thermal velocity dispersion for the studied core (thus including outflow feedback from the nascent protostars) from data in the literature (see Section 4.4 for references). σ noth−feedb was estimated using the temperature at radius 0.05 pc as inferred from the modeled temperature profile (Tables 1 and 2). M vir and α vir = M vir /M 0.05pc (with M 0.05pc being the mass within a radius of 0.05 pc) are assessed using the non-thermal velocity dispersion including feedback from outflows. c 'Turb. mode' refers to tentative turbulence mode according to works in the literature: 'c' indicates that a detailed study has been done comparing the properties of the region with simulations or studying the Probability Density Function; 'c?' indicates that the only hints of possible compressive turbulence are that the region lies at the border of an expanding bubble/HII region. References about turbulence mode: 1: Beltrán et al. (2002b) Crimier et al. (2010), and van der Tak et al. (2013) calculate the temperature profile self-consistently through radiative transfer codes, we assume that the temperature profile is related to the dust emissivity index (Section 2.2) and leave the temperature at a given radius as a free parameter. This is also different to Beuther et al. (2002) work, who assumed a fixed temperature power-law index of 0.4. Our method finds in general flatter temperature profiles compared to Crimier et al. (2010) and van der Tak et al. (2013), yielding higher density power-law indices. On the other hand, our method is perfectly consistent with the results of Mueller et al. (2002) for AFGL 2591. Therefore, once the different methodologies are taken into account, our methodology yields consistent results with previous works, and allows to fairly compare the fitted parameters from one region to the other.

Uncertainty in the fragmentation level and possible biases in the sample
Before comparing the inferred density structure for each massive dense core (estimated through the density power-law index p), to the fragmentation level (estimated through the N mm parameter, see Table 1 and Palau et al. 2013), we discuss here the main uncertainties in the assessment of these properties and possible biases. First of all we stress that the determination of both p and N mm has been carried out applying the same methodology and criteria to all the cores and thus, although the absolute values might be subject to discussion, the relative value between the regions is meaningful and valid for a comparison. In particular, N mm is a lower limit to the total number of fragments because of flux missed by the interferometers, and because we are missing objects emitting at wavelengths different to millimeter. However, the sample presented in Palau et al. (2013) was selected applying strict criteria for the filtered emission and uv-coverage (see column 8 of Table 3 of Palau et al. 2013), constraining this study to the most compact fragments ( 5000 AU). Since larger non-centrally peaked structures might be transient, this is not Fig. 6.-Relation between the fitted parameter p (density power-law index, panel "a"), ρ 0 (density at a radius of 1000 AU, panel "b"), and n 0.05 pc (density within 0.05 pc of radius, panel "c") with N mm (fragmentation level, see Table 1). In the last panel, the Jeans number is calculated as M 0.05pc × CFE/M Jeans , assuming a temperature of 20 K and a Core Formation Efficiency (CFE) of 17% (see main text), and the green stars indicate the measurements in the simulations of Commerçon et al. (2011) for the µ = 2 (strongly magnetized) and µ = 130 (weakly magnetized) cases. In all panels, red dots correspond to cores with hints of compressive turbulence, and blue dots correspond to cores with no assessment of the turbulence mode in the literature (see Table 2 and Section 4.4).
likely affecting our estimate of the fragmentation level. In addition, the number of infrared sources were compared to N mm and were found to be correlated, indicating that N mm is systematically underestimating the 'real' fragmentation level of the core by only a factor of ∼ 2. Finally, Palau et al. (2013) also used radiation MHD simulations to show that N mm is sensitive to the original fragmentation of the core. Altogether suggests that N mm is already a good indicator of the fragmen-tation level if used to compare the different regions and not as an absolute value.
A possible bias could come from the different distances in the sample. However, for 16 out of 19 cores distances range only from 0.8 to 2.5 kpc, and no trend of N mm , p, or any of the fitted parameters vs distance is found in our data (Table 1), indicating that this is not affecting our conclusions. Another bias could arise from different evolutionary stages of the cores of our sample, as simulations show that the degree of substructure in a cloud decreases with age (e. g., Smith et al. 2009). Our sample is not likely dramatically affected by this because the targets were selected to be in similar evolutionary stages: all targets harbor an intermediate/high-mass protostar detected in the far-infrared, are associated with a strong and compact millimeter core, are associated with molecular outflows, and have not developed an (ultra-)compact HII region yet. The small differences in evolutionary stage of the cores of this sample were also studied in Palau et al. (2013) by estimating evolutionary indicators, and no relation between N mm and these indicators was found (see e.g., Fig. 8B of .
Therefore, no biases in distance, evolutionary stage or missing fragments in N mm should be strongly affecting the results presented in this work.

Trend of fragmentation level with density power-law index and density
In this Section we study the relation between the core parameters inferred from the model and the fragmentation level, quantified through the number of millimeter sources within a field of view of 0.1 pc (in diameter), N mm , as in Palau et al. (2013). We find no trend of β, T 0 and total mass of the core with the fragmentation level (Tables 1 and 2). Conversely, Fig. 6 presents the relation between the density power-law index p, the central density (at 1000 AU) ρ 0 , and the density within a radius of 0.05 pc, n 0.05 pc , with the fragmentation level. Fig. 6-a reveals a possible weak trend of N mm with the density power-law index, with low fragmentation levels found mainly in cores with steeper density profiles. Concerning N mm vs ρ 0 , our sample splits up into two groups, cores with high fragmentation and high ρ 0 , and cores with lower fragmentation and low ρ 0 , suggesting that high central densities (in opposition to average density, as assessed in , could be related to the fragmentation level. Finally, there is a clear trend of high fragmentation levels in cores with high n 0.05 pc (Fig. 6-c). A closer inspection of the data using statistical tests, like the Spearman's ρ and Kendall's τ rank correlation coefficients, indicates that N mm and p show a faint anti-correlation (ρ ∼ −0.49, τ ∼ −0.3, Pearson's linear correlation coefficient R ∼ −0.3). As for n 0.05 pc , it is strongly correlated with N mm (ρ ∼ 0.70, τ ∼ 0.54), and the two parameters are well described by a linear relationship (R ∼ 0.90). These relations indicate that the high n 0.05 pc likely arise in cores with rather flat p and also high ρ 0 . A trend of high fragmentation level for flat density profiles was already suggested theoretically by Myhill & Kaula (1992). More recently, Girichidis et al. (2011) show, by studying the evolution of collapsing non-magnetized turbulent massive cores for different initial density profiles, that cores with a flat density profile produce a high number of protostars, while concentrated density profiles rather yield one massive protostar in the core center. This can be understood in the framework of Jeans fragmentation, as the larger the density, the smaller the Jeans mass (for a given tem-perature) and the larger the number of expected fragments N Jeans . In Fig. 6-c we overplot the expected number of fragments estimated as N Jeans = (M 0.05pc × CFE)/M Jeans , with M Jeans = 0.6285 [T/10K] 3/2 [n H 2 /10 5 cm −3 ] −1/2 M ⊙ , for which we used T ∼ 20 K, and M 0.05pc being the mass within a radius of 0.05 pc (calculated directly from the density law inferred from our model). The Core Formation Efficiency (CFE) was adopted to be 17%, the average of the CFEs measured in Palau et al. (2013). Thus, the trend of higher fragmentation for higher densities within a radius of 0.05 pc is consistent with what is expected for Jeans fragmentation. The relation found in this work between fragmentation level (which could be regarded as a proxy to cluster richness) and gas density is reminiscent of the relation found on larger spatial scales between cluster mass or protostellar density and gas surface density (e. g., González-Lópezlira et al. 2012;Lombardi et al. 2013;Pflamm-Altenburg et al. 2013) and further studies should be conducted in this direction.

Physical processes affecting the density structure of
Massive Dense Cores With the aim of shedding light on which physical conditions promote fragmentation in massive dense cores, in this Section we study if there is any relation between different physical processes and the density structure (density profile and central density) of massive dense cores.
Time evolution -It has been shown both from observations and numerical simulations that structures within collapsing molecular clouds become more centrally peaked with time (e. g., Smith et al. 2009;Giannetti et al. 2013). Therefore, the cores for which we measure higher values in the density power-law index p could just correspond to the most evolved ones. However, we searched for any relation between p and the ratio of the bolometric luminosity to the single-dish mass parameter (an evolutionary indicator used in Palau et al. 2013 18 ), and found no trend in our data ( Fig. 7-a). This lack of correlation of p with time probably results from the fact that our sample consists of cores at very similar evolutionary stages (Section 4.2), not properly sampling broad enough timescales, and suggests that other processes in addition to time evolution are probably also playing a role in determining the density structure of a dense core.
Rotational-to-gravitational energy (β rot ) -The rotation of a massive dense core is expected to promote flattening of its density profile. In order to test this in our sample, we completed the values of rotational-to-gravitational energy ratios already estimated in Palau et al. (2013) with the value for the core added in this work (DR21-OH), observed by  with similar conditions of sensitivity and spatial resolution as the sample of Palau et al. (2013). Here we estimated β rot for DR21-OH following exactly the same method as in Palau et al. (2013), i. e., using NH 3 (1,1) data (Mangum et al. 1992), and obtained a value of 0.24. We also estimated β rot for OMC-1S from the NH 3 (1,1) data by Wiseman & Ho (1998). It is worth mentioning that these measurements of β rot must be regarded with caution as the velocity gradients observed might not be necessarily associated with rotation, but could arise from multiple velocity components or infall motions (e. g., Lee et al. 2013). Keeping this in mind, in Fig. 7  Relation between the bolometric luminosity to the single-dish mass parameter L bol /M sd and the density power law index. b) Relation between rotational-to-gravitational energy β rot and the density power law index. c) Relation between non-thermal velocity dispersion (including outflow feedback) and density power-law index. Both plots include new cores after Palau et al. (2013, see Section 4.4). In all panels, red dots correspond to cores with 'confirmed' compressive turbulence or with hints of compressive turbulence, and blue dots correspond to cores with no assessment of the turbulence mode in the literature (see Table 2 and Section 4.4).
we present the relation between β rot and p, showing a weak inverse trend. This suggests that the initial angular momentum alone might not be the only agent shaping the density profile of a core. We also present an updated version of the plot of N mm vs β rot (after Palau et al. 2013) including these new measurements, and the updated plot, shown in Fig. 8, indicates a possible trend of β rot with fragmentation level, which deserves further studies.  and the rotational-togravitational energy ratio, after Palau et al. (2013, i.e., including new measurements for OMC-1S and DR21-OH). Red dots correspond to cores with 'confirmed' compressive turbulence or with hints of compressive turbulence, and blue dots correspond to cores with no assessment of the turbulence mode in the literature (see Table 2 and Section 4.4).
Non-thermal velocity dispersion (σ non−th ) - Palau et al. (2013) estimate the non-thermal velocity dispersion from the linewidth measured in NH 3 (1,1) of a quiescent massive dense core in the immediate surroundings of the core where fragmentation was assessed. This can be regarded as a first approach to the initial turbulent dispersion. We completed the values of non-thermal 'initial' velocity dispersions by adding two new measurements (CygX-N48 and DR21-OH by Mangum et al. 1992, Table 2). However, the non-thermal velocity dispersions cover a small range of values, from 0.33 to 0.67 km s −1 , and did not show any relation with the fragmentation level. We additionally study the relation of nonthermal velocity dispersion of the fragmenting cores (thus including feedback from outflows), covering a broader range from 0.41 to 2.16 km s −1 , with the density profile (Table 2, where we used the values from NH 3 (1,1) observed with the VLA: Torrelles et al. 1989;Zhou et al. 1990;Mangum et al. 1992;Tieftrunk et al. 1998;Wiseman & Ho 1998;Gómez et al. 2003;Sánchez-Monge et al. 2013). Since no clear trend is found (see Fig. 7-c and Table 2), this indicates that probably the density profile and the turbulent velocity dispersion are independent quantities. Finally, we estimated the virial masses and virial parameter (M vir /M 0.05pc ) using the non-thermal velocity dispersions of the fragmenting cores (i. e., including feedback), and found no relation of these quantities with the density power-law index or the fragmentation level.
Turbulence mode: compressive vs solenoidal -In order to study if different turbulence modes could be affecting the density structure of massive dense cores, we searched the literature for any evidence of our massive dense cores being affected by expanding shells/bubbles created by nearby O-type stars, and our findings and references are listed in the two last columns of Table 2. Regions with hints in the literature of being associated with expanding shells or HII regions, such as brightrimmed clouds, etc., or located at the meeting point of converging flows, are classified as 'possible compressive', while regions for which the probability density function has been found to be consistent with compressive turbulence are classified as 'confirmed compressive'. Both 'possible' and 'confirmed' compressive-turbulent cores (see Table 2 to distinguish among both types) are assigned a 'red' color in the figures, while the 'blue' color corresponds to cores where  Commerçon et al. (2011). Note that for the magnetized case the core has a more concentrated profile than for the weakly magnetized case. The r −3 and r −1.5 profiles (dashed) are shown to guide the eye. the turbulence mode cannot be assessed with the current information (thus a solenoidal mode cannot be ruled out for 'blue' cores). Fig. 6 shows that 'compressive-turbulent' cores tend to have large N mm , and cover a wide range of density power-law indices (from 1.5 to 2). In addition, cores with ρ 0 > 10 −16 g cm −3 and n 0.05 pc > 6 × 10 5 cm −3 fall all within the 'compressive-turbulent' group (see Fig. 6-b,c). This could be consistent with compressive turbulence because this mode yields higher density contrast structures than the solenoidal mode (e. g., Federrath et al. 2008Federrath et al. , 2010. Recently, Tremblin et al. (2013) find that expanding HII regions or bubbles tend to induce steep density profiles in the cores at the bubble edges. Thus, the steep density profiles obtained in DR21-OH, AFGL 5142, CB3 (p ∼ 2) could have been produced by the expanding ionization fronts from the massive stars in these regions (Hunter et al. 1995;Schneider et al. 2011;Sakai et al. 2013). A more detailed and uniform study of the turbulence mode in each core would be required to definitely assess the role of the different turbulence modes in setting the density structure of massive dense cores.
Magnetic field -Magnetic field could affect the density structure of a core because it helps removing angular momentum yielding rather steep density profiles. We tested this by using the magnetohydrodynamic simulations from Commerçon et al. (2011) of strongly (µ = 2, with µ being the mass-toflux over the critical mass-to-flux ratio) and weakly (µ = 130) magnetized cores. These simulations were found to reproduce remarkably well our data . It is important to note that in these simulations the initial density profile is the same and rather flat for both µ = 2 and µ = 130 cases. We extracted the density profiles of Commerçon et al. (2011) simulations for µ = 2 and µ = 130 cases, by making an average of the density on concentric spherical shells weighted by the mass for each cell in the shell. The result is displayed in Fig. 9, and shows that the µ = 2 case evolved into a more concentrated profile than the µ = 130 case, for radii up to ∼ 0.05 pc or 10000 AU. This results from the fact that for the low magnetic field case (µ = 130), the fragmented area is large because of conservation of angular momentum, which gives a flatter density distribution at the center. And the opposite, strong magnetic braking (in the µ = 2 case) helps concen-trate matter in the center and gives a steeper profile. In addition, previous theoretical work including ambipolar diffusion (e. g., Lizano & Shu 1989) also predict steep density profiles (density power-law indices of ∼ 2). Thus, the current theoretical and numerical works suggest a relation between the inferred density power-law index of massive dense cores and the magnetic field properties. Finally, we also used the simulations from Commerçon et al. (2011) to estimate the density within a radius of 0.05 pc for both µ = 2 and µ = 130, and find n 0.05 pc ∼ 2.3 × 10 5 cm −3 and 5.8 × 10 5 cm −3 , respectively. This, together with the number of millimeter sources that a typical interferometer would detect in each simulation (estimated in Palau et al. 2013), allowed us to compare the results of these simulations to our N mm vs n 0.05 pc plot showing observational data (see Fig. 6-c). The values predicted by the simulations match remarkably well the trend seen in the observational data, suggesting that the magnetic field strength might be related to n 0.05 pc of massive dense cores, with small densities found in the strongly magnetized cores and large densities found in weakly magnetized cores. Overall, the agreement between Commerçon et al. (2011) simulations and our data is very good and suggests that magnetic field might be affecting the shape and density of massive dense cores, which in turn seems to be related to the fragmentation level.

CONCLUSIONS
In order to study the relation between the density profile of massive dense cores and the fragmentation level, we have fitted simultaneously the radial intensity profiles at 450 and 850 µm and the SED of 19 massive dense cores for which their relative fragmentation level was known from previous works. Although the parameter used to estimate the fragmentation level is a lower limit, its relative value for the different cores is significant because it is assessed in a uniform way for all of them. The sample consists of massive dense cores undergoing active star formation and being at similar evolutionary stages. We modeled the cores as spherical envelopes with density and temperature decreasing with radius following power-laws, and found a weak trend of lower fragmen-tation for cores with steeper (p 2.0) density profiles. A correlation is found between the fragmentation level and the density within a radius of 0.05 pc (the radius where fragmentation was assessed), suggesting that high fragmentation levels might be related with flat cores and cores with high central densities. This correlation is consistent with Jeans fragmentation, as high densities yield lower Jeans masses.
While we find no clear relation between velocity gradients or velocity dispersion and the density power law index of our massive dense cores, cores with compressive turbulence seem to be associated with high central densities, thus promoting fragmentation. On the other hand, a comparison with radiation magnetohydrodynamical simulations suggests that a strong magnetic field could help shape steep density profiles in massive dense cores, for which low fragmentation levels are found. Thus, it is of crucial importance to study the relation of density profile, turbulence mode and magnetic field from both an observational an numerical point of view.  (Jenness et al. 1995;Di Francesco et al. 2008). For each source, the two green squares are ∼ 0.1 × 0.1 pc 2 , the size of the mm interferometer map (right panel, Palau et al. 2013) revealing the fragmentation in each core. For IC1396N, NGC 7129-FIRS2, and CB3-mm, contours of 850 (450) µm emission (left (middle) panels) are 8% to 99%, increasing in steps of 10% of the peak value, 5.0 (25.3), 4.1 (48.9), and 2.1 (12.8) Jy beam −1 , respectively. For I22198, contours of 850 (450) µm emission (left (middle) panel) are 20% to 99%, increasing in steps of 20% of the peak value, 5.5 (9.8) Jy beam −1 . For NGC 2071, contours of 850 (450) µm emission (left (middle) panel) are 5% to 99%, increasing in steps of 10% of the peak value, 13.1 (97.3) Jy beam −1 . Contours for interferometric images (right panels) are −4, 4, 8, 16, 32 and 64 times the rms noise of each map as given in Palau et al. (2013).