© The Authors 2024

1 Introduction

Most high-mass stars are known to form in clusters; however, it is not yet clear to what extent different physical processes contribute to the fragmentation of the parental star-forming regions. One wants to know, for example, whether thermal Jeans fragmentation is the main mechanism (e.g., Palau et al. 2013, 2014, 2015) or if turbulence plays a significant role (e.g., Wang et al. 2014). Other questions are understanding the importance of the initial conditions, in particular the density structure (e.g., Girichidis et al. 2011); whether thermal feedback can suppress fragmentation (e.g., Krumholz et al. 2009); and in what way the magnetic field influences the fragmentation of the cluster-forming regions (e.g., Commerçon et al. 2011, 2022; Peters et al. 2011; Myers et al. 2013, 2014; Zhang et al. 2014; Federrath 2015; Beuther et al. 2018b, 2020; Palau et al. 2021).

To study the fragmentation of high-mass star-forming regions, we embarked on the large program CORE¹ with the IRAM NOrthern Extended Millimeter Array (NOEMA; formerly called the Plateau de Bute Interferometer, PdBI) that is investigating 20 well-known high-mass star-forming regions in the northern hemisphere at high spatial resolution (~0.3″–0.4″) in the 1.3 mm continuum and spectral line emission (Beuther et al. 2018a). The sample was selected to have high luminosities > 10⁴ L_⊙, distances of less than 6 kpc, and high declination to ensure good uv-coverage with northern hemisphere interferometers. While the continuum data are used for the fragmentation analysis, the spectral line data allow us to investigate the turbulent, kinematic, and physical properties, as well as the chemical characteristics of the regions (e.g., Beuther et al. 2018a; Gieser et al. 2021; Ahmadi et al. 2023).

Investigating the 1.3 mm continuum data of the sample, Beuther et al. (2018a) distinguished a large variety of fragmentation properties ranging from sources that are dominated by a single massive core to those that fragment into as many as 20 cores (independent of distance). The sample was selected uniformly with respect to their evolutionary phase, they are all high-mass protostellar objects (HMPOs; see Beuther et al. 2018a for the indicators); some also harbor hot molecular cores and ultracompact HII regions. Hence, evolutionary effects are unlikely to explain the results, and the diversity from highly fragmented to barely fragmented regions is a real outcome of this study (Beuther et al. 2018a).

The challenge is to isolate which of the above physical processes dominate the fragmentation diversity. Conducting a minimum spanning tree analysis, Beuther et al. (2018a) find that typical nearest neighbor separations are below the thermal Jeans fragmentation scale, indicating that turbulence cannot explain the observed core separations. Since all regions are at similar evolutionary stages, thermal feedback is also unlikely to cause the fragmentation differences. Two possible explanations are differences in the magnetic field and/or differences in the density structure of the parental gas clump. We also investigated whether the interferometric spatial filtering varies for regions with different fragmentation properties because steep density profiles should result in less missing flux than flatter density distributions. However, no considerable spatial filtering differences were found within the sample (Beuther et al. 2018a). In the following we refer to clumps as the parental star-forming regions on parsec scales, and to cores as the fragments on sub-0.1 pc scales.

In this study we investigate the density distribution of the parental gas clumps and the magnetic field structures of the regions. While the density structures are studied by means of large-scale single-dish 1.2 mm dust continuum mapping of the regions with the IRAM 30 m telescope, the magnetic fields are investigated with polarization observations of the dust emission with the Submillimeter Array (SMA) at 875 µm. Overall parameters for the investigated CORE sample are presented in Table 1.

Density structure. Density distributions ρ of star-forming regions can be described via radial r power-law distributions, such as ρ ∝ r^−p. To first order, flatter density distributions may result in more fragments, whereas steeper density distributions favor less fragmentation (e.g., Girichidis et al. 2011). Observational studies of low-mass cores typically found density distributions with p varying between 1.5 and 2, resembling finite-sized Bonnor-Ebert spheres (e.g., Motte & André 2001; Alves et al. 2001). While early single-dish studies in the high-mass regime of ultracompact HII regions indicated potentially shallower slopes that may have favored logatropic equations of state (e.g., van der Tak et al. 2000a; Hatchell et al. 2000), later studies of younger evolutionary stages again found density distributions with exponents p mostly between 1.5 and 2 (Beuther et al. 2002; Mueller et al. 2002; Hatchell & van der Tak 2003; Palau et al. 2014). Hence, these studies indicated that the density distributions between low- and high-mass star-forming regions should not be very different. Since then, (very) few observational studies have focused on clump scales to constrain the density structures of the whole cluster-forming regions. Although the interstellar medium is known to be filamentary (e.g., André et al. 2014; Molinari et al. 2016), investigating structures on clump scales still allows us to derive the physical properties of individual (high-mass) star-forming regions. On the theoretical side, Shu (1977) already determined that the singular isothermal sphere solution results in density profiles with p = 2. Analytic and numerical modeling showed that p = 2 is an attractor, which means that initially flatter density distributions approach p = 2 during the star formation process (e.g., Naranjo-Romero et al. 2015; Gómez et al. 2021).

While the density structures of fragmented cores for the CORE sample have been studied in detail by Gieser et al. (2021), their large-scale density structure from the parental clumps is not well constrained yet (with the exceptions of IRAS 23033 and IRAS 23151 studied in Beuther et al. 2002). Studying now the large-scale distributions, we can directly set them into context with the observed fragmentation properties (e.g., number of fragments) and with the small-scale density structures reported in Gieser et al. (2021).

Magnetic field structure. Six of the CORE sample regions have been observed at single-dish resolution (20″) with the SCUBA polarimeter SCUPOL (Matthews et al. 2009). However, all these data show comparably weak polarization signal toward the Stokes I peak positions (called polarization holes), mainly caused by unresolved magnetic field structures in the single-dish beam (cf. W3(H₂O) in Matthews et al. 2009 and Chen et al. 2012). Thus, higher spatial resolution observations are imperative. For three of the regions (W3(H₂O), NGC 7538IRS1, and NGC 7538S), SMA magnetic field studies have already been published (Chen et al. 2012; Frau et al. 2014; Palau et al. 2021). While NGC 7538IRS1 shows barely any fragmentation in the CORE study (Beuther et al. 2018a, Table 1), and high magnetic field strengths on the order of ~2 mG in the SMA data (Frau et al. 2014) are consistent with reduced gas fragmentation (e.g., Commerçon et al. 2011, 2022; Peters et al. 2011; Myers et al. 2013, 2014), W3(H₂O) and NGC 7538S show intermediate levels of fragmentation (seven and six fragments, respectively; see Table 1), again with magnetic field strengths in the milligauss regime (Chen et al. 2012; Palau et al. 2021). Therefore, with only three regions of the CORE sample observed in the past, the effect of magnetic fields remains inconclusive. Zhang et al. (2014) observed a sample of intermediate- to high-mass star-forming regions with the SMA, and inferred that magnetic fields are indeed important during the star-formation process. Recently, Palau et al. (2021) compared the same magnetic field data to the fragmentation properties of these regions, and they found a tentative correlation between the number of fragments and the mass-to-flux ratio. Summaries of the current state of high-resolution interferometric studies of star-forming regions can be found in Hull & Zhang (2019) or Liu et al. (2022). Clearly this is only the beginning of an exciting avenue to follow. Therefore, to further investigate the fragmentation and magnetic field properties in context, complementary high-spatial-resolution density and magnetic field data are needed. Since the 0.3″–0.4″ mm continuum data for the CORE sample have already been acquired, we now complement these with the magnetic field information from the SMA.

2 Observations

2.1 IRAM 30 m telescope: large-scale density structure

All 20 CORE regions were observed with the New IRAM KIDs Array 2 (NIKA2, Adam et al. 2018; Perotto et al. 2020²) on the IRAM 30 m telescope between December 2019 and February 2020 (project 143-19). NIKA2 consists of 2900 kinetic inductance detectors (KIDs) operating at 150 and 260 GHz. Since the main goal of this study is to derive reliable intensity profiles (and in turn the density profiles) of the central star-forming regions, the inherent spatial filtering in single-dish continuum observations has to be considered (e.g., Perotto et al. 2020, see also Sect. 3.1). Therefore, we followed the IRAM recommendation to observe extended map sizes of roughly 9.8′ per region. While smaller maps were also observed, here we use only the extended maps.

Data reduction was conducted with the Pointing and Imaging In Continuum (PIIC) software package³. Since all the regions have extended emission, for the best recovery of spatial structures, we extensively tested various source structures as well as mapping parameters within PIIC. For this purpose we first created artificial sources with different power-law profiles. These were then smoothed to the 12″ beam of NIKA2 at 1.2 mm wavelengths and run through the PIIC pipeline. We explored, for example, how much source masks affect the retrieval of the intensity profiles. We found that for our application the best results were obtained without applying any source masks since the emission is extended with complex substructures, such as dense clumps and filaments. Therefore, for our final dataset no source masks were defined and we applied baselines of order one. At most 35 iterations were used, but we always stopped further iterations when the peak flux density difference between the previous two iterations was less than 0.1% (for more details about source structure recovery, see also Sect. 3.1).

NIKA2 always observes simultaneously at 1.2 and 2 mm wavelengths. For our study, we require the highest spatial resolution. Therefore, we only present the 1.2 mm data. The angular resolution of the final maps is 12″ and the 1σ rms varies between 4 and 16 mJy beam⁻¹ (Table 2).

2.2 Submillimeter Array (SMA): magnetic field structure

The CORE sample was observed with the Submillimeter Array (Ho et al. 2004) over a period of several years. The first few regions were observed in the winter term 2018–2019 and the last sources in September 2021. Altogether, our observing campaign covered 17 CORE regions; the remaining three CORE regions had already been observed with the SMA in earlier campaigns (Sect. 1 and Table 1). Typically, we observed two regions together in track-sharing mode. Using the compact configuration in the 875 µm wavelength band, a spatial resolution between 2″ and 4″ can be achieved. The final synthesized beams for all regions are listed in Table 2.

The target sources were always observed in loops together with gain calibrators that are also sensitive to the polarization. Typical loop lengths were 15 min. Bandpass and flux calibrations were conducted either at the beginning or the end of the tracks. The bandpass calibrator was also used to derive the instrumental polarization of the Array. The SMA observations for the 17 newly observed targets cover in total 16 GHz of spectral band-with, 8 GHz each in the lower and upper sideband. For details of the remaining three regions (W3(H₂O), NGC 7538IRS1, and NGC 7538S) we refer to the respective publications by Chen et al. (2012), Frau et al. (2014), and Palau et al. (2021). The exact frequency coverage of the new observations was 338.85–346.88 GHz in the lower sideband and 354.84–362.87 GHz in the upper sideband. The native spectral resolution was 140 kHz per channel and we binned the data to a spectral resolution of 559 kHz, corresponding at 345 GHz to a velocity resolution of ~0.5 km s−1. The continuum data were created by collapsing the entire bandpass after flagging the strong CO(3–2) line. While the polarized Stokes Q and U data exhibit barely any line emission beyond the CO(3–2) line, the Stokes I continuum emission can have some additional line contamination, in particular toward the peak positions of the hot cores, such as NGC 7538IRS1, W3(H₂O), or AFGL2591. However, past analyses of hot cores typically found line contamination even toward the peak positions of less than 10–20% (e.g., Beuther et al. 2017), and it is even lower for the more extended filamentary structures. Hence, line contamination does not significantly affect our results. The typical rms in Stokes I, Q, and U continuum images is a few mJy beam−¹, a bit higher for Stokes I because of larger side-lobes during the cleaning of the strong submillimeter Stokes I continuum emission. All rms values are listed in Table 2.

3 Results

While the density structure and magnetic field analysis relies on two very different datasets, we start the presentation of results on the larger scales with the single-dish dust continuum observations in Sect. 3.1, and then move into the more central regions studied with the SMA observations in Sect. 3.2.

3.1 Intensity and inferred density profiles

The main goal of the single-dish 1.2 mm continuum mapping of the CORE sample is to derive intensity profiles, and based on these profiles to derive density profiles of the parental gas clumps of the CORE regions. Typical mean densities for such regions based on millimeter continuum data are in the regime of 10⁵ cm⁻³ (e.g., Beuther et al. 2002). To obtain reliable source structures from single-dish continuum mapping, comparatively large maps are needed to account for the inherent spatial filtering during continuum mapping (see also Sect. 2.1). In such scanning map approaches (e.g., Adam et al. 2018), the correlated noise between KIDs has to be considered, and the largest mapped scales are regarded as emission-free and are used to subtract the background (e.g., Sect. 4.4 in Perotto et al. 2020). Therefore, large-scale emission is filtered out in the final map. This approach can also affect source structures. Hence, the larger the maps, the more reliable the measured source structures (e.g., Motte et al. 1998; Beuther et al. 2002; Hatchell & van der Tak 2003; Perotto et al. 2020; Rigby et al. 2021). As outlined in Sect. 2.1, the mapping was conducted with map sizes of ~9.8′, and the final maps are presented in Figs. 1, A.1, and A.2. The magenta circles in these figures outline the primary beam size of the corresponding SMA observations of 36″ discussed in Sect. 3.2. The primary beam size of the original NOEMA 1.3 mm CORE observations was 22″ (Beuther et al. 2018a).

These larger maps outline the general environment of the regions and reveal quite a diversity. While some of the star-forming regions are rather isolated (e.g., IRAS 23151+5912, Fig. 1, or G100.3779-3.578, Fig. A.1), others are part of large-scale structures (e.g., W3IRS4, G75.78+0.34, Fig. A.1, or the NGC7538 complex, Fig. A.2). While these large-scale structures are interesting in themselves and deserve a separate analysis, here we focus only on the intensity and density structures of the main central regions indicated by the magenta ellipses in Figs. 1, A.1, and A.2. Nevertheless, all these 1.2mm continuum maps are provided in electronic form via CDS (Centre de Données astronomiques de Strasbourg) for further analysis.

For our intensity profile analysis, we largely follow the approach described in Beuther et al. (2002) for a sample of 69 high-mass protostellar objects (HMPOs). The two-dimensional intensity distributions are a convolution of the intrinsic source intensity with the 12″ beam, and it is also affected by the scanning mapping technique (see discussion above). Regarding the convolution with the beam, theoretical work and simulations have shown that the beam convolution does not significantly affect scales larger than the beam (e.g., Adams 1991; Motte & André 2001). Structures smaller than the beam size cannot be spatially resolved, but the data still give the integrated continuum flux densities of the inner regions.

Since most central regions resemble to first order spherical structures (Figs. 1, A.1, and A.2), we derive the radial intensity profiles in circular annulli of 3" steps starting at the continuum peak position. Figures 2, A.3, and A.4 present the derived intensity profiles. As already seen in the data of the 69 regions presented in Beuther et al. (2002), outside the spatial resolution limit of ~12″, the profiles typically exhibit broken power-law slopes, steepening to the outside. This steepening to the outside is most likely an observational artifact of the spatial filtering caused by the scanning maps conducted with NIKA2, and the subsequent background subtraction during the calibration. The steepening typically occurs at ~30″ from the center. Following Beuther et al. (2002), we fit the intensity profiles with an inner power-law slope $I\propto r^{-m_i}$ out to a radius of 30″ and an outer power-law slope $I\propto r^{-m_o}$ from 30" to 54″. Since the integrated flux toward the center is finite, the inner power law $I\propto r^{-m_i}$ has to break somewhere, and we model it with an inner flat region. For the combination of m_i and the inner breakpoint we create models, smooth them with a Gaussian beam of 12″, and fit m_i and inner breakpoint simultaneously to the data. That inner flat but unresolved region typically corresponds roughly to the area where fragmentation is observed in the higher-resolution interferometric NOEMA and SMA data. Hence, with the single-dish continuum observations, we trace the density structure of the paternal gas clump and its potential implication for the fragmentation observed at smaller spatial scales with NOEMA and the SMA. In total, our model fitting contains three free parameters m_i m_o, and the inner breakpoint. In practice, the outer profile m_o is directly fitted to the obtained profiles since the beam at those scales is negligible.

As outlined in Sect. 2.1, by conducting tests with artificially created sources, convolving them with the 12″ beam, and injecting typical noise parameters, we explored how much the fitting deviates from the model profiles. For test inputs of I ∝ r^−m with m equal to 2, 1.5, and 1.0, the fitting resulted in corresponding values of m_i of 2.1, 1.7, and 1.3, respectively. The fits to the artificially created profiles also steepen at larger radii, similar to our observed data and confirming the spatial filtering of the scanning mapping approach with NIKA2. While the fitted m_i of the steeper profile is barely affected by the correlated noise and spatial filtering, flatter profiles have increased deviations. This is reasonable since for flatter initial profiles, more larger-scale flux is filtered out in the observations, which steepens the fitted profiles. Based on these tests, we estimate the uncertainties of ∆m_i to ~0.3 (similar results were also found in the past; e.g., Motte & André 2001). Since the outer profiles suffer more from such correlated noise and filtering effects, and the inner breakpoint is unresolved, in the following analysis only the inner power law m_i is considered. All resulting fits are presented in Figs. 2, A.3, and A.4. Furthermore, Table 2 lists the fitted values for the inner power-law profiles m_i. The range of fitted inner intensity profiles m_i between 0.5 and 1.3 is slightly narrower than for the larger HMPO sample of 69 regions in Beuther et al. (2002) where m_i values were found in the range 0.4–2.1. Two regions are in common between our CORE sample and the Beuther et al. (2002) study (IRAS 23151 and IRAS 23033), and while we derive here values of m_i of 1.1 for both sources, Beuther et al. (2002) found m_i of 1.3 for the same regions. Considering that these observations were made with two entirely different continuum instruments decades apart, the similarity of the profiles, consistent within the error budget, is reassuring. Furthermore, only four regions or 20% of the sample exhibit extremely flat intensity profiles with m_i below 1.0. The remaining 80% of the sample cluster in the narrow range of intensity power-law profiles between 1.0 and 1.3.

Figure 3 compares the fitted intensity profiles to the luminosities of the regions (Table 1). While the data may indicate a weak tendency that the steeper profiles could be related to more luminous regions, considering the error bars and the scatter, no reliable correlation between these two quantities can be discerned. This is consistent with the recent compilation of density indices from the literature by Gómez et al. (2021), who also did not find significant differences between samples of low-and high-mass star-forming regions. We also checked whether the intensity profile slopes may be related to the distances of the sources, but no correlation between these parameters can be identified.

A more important question is how the fitted intensity profiles relate to the fragmentation properties of the regions. Although the fragmentation in the CORE interferometer data is observed over a primary beam size of 22″ and the single-dish intensity profiles cover larger radii of ~60″, the difference between barely fragmented to highly fragmented regions may show signatures in the large-scale intensity distributions. In the high-spatial-resolution NOEMA data, low-fragmentation regions typically show emission concentrated in the center of the field, whereas highly fragmented regions are well distributed over the size of the primary beam (Beuther et al. 2018a). Hence, although the absolute number of fragments depends on the size of the primary beam, the sensitivity and the spatial resolution (all similar for the entire CORE sample, Beuther et al. 2018b), the observed variations between barely fragmented to highly fragmented regions is real. To investigate the relation between large-scale intensity distribution and number of cores per region, Fig. 4 shows the power-law indices m_i against the number of cores identified in Beuther et al. (2018a) within each region. While the regions cover a range of distances (Table 1), we compare the intensity structures to the number of fragments within the same regions. Therefore, the distance is not critical for this direct comparison.

No clear trend is visible here, and it is not possible to identify any relation of the single-dish intensity profile with the number of cores. Although the data point with 20 cores and an intensity index m_i of 1.3 appears to be an outlier in Fig. 4, we point out that this location is occupied twice. Both highly fragmented regions with 20 identified cores (IRAS 20178 and G100.3779) exhibit the same steep intensity distribution. With two similar regions, an outlier "rejection" argument seems less likely. This is interesting also when comparing their large-scale distributions in Fig. A.1. While IRAS 21078 is part of a more extended region with filamentary large-scale structures, G100.3779 appears more isolated. However, there may also be more extended low-intensity structures in G100.3779 that could be below our sensitivity limits. Looking at the whole sample, the large-scale intensity distributions seem to be uncorrelated with the number of observed cores or fragments. We return to this point in Sect. 4.

The millimeter continuum emission profiles can also be used to estimate the underlying density profiles. Here, we concentrate only on the inner profiles m_i. Following Adams (1991), Motte & André (2001), and Beuther et al. (2002), the millimeter intensity profile is related to the density profile as $$m=-1+p+Qq+{ϵ}_f,$$(1)

with the intensity profile power-law index m, the density profile n ∝ r^−p, and the temperature profile Τ ∝ r^−q; Q is a temperature and frequency dependent correction factor that is ~1.2 at 30 Κ and 1.2 mm wavelength (Adams 1991; Beuther et al. 2002). The de-projection term ∈_f takes into account that the relation between intensity, density, and temperature profile was originally derived for infinite power-law distributions, but our regions typically have finite sizes (Yun & Clemens 1991). However, ∈_f is rather small, and Motte & Andro (2001) estimated that ∈_f should typically be around 0.1. Regarding the temperature distribution Τ ∝ r^−q, for centrally heated regions like those in our CORE sample, radiative transfer calculations typically find power-law indices q around 0.4 (e.g., Emerson 1988; van der Tak et al. 2000b). Furthermore, Gieser et al. (2021) measured the temperature distributions in the CORE sample from H₂CO and CH₃CN emission line data, and they found an average temperature power-law index q = 0.4 ± 0.1. Hence, we use this q value for the estimates of the density distributions based on the single-dish 1.2 mm data. Using the discussed parameters in Eq. (1), the power-law density index ρ can be approximated as $$p\approx m+1-1.2\times 0.4-0.1\approx m+\mathrm{0.4.}$$(2)

The corresponding estimated inner density power-law indices Pi for all regions are also listed in Table 2. With Gaussian error propagation, the error ∆p_i can be approximated as $\text{Δ}{p}_i\approx \sqrt{\text{Δ}{m}_i^2+\text{Δ}{q}^2}$. With ∆m_i ≈ 0.3 and ∆q ≈ 0.1 we obtain ∆p_i ≈ 0.32 ≈ 0.3. Similarly to the intensity profiles, 80% of the regions exhibit large-scale density distributions with a narrow power-law exponent range between 1.4 and 1.7, clustering around 1.5. Only four regions, or 20% of the sample, have flatter profiles closer to power-law indices of 1.0.

We now compare the density power-law distributions derived on large scales from the single-dish data with the corresponding small-scale density structures estimated from the NOEMA data in Gieser et al. (2021). They estimate the density profiles by fitting the interferometric visibilities directly in the uv-plane with a power-law slope α. This α relates to the power-law slopes p and q of the density and temperature distributions as α = p + q − 3 (e.g., Adams 1991; Looney et al. 2003; Gieser et al. 2021). For a consistent comparison with our data, we derive the small-scale density power-law slope p from their α using again a temperature power-law slope of q = 0.4 (Table 2). We approximate the same Δp ≈ 0.3 for these interferometrically derived small-scale density profiles.

While the large-scale density profiles cluster around a power-law slope of 1.5, the small-scale density power-law slopes have a mean value of 2.0 (Gieser et al. 2021). Although the error bars are large on these measurements, on average there seems to be a steepening of the density profiles from large to small spatial scales. Figure 5 presents the direct comparison of the large-and small-scale density structures. While for many regions the difference between large-scale and small-scale profiles is rather small, interestingly, the two flattest large-scale profiles p_i below 1.0 both exhibit comparably steep small-scale profiles above 2.0. Figure 5 also gives the impression of increasing small-scale p with decreasing large-scale p_i. While this trend is only tentative considering the given uncertainties in the density profiles, we note that Gieser et al. (2023) found in a recent analysis of ALMA data toward high-mass star-forming regions at different evolutionary stages a similar flattening of the density profiles from small core to large clump spatial scales. The density profiles are further discussed in Sect. 4.

In addition to the intensity and density profiles, we can also use the 1.2 mm continuum data to estimate column densities and mean densities toward the peak position. This will be important for the magnetic field analysis following below. Table 2 lists the peak flux densities S_peak derived toward the 20 regions. For those regions that also encompass an ultracompact HII region, the S_peak values are corrected for free-free emission based on the studies listed in footnote (c) of Table 2. Following Hildebrand (1983) or Schuller et al. (2009), assuming optically thin dust emission, a mean dust temperature of 25 K, a gas-to-dust mass ratio of 150 (Draine 2003), and a dust opacity κ ≈ 1.0 cm²g⁻¹ at 1.2 mm (Ossenkopf & Henning 1994), we can estimate the mean column densities in the regions where the polarized dust emission is also detected (see following section). The derived column densities range between ~4 × 10²² and almost 10²⁴ cm⁻², as listed in Table 2. As a next step, assuming a spherical structure of the parental star-forming clump, we can use these column densities to estimate the mean densities by dividing the peak column densities by the corresponding linear extent of the 12″ beam. These mean densities are listed in Table 2, and range between ~0.5 × 10⁵ and several times 10⁶ cm⁻³.

Fig. 1 NIKA2 1.2 mm dust continuum images toward three example regions of the CORE sample. The color scale shows the flux densities, and the contours are always from 3σ to 12σ in 3σ steps. The 12″ beam and a 1 pc scale bar are shown at the bottom left of each panel. The magenta circles outline the ~36″ primary beam size of the corresponding SMA observations. The remaining maps are shown in Figs. A.1 and A.2.

Fig. 2 Radial intensity profiles derived for the main 1.2mm dust continuum sources in the CORE sample. The power-law slopes of the inner and outer fits (mi and m0), and the inner breakpoint in arcsecond are labeled in each panel (see main text for fitting details). The beam size of 12″ corresponds to radial distances of 6″. The remaining fits are presented in Figs. A.3 and A.4.

Fig. 3 Intensity power-law index mi plotted against the luminosity of the region (Table 1). A representative error bar for mi is shown at the bottom right.

Fig. 4 Intensity power-law index mi plotted against the number of cores (Table 1) from Beuther et al. (2018a). The point for 20 cores and a slope of 1.3 consists of two sources. A representative error bar for mi is shown at the bottom right.

Fig. 5 Comparison of density distributions measured on large scale with the IRAM 30 m telescope (pi) with the small-scale density power-law index p derived from the NOEMA data by Gieser et al. (2021). The representative error bars are shown at the bottom left.

Fig. 6 Example SMA polarization maps. The color-scale presents the Stokes I total intensity data. The contours show the same data starting at the 4σ level and continue in 8σ steps. The green constant-length line segments present the magnetic field orientation (polarization angles rotated by 90 deg) derived from the linearly polarized continuum data above the 2σ level (independent of the polarization fraction). The synthesized beam and a linear scale bar are shown at the bottom of each panel. The corresponding maps for the remaining sample are presented in Figs. B.1 and B.2.

3.2 Polarization and magnetic field properties

3.2.1 Alignment of magnetic field and submillimeter continuum emission

Polarized submillimeter continuum emission is detected toward the entire CORE sample. Figures 6, B.1, and B.2 show the polarized emission angles rotated by 90 deg, outlining the plane-of-the-sky magnetic field structure, overlaid on the Stokes I 875 μm continuum emission. In many regions we are able to map the magnetic field structure toward extended, often filamentary structures. Visual analysis indicates that the magnetic field structure is often aligned with filamentary structures leading toward the central source.

A way to quantify the relative orientation between the magnetic field and the gas column density structure is the histogram of relative orientations (HRO), introduced for astrophysical magnetic field studies by Soler et al. (2013; for comparable methods, see also Li et al. 2013 or Liu et al. 2022). The HRO method measures the relative orientation between the plane-of-the-sky magnetic field component and the corresponding column density structure using its gradient. In our case we measure the relative orientations of the magnetic field vectors against the gradient vectors of the Stokes I submillimeter continuum emission, which traces the optically thin dust emission and hence gas column density (e.g., Hildebrand 1983). Figure 7 presents the corresponding results for two examples where the Y-axis quantifies the relative orientation between the magnetic field and continuum intensity gradient (Z_x, also known as the projected Rayleigh statistic V) and the X-axis shows the Stokes I 875 μm continuum emission. Positive Z_x values correspond to preferentially parallel orientation[s] and negative Z_x to largely perpendicular orientations between the magnetic field and gas column density. The null hypothesis in this approach implies a uniform distribution. Values of Z_X of around |1.64| and |2.57| correspond to rejections of the null hypothesis with probabilities of 5 and 0.5%, and Z_X ≈ |2.87| corresponds roughly to a 3σ limit (e.g., Soler et al. 2022; Batschelet 1972). Assuming the same temperature, increasing Stokes I intensities correspond to increasing column densities.

Although the visual inspection of the data indicates that on larger scales of filamentary structures, the magnetic field appears aligned with some of the filaments (e.g., AFGL2591 or G75.78 in Fig. 6), the HRO analysis is less straightforward. While for IRAS 23151 negative Z_X values at low Stokes I intensities are consistent with the magnetic field vectors roughly perpendicular to the densest inner regions (see also Fig. 6), even that is not highly significant. In most other cases, the Z_X values lie around 0. While this HRO analysis does not exclude alignment or misalignment, as indicated by the visual analysis, the number of polarization vectors is typically too low for a proper statistical analysis, as was done for example on larger scales with the Planck data by Soler et al. (2013) or Jow et al. (2018), among others.

Nevertheless, the visual analysis above can be interpreted in a framework where gravity dominates the dynamics of the gas flows, and the gas is channeled along filamentary structures toward the main gravitational centers where the most massive protostars are forming. Magneto-hydrodynamic (MHD) simulations result in comparably aligned structures (e.g., Klassen et al. 2017; Gómez et al. 2018). At even higher spatial resolution, similar results have been found by Beuther et al. (2020), Sanhueza et al. (2021), and Cortés et al. (2021).

Fig. 7 Histogram of oriented gradient (HRO) plots for two example regions, IRAS 23151 (top) and G75.78 (bottom). The X-axis shows the Stokes I intensity, and the Y-axis the projected Rayleigh statistic Zx (also known as V). Positive Zx values correspond to largely parallel orientation and negative Zx to largely perpendicular orientations between magnetic field and gas column density. The dashed lines at ~12.871 correspond to 3σ significance in circular statistics.

3.2.2 Polarization fraction and relation to Stokes I

A different quantity to analyze is the polarization fraction, which is the fraction between linearly polarized and Stokes I intensities. Based on the typical flux calibration uncertainty of ~10% and potential line contamination toward the hot core peak positions discussed in Sect. 2.2, we estimate the uncertainty on the polarization fraction to ~20%. Figure B.3 presents the maps of polarization fractions for the entire CORE sample, and in Fig. 8 we show a histogram of the polarization fractions for all regions. More than half of the observed polarization fractions are below 4%, and 88% are below 10% polarization fraction.

Over the last decades these polarization holes, which may indicate a decrease in fractional polarized intensity toward the Stokes I peak intensities, have been discussed (e.g., Dotson 1996; Fiege & Pudritz 2000; Henning et al. 2001; Matthews & Wilson 2002; Wolf et al. 2003; Girart et al. 2006; Tang et al. 2009; Liu et al. 2013; Hull et al. 2014; Soam et al. 2018; Fernández-López et al. 2021). Physical reasons for measured lower polarization fractions at increased column densities range from unresolved magnetic field structures that may be smoothed out by too low angular resolution to less efficient radiative torque alignment of the dust grains at high densities and high optical depth (e.g., Lazarian & Hoang 2007; Hull et al. 2014; Soam et al. 2018; Girart et al. 2018). Inspecting the CORE sample, we see various structures (Figs. 6, B.1, B.2, and B.3). While we do detect polarized emission toward many of the continuum peak positions (e.g., IRAS 23033, G75.78, and NGC 7538IRS1), toward other regions we find barely any polarized emission toward the 875 μπι peak positions (e.g., IRAS 23151, NGC 7538IRS9, and NGC 7538S; see Fig. B.3). Quantifying this behavior further, Fig. 9 shows the polarization fraction plotted against the Stokes I total intensity at 875 μm wavelength. We clearly see that the polarization fraction decreases with increasing Stokes I intensities. Even for sources where the polarized emission is detected toward the peak position, the fractional polarized intensities still decrease with increasing Stokes I intensities (Fig. B.3). Fitting a power law between the fractional polarized intensities (pol_frac) and the Stokes I intensities (S_I), we find a relation as ${\text{pol}}_{\text{frac}}\propto S_I^{-0.62}.$. This slope is similar to the value of −0.72 recently reported for the high-mass region G5.89-0.39 (Fernández-López et al. 2021).

Fig. 8 Histogram of polarization fractions for the whole sample.

Fig. 9 Polarization fraction plotted against the Stokes I total intensity. Uncertainties for Stokes I are around the calibration uncertainty of ~10%. Uncertainties on the polarization fraction are ~20%. The red line shows a power-law fit to the data with the polarization fraction (polfrac) proportional to the Stokes I intensity $\left(S_I\right):{\mathrm{pol}}_{\text{frac }}\propto S_I^{-0.62}.$.

Fig. 10 Dispersion angles vs. number of cores. The black stars correspond to the directly measured σψ, whereas the red stars correspond to the b(0) value measured with the structure function analysis. The approximate error bars are shown at the top right in black and red respectively for the σψ and b(0) estimated magnetic field strengths.

3.2.3 Magnetic field and fragmentation

The most important thing we wanted to investigate with the polarisation data is how the magnetic field relates to the fragmentation of these high-mass star-forming regions. Figure 10 presents the dispersion angles of polarized emission (i.e., a proxy of the magnetic field) plotted against the number of cores or fragments from Beuther et al. (2018a). Following Appendix D in Planck Collaboration Int. XXXV (2016), we estimate the dispersion of polarization angles σ_ψ directly from the Stokes Q and U images via $${\sigma }_{\psi }=\sqrt{\langle {(\text{Δ}\psi )}^2\rangle },$$(3)

with $$\text{Δ}\psi =0.5\times \arctan \left(\frac{Q\langle U\rangle -\langle Q\rangle U}{Q\langle Q\rangle -\langle U\rangle U}\right),$$(4)

where 〈…〉 denotes the average values in the selected areas of the maps. To estimate reliable dispersion angles, we impose a threshold of the polarized emission to be detected in at least nine pixels, corresponding roughly to one beam size. This way, we can estimate the dispersion angle σ_ψ for 16 regions (Table 3).

To better account for the fact that the entire angle dispersion distribution may not be attributed to MHD waves and turbulence, Houde et al. (2009) and Hildebrand et al. (2009) developed a different approach to estimate the approximate angle dispersion associated with the magnetic field via the second-order structure function analysis S₂(l), where I is a length scale. This structure function S₂(l) measures the differences between polarization angles depending on the distances I between them. The structure function analysis allows the estimation of the angle dispersion, while reducing the potential contributions of large-scale turbulence (Hildebrand et al. 2009). In practice, the angle dispersion σ_ψ is estimated by fitting the square of the structure function with a polynomial of second order $S_2^2(l)=b(l)+a_2{l}^2$ above length scales corresponding to the beam size (see Fig. 11 for an example). The intercept of b(0) is then an alternative measure for the dispersion angle σ_ψ. We are able to derive reasonable structure functions for 13 regions. The results are listed in Table 3 and are shown as red markers in Fig. 10. In many cases the measured b(0) is below the directly inferred dispersion angles σ_ψ from Eq. (3). However, some regions also exhibit b(0) larger than σ_ψ. These are typically at high angle dispersion values, close to ~52 deg, which corresponds to random distributions. In that regime measuring reliable dispersion values is clearly a difficult undertaking.

Even without estimating quantitative magnetic field strengths, the measurement of the dispersion angles σ_ψ already allows a qualitative estimate of the magnetic field strengths. A larger value of σ_ψ indicates a smaller magnetic field strength (see Eq. (5) below). Using this qualitative estimate, we plot the dispersion angle σ_ψ against the number of cores in Fig. 10. This figure shows an interesting feature in the way that the measurements cluster in the top left half of the plot with no entries in the bottom right of the figure. In other words, we indeed see large dispersion angles σ_ψ dominating, but not strictly correlated with, the number of cores. However, the emptiness of the bottom right part of Fig. 10 implies that a larger number of cores is not seen with low dispersion angles σ_ψ. Formulated with respect to the magnetic field, a large number of fragments are not found for large magnetic field values. This indicates that the magnetic field can indeed prevent the parental gas clump from fragmenting, consistent with simulations (e.g., Commerçon et al. 2011, 2022; Peters et al. 2011; Myers et al. 2013, 2014).

Quantifying the magnetic field strength and following Davis (1951) and Chandrasekhar & Fermi (1953), in the Davis-Chandrasekhar-Fermi method (DCF) the angle dispersion σ_ψ is inversely related to the magnetic field in the plane of the sky. If the magnetic field is closely coupled ("frozen") to the gas, and the dispersion of the local magnetic field orientation angles is caused by transverse and incompressible Alfven waves, the magnetic field in the plane of the sky can be estimated via $$B_{\text{pos}}^{\text{DCF}}=\sqrt{4\pi \rho }\frac{{\sigma }_v}{{\sigma }_{\psi }},$$(5)

with ρ the gas density and σ_v the ID velocity dispersion of the gas. The accuracy of the DCF method has been discussed frequently (e.g., Ostriker et al. 2001; Heitsch et al. 2001; Liu et al. 2021, 2022; Skalidis & Tassis 2021; Skalidis et al. 2021), and different correction factors to Eq. (5), typically around 0.5, have been proposed (e.g., Ostriker et al. 2001; Liu et al. 2022). Recently, Skalidis & Tassis (2021) proposed a method that depends on the square root of the angle dispersion $\sqrt{{\sigma }_{\psi }}$ and takes into account compressible motions in the interstellar medium. In the following, we estimate the magnetic field strength following this approach (Skalidis & Tassis 2021): $$B_{\text{pos}}=\sqrt{2\pi \rho }\frac{{\sigma }_{\upsilon }}{\sqrt{{\sigma }_{\psi }}}.$$(6)

While the central densities estimated at the highest resolution (~0.3−0.4″) of the CORE data range between 10⁶ and 10⁸ cm⁻³ (Beuther et al. 2018a), the polarized emission of the SMA data at lower angular resolution (~3″, Table 2) is observed more in the environmental structures at lower densities. Therefore, for the density ρ we use the mean densities derived from the single-dish 1.2 mm data at 12″ resolution in Sect. 3.1 (Table 2), roughly encompassing the areas of detected polarized emission in the SMA data. Furthermore, we use a mean molecular weight μ = 2.8m_p with m_p the proton mass. The ID velocity dispersion is estimated from the single-dish IRAM 30m C¹⁸O(2−l) data (Beuther et al. 2018a; Mottram et al. 2020) over the beam size of ~11″, typically encompassing large parts of the detected polarized emission (Figs. 6, B.1, and B.2). The full width at half maximum ∆v is given in Table 3 (typically a few km s⁻¹), and the corresponding ID velocity dispersion is ${\sigma }_v=\text{Δ}v/\sqrt{8\ln (2)}$.

The estimated magnetic field strengths Β_pos(σ_ψ) and B_vos(b), using either the direct angle dispersion σ_ψ or the structure function intercept b(0), vary roughly between ~0.2 and ~4.5 mG (Fig. 12 left and right panels, respectively, and Table 3). The mean error shown in Fig. 12 is estimated from Gaussian error propagation with a relative error on the number density of a factor of 2 (which could easily be higher, even up to a factor 10; e.g., Liu et al. 2021), the measured line width error of ~0.1 km s⁻¹, and the dispersion angle error of ~3 and ~17 deg from the direct σ_ψ and structure function intercept b(0), respectively. In addition, it should be kept in mind that the absolute error of the Davis-Chandrasekhar-Fermi method could be even larger, as discussed in Planck Collaboration Int. XXXV (2016) or Liu et al. (2021). One caveat is that it is difficult to determine whether the angle dispersion is only caused by magneto-hydrodynamic (MHD) waves and turbulence. Furthermore, the angle dispersion is an average measurement within the beam along the line of sight, potentially reducing the measured projected dispersion again. Therefore, the absolute derived magnetic field strength should be considered as an order-of-magnitude estimate, but does not capture the entire complexity of the magnetic field structure. However, since the assumptions are the same for our entire sample, the relative differences between the sources are considered more reliable. We discuss possible dependences between the magnetic field and fragmentation, as well as further insights about the magnetic field in Sect. 4.2.

Fig. 11 Example structure function S2(l) for CepA. The black line shows the data, and the orange line presents the fit performed for length scales above the beam size that is marked by the two vertical lines (major and minor beam axis). The orange dashed line shows the intercept b(0) and the blue dashed line shows the dispersion angles derived directly from the data, as discussed above. The dashed line at ~52 deg shows the extreme for a random distribution.

Fig. 12 Magnetic field Bpos vs. number of cores. The left panel uses the directly measured σψ, whereas the right panel uses the b(0) value measured with the structure function analysis. Approximate mean errors are shown at the right of each panel.

4 Discussion

4.1 Density structures and fragmentation

Although only tentative, we find a trend of flatter-density power-law slopes on large scales compared to those derived with NOEMA on smaller spatial scales (Fig. 5). Interestingly, Gieser et al. (2023) found a similar trend in recent ALMA observations of high-mass star-forming regions in different evolutionary stages. While power-law indices of around 2 on the small scales observed by Gieser et al. (2021, 2022, 2023) correspond well to classical collapse solutions (e.g., Larson 1969; Shu 1977; Stahler & Palla 2005; Bhandare et al. 2020), the flattening on larger scales to first order does not fit into that picture. However, high-mass star-forming regions are typically parts of larger cloud structures, as also evident in our 1.2 mm IRAM 30 m continuum maps (Figs. 1, A.1, and A.2). Therefore, the fitted structures of the star-forming regions cannot resemble infinite power-law profiles, but they may flatten out when they merge with the environmental molecular clouds. Such a merging of collapsing star-forming clumps with the environmental cloud could then result in a flattening of the observable density structures on large spatial scales. In addition to this, as shown in analytical and numerical studies, even starting with initially flatter density distributions, during the collapse process, the density structures should approach p = 2 (e.g., Naranjo-Romero et al. 2015; Gómez et al. 2021). Assuming that the flatter large-scale profiles found in our study resemble more closely the initial conditions of the region, the steeper interferometrically derived inner density profiles would then correspond to the collapse-induced density structure.

Regarding the finding that we do not see a correlation between the large-scale intensity distribution (and hence density distribution; see Sect. 3.1) and the number of fragments (see also Palau et al. 2014), Girichidis et al. (2011) argued that the initial density profile of a region should be related to the fragmentation of a region in a way that flatter density profiles would result in more fragments, whereas steeper density profiles would result in less fragmentation. Considering the millimeter continuum intensity distribution as a proxy of the density distribution, our data do not allow such conclusions. However, it should be kept in mind that theoretical studies (e.g., Girichidis et al. 2011) consider the density distribution of the initial conditions of the parental gas clump, whereas we observe the gas clump at a later time with already ongoing star formation and fragmentation. Furthermore, in the simulations by Girichidis et al. (2011), the density distribution is the only parameter that varies. If other parameters (e.g., the magnetic field) vary as well, less clear trends in the observational data can be expected. Following the chemical analysis by Gieser et al. (2021), the mean approximate ages of our regions are roughly 6 × 10⁴ yr. Hence, we are not observing the initial intensity distribution, but one that has already evolved over several 10⁴ yr. Nevertheless, if the main collapse were largely progressing in an inside-out fashion, the larger-scale density distribution would not be affected that much during the early evolution. Independent of that, while our data do not allow us to draw a conclusion about a potential relation between the initial density structure and fragmentation properties, the number of fragments apparently does not directly correlate with the corresponding currently observed large-scale intensity and density structures.

To investigate the previously suggested relations between column density and density versus the level of fragmentation (e.g., Lombardi et al. 2013; Palau et al. 2014, 2015), we checked whether our derived column densities and mean densities (Table 2) could be correlated with the number of cores in the respective regions, but we did not find any such correlation.

This is consistent with the findings of Palau et al. (2013), where the mean densities were derived following the same approach as used in this work.

A future step to address a potential relation between the density distribution versus fragmentation properties will be to conduct similar millimeter dust continuum observations toward younger regions (e.g., infrared dark clouds). These should resemble the initial conditions better, and hence allow a clearer association of the initial density structure with the early fragmentation of high-mass star-forming regions.

4.2 Magnetic field structures and fragmentation

While the angle dispersion versus number of cores plotted in Fig. 10 suggests that there may be a relation between the polarization angle dispersion with the number of fragments, after conversion to magnetic field strength in Fig. 12 this trend is not obvious. For a low number of fragments (below 10) there is a broad spread of angle dispersion and magnetic field strength; instead, the one region with 20 cores, for which we can also estimate a magnetic field strength⁴, exhibits an intermediate magnetic field strength around ~1 mG. Hence, the strength of the magnetic field itself does not allow us to predict the number of fragments.

Additional important magnetic field parameters can be estimated from the data, in particular the Alfvenic velocity, the ratio of turbulent to magnetic energy and the mass-to-flux ratio. The ID Alfvenic velocity σ_A can be estimated from the magnetic field strength Β and the density ρ via ${\sigma }_A=B/\sqrt{4\pi \rho }$. Using a mean Β(σ_ψ) ≈ 1.0 mG (Table 3 and Fig. 12) and a mean density of ρ ≈ 0.9 × 10⁶ cm⁻³, the ID Alfvenic velocity is σ_A ≈ 1.37 km s⁻¹. With the uncertainties in Β and p, the σ_A uncertainty is roughly a factor of 2. For comparison, based on the mean line width of ≈3.2 km s⁻¹ (Table 3, based on single-dish C¹⁸O(2−1) data at 11″ resolution, Gieser et al. 2021⁵), the mean ID velocity dispersion is σ_υ ≈ 1.4 km s⁻¹. Considering the given uncertainties, these values are surprisingly close and indicate that Alfvenic and turbulent velocities are roughly on the same order of magnitude.

In a similar vein, one can estimate the ratio of turbulent to magnetic energy β ~ 3(σ_v/σ_A)² (e.g., Girart et al. 2009). Using the ID Alfvenic and turbulent velocity dispersion values discussed above, one finds a β ≈ 3.1. With the velocities being squared for the estimate of β, the uncertainties on β are around a factor of 4. Nevertheless, even such rough estimates indicate the similar importance of magnetic and turbulent energies in these regions, with a tendency for slightly higher turbulent energies. This is different to the turbulent-to-magnetic field energy ratios found in younger still starless regions where β was estimated << 1 (e.g., Beuther et al. 2018b), indicative of a greater importance of the magnetic field at earlier evolutionary stages.

In addition to sample estimates, we can also derive β for each region and plot that against the number of fragments, which is shown in Fig. 13. While there is a scatter in β between ~1 and ~5.4 for fragmented regions with fewer than ten cores, the two regions with more than ten cores, for which we again also have a β measurement, both have turbulent-to-magnetic energy ratios β of around 4. Such a trend could be interpreted as more kinetic energy favoring more fragments. However, this trend remains inconclusive with the given uncertainties.

A different way to assess the stability of the regions against collapse is estimating the mass-to-flux ratio $M/{\Phi }_B~7.6\times {10}^{-24}\frac{N_{\text{H}2}}{B}$ with the mass M, the magnetic flux Φ_Β, the column density N_H2, and the magnetic field strength B. The latter two are in units of cm⁻² and mG, respectively (e.g., Crutcher 1999; Troland & Crutcher 2008). The mass-to-flux ratio Μ/Φ_Β is then given in units of the critical mass-to-flux ratio (Μ/Φ_Β)_crit where values greater or smaller than 1 correspond to collapsing or more stable configurations, respectively. For the magnetic field B, we use Β_pos(σ_ψ) (Table 3). For the column density N_H2, we use the mean column densities derived from the 1.2 mm single-dish data in Sect. 3.1 (Table 2). To obtain an estimate for the uncertainties of Μ/Φ_Β, we apply Gaussian error propagation assuming a mean error on Β of ~0.63 mG (Fig. 12) and an uncertainty on N_H2 of the same order as the measured mean column density. The estimated mass-to-flux ratios Μ/Φ_Β are listed in Table 3. All mass-to-flux ratios are greater than 1, ranging between ~2 and ~7. Figure 14 plots Μ/Φ_Β versus the number of cores within each region. Even considering the given uncertainties on Μ/Φ_Β, for all regions the mass-to-flux ratio is ≥2, consistent with all regions already collapsing and actively forming stars. Recently, Palau et al. (2021) proposed a tentative positive correlation between the mass-to-flux ratio and the fragmentation level. For comparison, we do not see any clear trend between Μ/Φ_Β and the fragmentation of the regions, and hence can neither confirm nor reject that tentative relation found before.

We also checked whether there could be any correlation between the magnetic field strength and the density structure measured on large and small spatial scales with IRAM 30 m and NOEMA, respectively (Fig. 5). However, no correlation could be identified, which indicates that the clump and core density structures of the regions are largely independent of the magnetic field strength.

Fig. 13 Ratio of turbulent to magnetic energy β vs. number of cores. Approximate uncertainties on β are around a factor 4.

Fig. 14 Mass-to-flux ratio vs. number of cores. For clarity in the logarithmic plotting, the red error bars are only shown one-sided to higher values.

5 Conclusions and summary

With the goal of identifying and characterizing the main physical processes that determine the fragmentation properties of high-mass star-forming regions, we investigated the large-scale density distributions and the magnetic field properties of the 20 high-mass star-forming regions constituting the CORE sample by means of IRAM 30 m continuum and SMA polarization observations. Although the data allow additional possible investigations, for example large-scale environmental studies with the IRAM 30 m data and chemical and physical analyses of the corresponding SMA spectral line data, here we concentrate on the Stokes I and linearly polarized continuum emission to characterize the density and magnetic field structure and set that into context with the fragmentation properties previously derived by the high-angular-resolution PdBI/NOEMA data of the CORE program (Beuther et al. 2018a).

The IRAM 30 m 1.2 mm dust continuum data allow us to infer the density structure of the regions. While the measured intensity profiles on larger scales (>35″) exhibit a steepening, this is an observational artifact caused by the continuum mapping approach filtering out emission beyond those scales. However, analysis of the inner intensity profiles reveals that they should be accurate within ∆m_i ~ 0.3, resulting in similar uncertainties for the inner density structure profiles p_i. To the first order, we find no correlation between the estimated density structures and the number of fragments (see also Palau et al. 2014). This may partially be due to the evolutionary stage: the regions are typically in the high-mass protostellar object stage and may not reflect the initial conditions anymore. However, it remains to be investigated whether the density structures do change significantly on the given timescales of a few times ~10⁴ yr, or whether they stay relatively constant in the framework of an inside-out collapse. In the latter cases, the given observations could still resemble the initial conditions comparably well.

Additionally, we compared the large-scale density structures with the small-scale density distributions previously derived from the interferometer data by Gieser et al. (2021). Interestingly, we find that the large-scale density distributions are typically flatter (power-law slopes around ~1.5) than the small-scale density structure (typical power-law slopes around ~2.0). A similar steepening of the density structure toward smaller scales was recently also identified for other regions with ALMA observations by Gieser et al. (2023). Possible reasons for this behavior are that the star-forming regions are still embedded in larger cloud structures, causing a density profile flattening on larger scales. Furthermore, as discussed for simulations by Gómez et al. (2021), among others, even starting with initially flatter density distributions, the collapse of the regions steepens the profiles typically approaching slopes of around 2.0. Hence, the measured large-scale flatter intensity profiles may resemble the initial structures more closely, where the smaller-scale profiles are caused by the collapse motions.

Regarding magnetic field structure, we detect a significant polarization signal in 16 out of 20 regions. This is one of the largest high-resolution magnetic field studies today (for comparable SMA samples, see Zhang et al. 2014 or Palau et al. 2021). With ALMA, larger samples with high angular resolution are now starting to be studied as well (e.g., Sanhueza et al. 2021; Cortés et al. 2021). Within our CORE sample we find magnetic field structures appearing aligned with filaments leading toward the dense central cores. However, because of the still relatively low number of individual polarization measurements within individual regions, a statistical quantification of this behavior (e.g., with the histogram of oriented gradients,) is still limited.

More than half of the observed polarization fractions are below 4%; 88% are below 10%. While some regions exhibit polarization holes toward the Stokes I peak positions, as previously discussed in the literature, others show polarized emission also toward the intensity peaks. Nevertheless, the polarized intensities are inversely related to the Stokes I intensities and roughly follow a power-law slope ${\text{pol}}_{\text{frac}}\propto S_I^{-0.62}$.

Estimating the magnetic field strength via a modified DCF method (Skalidis & Tassis 2021), we find a range of magnetic field strengths between ~0.2 and ~4.5 mG. While the original dispersion of polarization angles may indicate a weak trend between angle dispersion and number of fragments, the derived magnetic field strength does not exhibit a clear trend with the number of fragments. The mass-to-flux ratio varies between ~2 and ~7, consistent with all regions being collapsing and actively forming stars. Comparing the estimated 1D turbulent and Alfvenic velocity dispersion as well as the turbulent to magnetic energy ratio, we find that the turbulent and magnetic energies in this sample are of similar importance.

Coming back to our original question of whether the level of fragmentation depends on the density and/or magnetic field structure, for this CORE sample we find no clear correlation between these parameters. This result may allow the tentative conclusion that apparently high-mass star-forming regions fragment almost independently of their initial density distribution or magnetic field properties.

However, the presented density and magnetic field analysis of the CORE sample is obviously limited by several parameters, for example evolutionary stage or spatial scales. Hence, this points to several different study approaches for forthcoming projects. Regarding the density structure analysis, similar studies of younger regions reflecting potentially better the initial conditions (e.g., infrared dark clouds) are needed to better study the connection of the initial conditions with the magnetic field properties. With respect to the magnetic field studies, the SMA data presented here investigate intermediate spatial scales for a source sample covering only a comparatively narrow range in evolutionary stages. Extending the sample size to younger and older regions will allow an even better investigation of evolutionary changes. Furthermore, it will be important to connect these intermediate spatial scales to larger cloud scales and to smaller core and disk scales. For the larger cloud scale, the forthcoming NIKA2POL polarization capability at the IRAM 30m will allow us to directly connect the SMA magnetic field measurements with the larger cloud-scale magnetic field structure. Regarding smaller spatial scales, higher-resolution SMA observations as well as potential polarization at NOEMA are possible in the northern hemisphere. While the CORE sample is only accessible from the northern hemisphere, in the south ALMA allows much higher angular resolution studies of the innermost regions. All these opportunities outline the exciting prospects that are on the horizon.

Acknowledgements

The authors are grateful to the staff at Pico Veleta and Submillimeter Array observatories for their support of these observations. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. This work is based on observations carried out under project number 143-19 with the IRAM 30m telescope. We particularly like to thank Stefano Barta and Robert Zylka with helping in the NIKA2 data reduction process. R.K. acknowledges financial support via the Heisenberg Research Grant funded by the German Research Foundation (DFG) under grant no. KU 2849/9. A.P. is grateful to Gilberto Gómez and Enrique Vázquez-Semadeni for very insightful discussions. A.P. acknowledges financial support from the UNAM-PAPIIT IG100223 grant, the Sistema Nacional de Investigadores of CONAHCyT, and from the CONAHCyT project number 86372 of the ‘Ciencia de Frontera 2019’ program, entitled ‘Citlalcóatl: a multiscale study at the new frontier of the formation and early evolution of stars and planetary systems’, México. A.S.M. acknowledges support from the RyC2021-032892-I and PID2020-117710GB-I00 grants funded by MCIN/AEI/10.13039/501100011033 and by the European Union ‘NextGenerationEU’/PRTR, as well as the program Unidad de Excelencia María de Maeztu CEX2020-001058-M. J.D.S. acknowledges funding by the European Research Council via the ERC Synergy Grant “ECOGAL – Understanding our Galactic ecosystem: From the disk of the Milky Way to the formation sites of stars and planets” (project ID 855130).

Appendix A IRAM 30m data

Fig. A.1 NIKA2 1.2 mm dust continuum images toward the CORE sample. The color-scale shows the flux densities, and the contours are from 3σ to 12σ in 3σ steps. The 12″ beam and a 1 pc scale-bar are shown at the bottom left of each panel. The magenta circles outline the ~ 36″ primary beam size of the corresponding SMA observations.

Fig. A.2 NIKA2 1.2 mm dust continuum images toward the CORE sample. The color-scale shows the flux densities, and the contours are from 3σ to 12σ in 3σ steps. The 12″ beam and a 1 pc scale-bar are shown at the bottom left of each panel. The magenta circles outline the ~ 36″ primary beam size of the corresponding SMA observations.

Fig. A.3 Radial intensity profiles derived for the main 1.2 mm dust continuum sources in the CORE sample. The power-law slopes of the inner and outer fits (mi and mo) as well as the inner breakpoint in arcsecond are labeled in each panel (see main text for fitting details).

Fig. A.4 Radial intensity profiles derived for the main 1.2 mm dust continuum sources in the CORE sample. The power-law slopes of the inner and outer fits (mi and mo) as well as the inner breakpoint in arcsecond are labeled in each panel (see main text for fitting details).

Appendix B SMA data

Fig. B.1 Remaining SMA polarization maps. The color-scale presents the Stokes I total intensity data. The contours show the same data starting at the 4σ level and continue in 8σ steps. The green constant-length line segments present the magnetic field orientation (polarization angles rotated by 90 deg) derived from the linearly polarized continuum data above the 2σ level (independent of the polarization fraction). The synthesized beam and a linear scale bar are shown at the bottom of each panel.

Fig. B.2 Remaining SMA polarization maps continued. The color-scale presents the Stokes I total intensity data. The contours show the same data starting at the 4σ level and continue in 8σ steps. The green constant-length line segments present the magnetic field orientation (polarization angles rotated by 90 deg) derived from the linearly polarized continuum data above the 2σ level (independent of the polarization fraction). The synthesized beam and a linear scale bar are shown at the bottom of each panel.

Fig. B.3 Polarization fraction maps for all regions are shown in color-scale. The color bar is the same for all and is shown next to the top right panel. The contours show the Stokes I total intensity data starting at the 4σ level and continue in 8σ steps.

References

All Tables

All Figures

Fig. 1 NIKA2 1.2 mm dust continuum images toward three example regions of the CORE sample. The color scale shows the flux densities, and the contours are always from 3σ to 12σ in 3σ steps. The 12″ beam and a 1 pc scale bar are shown at the bottom left of each panel. The magenta circles outline the ~36″ primary beam size of the corresponding SMA observations. The remaining maps are shown in Figs. A.1 and A.2.

In the text

	Fig. 2 Radial intensity profiles derived for the main 1.2mm dust continuum sources in the CORE sample. The power-law slopes of the inner and outer fits (mi and m0), and the inner breakpoint in arcsecond are labeled in each panel (see main text for fitting details). The beam size of 12″ corresponds to radial distances of 6″. The remaining fits are presented in Figs. A.3 and A.4.
In the text

	Fig. 3 Intensity power-law index mi plotted against the luminosity of the region (Table 1). A representative error bar for mi is shown at the bottom right.
In the text

	Fig. 4 Intensity power-law index mi plotted against the number of cores (Table 1) from Beuther et al. (2018a). The point for 20 cores and a slope of 1.3 consists of two sources. A representative error bar for mi is shown at the bottom right.
In the text

	Fig. 5 Comparison of density distributions measured on large scale with the IRAM 30 m telescope (pi) with the small-scale density power-law index p derived from the NOEMA data by Gieser et al. (2021). The representative error bars are shown at the bottom left.
In the text

Fig. 6 Example SMA polarization maps. The color-scale presents the Stokes I total intensity data. The contours show the same data starting at the 4σ level and continue in 8σ steps. The green constant-length line segments present the magnetic field orientation (polarization angles rotated by 90 deg) derived from the linearly polarized continuum data above the 2σ level (independent of the polarization fraction). The synthesized beam and a linear scale bar are shown at the bottom of each panel. The corresponding maps for the remaining sample are presented in Figs. B.1 and B.2.

In the text

Fig. 7 Histogram of oriented gradient (HRO) plots for two example regions, IRAS 23151 (top) and G75.78 (bottom). The X-axis shows the Stokes I intensity, and the Y-axis the projected Rayleigh statistic Zx (also known as V). Positive Zx values correspond to largely parallel orientation and negative Zx to largely perpendicular orientations between magnetic field and gas column density. The dashed lines at ~12.871 correspond to 3σ significance in circular statistics.

In the text

	Fig. 8 Histogram of polarization fractions for the whole sample.
In the text

	Fig. 9 Polarization fraction plotted against the Stokes I total intensity. Uncertainties for Stokes I are around the calibration uncertainty of ~10%. Uncertainties on the polarization fraction are ~20%. The red line shows a power-law fit to the data with the polarization fraction (polfrac) proportional to the Stokes I intensity $\left(S_I\right):{\mathrm{pol}}_{\text{frac }}\propto S_I^{-0.62}.$.
In the text

	Fig. 10 Dispersion angles vs. number of cores. The black stars correspond to the directly measured σψ, whereas the red stars correspond to the b(0) value measured with the structure function analysis. The approximate error bars are shown at the top right in black and red respectively for the σψ and b(0) estimated magnetic field strengths.
In the text

Fig. 11 Example structure function S2(l) for CepA. The black line shows the data, and the orange line presents the fit performed for length scales above the beam size that is marked by the two vertical lines (major and minor beam axis). The orange dashed line shows the intercept b(0) and the blue dashed line shows the dispersion angles derived directly from the data, as discussed above. The dashed line at ~52 deg shows the extreme for a random distribution.

In the text

	Fig. 12 Magnetic field Bpos vs. number of cores. The left panel uses the directly measured σψ, whereas the right panel uses the b(0) value measured with the structure function analysis. Approximate mean errors are shown at the right of each panel.
In the text

	Fig. 13 Ratio of turbulent to magnetic energy β vs. number of cores. Approximate uncertainties on β are around a factor 4.
In the text

	Fig. 14 Mass-to-flux ratio vs. number of cores. For clarity in the logarithmic plotting, the red error bars are only shown one-sided to higher values.
In the text

	Fig. A.1 NIKA2 1.2 mm dust continuum images toward the CORE sample. The color-scale shows the flux densities, and the contours are from 3σ to 12σ in 3σ steps. The 12″ beam and a 1 pc scale-bar are shown at the bottom left of each panel. The magenta circles outline the ~ 36″ primary beam size of the corresponding SMA observations.
In the text

	Fig. A.2 NIKA2 1.2 mm dust continuum images toward the CORE sample. The color-scale shows the flux densities, and the contours are from 3σ to 12σ in 3σ steps. The 12″ beam and a 1 pc scale-bar are shown at the bottom left of each panel. The magenta circles outline the ~ 36″ primary beam size of the corresponding SMA observations.
In the text

	Fig. A.3 Radial intensity profiles derived for the main 1.2 mm dust continuum sources in the CORE sample. The power-law slopes of the inner and outer fits (mi and mo) as well as the inner breakpoint in arcsecond are labeled in each panel (see main text for fitting details).
In the text

	Fig. A.4 Radial intensity profiles derived for the main 1.2 mm dust continuum sources in the CORE sample. The power-law slopes of the inner and outer fits (mi and mo) as well as the inner breakpoint in arcsecond are labeled in each panel (see main text for fitting details).
In the text

Fig. B.1 Remaining SMA polarization maps. The color-scale presents the Stokes I total intensity data. The contours show the same data starting at the 4σ level and continue in 8σ steps. The green constant-length line segments present the magnetic field orientation (polarization angles rotated by 90 deg) derived from the linearly polarized continuum data above the 2σ level (independent of the polarization fraction). The synthesized beam and a linear scale bar are shown at the bottom of each panel.

In the text

Fig. B.2 Remaining SMA polarization maps continued. The color-scale presents the Stokes I total intensity data. The contours show the same data starting at the 4σ level and continue in 8σ steps. The green constant-length line segments present the magnetic field orientation (polarization angles rotated by 90 deg) derived from the linearly polarized continuum data above the 2σ level (independent of the polarization fraction). The synthesized beam and a linear scale bar are shown at the bottom of each panel.

In the text

	Fig. B.3 Polarization fraction maps for all regions are shown in color-scale. The color bar is the same for all and is shown next to the top right panel. The contours show the Stokes I total intensity data starting at the 4σ level and continue in 8σ steps.
In the text