© ESO 2021

1 Introduction

The ionising radiation of early type (massive) stars generates a photon dominated region (PDR) where it hits a nearby molecular cloud. Molecular clouds are complex entities, having a self-similar structure over a wide range of scales (e.g. Ossenkopf et al. 2007). With increasing spatial resolution molecular clouds appear fragmented in substructures of smaller angular dimensions. Hence, the molecular material is typically very clumpy, implying that a large part of the molecular cloud can be exposed to far-ultraviolet (FUV) radiation (6 < hν < 13.6 eV), which can penetrate much deeper than if the cloud were homogeneous, making the PDR volume vastly larger. The radiative and mechanical (winds) action of massive stars can evaporate existing structures, but it may also trigger the formation of a new generation of stars in the surrounding material, for example through radiation-driven implosion (RDI; see e.g. Bisbas et al. 2011).

The importance of gas fragmentation in star-forming regions has been reinforced by the discovery that the smallest (~ 0.1 pc) and densest structures (≳10⁴−10⁵ cm⁻³), named cores, inside giant molecular clouds follow a well-defined mass distribution (core mass function: CMF) resembling the stellar initial mass function (IMF; see Rathborne et al. 2009; Ikeda & Kitamura 2009; André et al. 2014 and references therein). We note that cores are embedded in molecular clumps, which are larger (0.3−3 pc) and less dense (10³−10⁴ cm⁻³) gas concentrations (see e.g. Bergin & Tafalla 2007). Although the link between the CMF and the IMF is still debated, there seems to be growing consensus that the IMF derives from the CMF through further fragmentation and star formation (taking into account the formation ofbinaries and multiple systems; Goodwin et al. 2008), provided that the star-forming efficiency (SFE) of each core is less than ~30−50% (Alves et al. 2007). On the other hand, Motte et al. (2018) have recently found that the CMF in the mini-starburst galactic region W43-MM1 challenges this scenario of direct proportionality, in that its slope is “markedly shallower” than that of the IMF, for core masses greater than about a solar mass; variations in the CMF are indeed possible in different environments.

If the CMF does in fact evolve into the stellar IMF, then the distribution of stellar masses will be linked to the origin of dense cores, and different CMFs may significantly affect the mass distribution of the final stellar population. The characteristics of the CMF (its slopes as well as the presence of a turn-over and its location in mass) are likely to depend on conditions in the molecular clouds in which the cores are embedded, the properties of which, in turn, depend on the environment in which they were formed. Therefore, it has to be expected that if intense UV radiation has any effect on star formation, it must be by primarily affecting the CMF in the exposed molecular cloud(s).

The region of interaction between young massive stars and a molecular cloud is thus of great interest. Hydrodynamical simulations (e.g. Dale et al. 2007; Bisbas et al. 2011; Arthur et al. 2011; Wall et al. 2020; González-Samaniego & Vazquez-Semadeni 2020) have shown that star formation may be either triggered or halted in clouds exposed to intense FUV radiation. The FUV flux intensity is possibly a key element in determining one of the two scenarios. Feedback from young massive stars is believed to be instrumental in regulating the SFE in galaxies (Krumholz 2014; MacLachlan et al. 2015). Lucas et al. (2017) propose that the removal of molecular cloud envelopes by stellar feedback can significantly decrease a cloud’s SFE. Thermal feedback affects the way in which a molecular cloud fragments and can lead to a top-heavy IMF (Krumholz 2014). A way of observationally assessing the effects of intense FUV radiation fields on molecular gas is therefore to determine whether the CMF near to a PDR significantly differs from that observed in more quiescent regions (i.e. in regions less exposed to intense stellar FUV radiation).

NGC 6357 is a galactic complex of molecular clouds and HII regions that hosts several active star-forming sites. Massi et al. (2015) discussed the issue of distance, concluding that NGC 6357 is located 1.7 kpc away. A more reliable assessment can now be based on Gaia DR2 data (Gaia Collaboration 2018). Xu et al. (2018) retrieved the parallaxes of a large sample of OB stars from the Gaia DR2 catalogue. Those associated with NGC 6357 lie in the range 1.7− 2.0 kpc, with quoted uncertainties of ~10%. In particular, Ramírez-Tannus et al. (2020) derive a distance of 1.77 ± 0.12 kpc and Maíz Apellániz et al. (2020) a distance of ${1690}_{-110}^{+130}$ pc, both using Gaia data. This definitely rules out a distance as large as 2.56 kpc, which was derived by Massey et al. (2001). Thus, we will safely continue using d = 1.7 kpc as in Massi et al. (2015), which is likely to be correct within ~10%.

Figure 1 in Giannetti et al. (2012) shows that the molecular gas in NGC 6357 is arranged in a ring surrounding a large (~ 15 pc × 10 pc) cavity (or a collection of connected smaller cavities), possibly shaped by feedback from one or more star clusters (Wang et al. 2007). In the northern part of this complex is located G353.2+0.9 (cf. Fig. 3 in Massi et al. 2015), which is the most luminous region found in the ring at wavelengths from radio to optical. Just 55″ (0.45 pc at a distance of 1.7 kpc) south of G353.2+0.9 lies the massive, young open star cluster Pismis 24, which contains some of the brightest and bluest stars known (O3.5; Maíz Apellániz et al. 2007). G353.2+0.9 is a blister HII region seen nearly face-on, with an extended ionisation front (IF) or PDR on the side facing the ionising stars and with an elephant trunk seen prominently against the ionised gas. South-east of G353.2+0.9 is G353.1+0.6, another HII region with an IF or PDR, seen edge-on, associated with the cluster AH03J1725-34.4 (Massi et al. 1997, 2015). East of these two is another HII region, G353.2+0.7, associated with cluster B of Massi et al. (2015). G353.2+0.9, G353.1+0.6,and G353.2+0.7 are ideal targets for studying the molecular gas fragmentation expected in FUV-exposed gas, considering both the extent of the PDR around the giant molecular clouds observed and the relatively close distance.

Several tracers of molecular density and temperature were observed with the SEST (Massi et al. 1997; Giannetti et al. 2012) in a 3′ × 3′ area around G353.2+0.9 (total mass ~2000 M_⊙). An analysis of the complex molecular emission profiles led to the identification of at least 14 clumps in the area mapped, indicating that the gas is indeed fragmented. However, the angular resolution (21–54′′) of these single-dish observations is insufficient to resolve some clumps for which there are indications of smaller-scale structures. Near-infrared (J, H, and K) images of G353.2+0.9 obtained with the ESO-NTT telescope (Massi et al. 2015), indicate the presence of several embedded objects coinciding with the brightest part of the IF, which are younger than the stars of Pismis 24. This suggests that star formation in the region has been triggered by the action of the massive stars. Giannetti et al. (2012) also used the ATLASGAL map of the region at 870 μm (Schuller et al. 2009), showing the most prominent clumps seen at millimetre lines. However, although its spatial resolution (~ 18′′, ~ 0.2 pc) is enough toresolve gas at the typical scale of young small embedded clusters, this map has a poor sensitivity (100 mJy beam⁻¹), and most of the clumps found at millimetre lines went undetected.

NGC 6357/G353.2+0.9 represents a rather extreme environment, with 2-3 O3.5f stars and many later-type O stars at ~ 1 pc from the cloud. We can expect the intense FUV field to affect the nearby gas, for example by rapidly photoevaporating the smallest cores and thus producing a flatter CMF. Other parts of the complex are more quiescent and have no clear signs of star formation, and therefore no intense FUV-radiation field is expected there. Because of this contrast, NGC 6357 may be one of the best locations to study the influence of FUV radiation on star formation. In this paper we identify dense cores from dust-continuum data in an area of about 0.5 square degrees in NGC 6357. Combining these data with Herschel Hi-GAL (Molinari et al. 2010) images we construct the spectral energy distribution (SED) for each core. Fitting a greybody to those SEDs yields temperatures and masses for all cores, from which a mass function is derived. By dividing the field into a sub-field where a more intense FUV field is expected and a sub-field less exposed to FUV radiation, we carry out a comparison between the mass functions from the two sub-fields.

The layout of the paper is the following. In Sect. 2 we describe our observations and data reduction. Our results are presented in Sect. 3 and discussed in Sect. 4. Finally our conclusions are listed in Sect. 5.

2 Observations and data reduction

We used the Submillimetre Common-User Bolometer Array 2 (SCUBA-2: Holland et al. 2013) at the James Clerk Maxwell Telescope (JCMT) to observe 12 partially overlapping fields, simultaneously at 450 and 850 μm: on 13 July 2012 (fields 1–6; project m12an004) and on 6 and 23 April and 22 May 2015 (fields 7–12; project m15ai093). The observations were carried out using a so-called CV Daisy pattern¹ ; each daisy-field observation took 30 min. Each field has a radius of about 8 arcmin, and together they cover the entire ring-like structure of NGC 6357, including three star clusters, as illustrated in Fig. 1. The images were reduced using SMURF and associated packages contained in STARLINK², following the Cookbook version SC/21.1³ (Chapin et al. 2013). After reduction, the fields were mosaicked together. The resulting sensitivity varies with location, and except at field edges is ≲ 4− 15 mJy beam⁻¹ (at 850 μm), and 60–400 mJy beam⁻¹ (at 450 μm). We used the flux conversion factors (FCFs; see Dempsey et al. 2013) listed in Table 1. These were derived from measurements of CRL2688 taken close in time to each series of observations. The Starlink SDF files were then converted into FITS format for further analysis. The (effective) beam size of SCUBA-2 is 9′′ ± 1′′ at 450 μm and 14′′ ± 1′′ at 850 μm (i.e. ~ 0.07 pc and ~ 0.12 pc, respectively). To complement these observations and construct SEDs, we used the 70 μm and 160 μm images from the Herschel/Hi-GAL survey (Molinari et al. 2010), with beam sizes of 5′′ and 11.′′ 5, respectively.

Correction for line-contamination

The 850 μm emission is potentially contaminated by ¹²CO(3–2) line emission at 345.8 GHz. The continuum band is very wide (ca. 39 GHz, see Naylor et al. 2014) and lines have to be bright and broad in order to significantly affect the continuum emission. Several studies (e.g. Sadavoy et al. 2013) show that the contamination is typically less than ~ 70 mJy beam⁻¹, but may reach ~ 150 mJy beam⁻¹ (10− 30 times the rms level in our map), and would thus affect especially the lower-level emission and hence the lower-mass cores, which are important to define the low-mass turnover point in the CMF. As there are numerous locations with star formation in NGC 6357, outflow emission, and therefore contamination at relevant levels, is to be expected.

On the other hand, ¹²CO(6–5) falls at the edge of the 450 μm band. No significant contamination is expected (also considering the high noise level of the 450 μm continuum data).

We used the HARP receiver at the JCMT (see Buckle et al. 2009) to map the ¹² CO(3–2) emission in the ~1900′′ × 1700′′ region covered by the SCUBA-2 observations on 2, 7, and 20 May and 9 and 18 July 2014 (project m14au032). Observations were made in ‘basket-weave’ raster-scanning mode. HARP has 16 SIS receptors, 14 of which were operational at the time of the observations. The reference position (17:25:36.2 −33:55:17) was checked using frequency switching and showed only the telluric CO line; rms = 0.3 K ($T_{\text{A}}^{*}$) at 0.3 km s⁻¹ resolution. We used the ACSIS backend with a bandwidth of 1000 MHz, which provided 2048 channels with a spectral resolution of 0.488 MHz (0.423 km s⁻¹ at 345.796 GHz). The data were reduced using the ORAC-DR pipeline (Jenness et al. 2015); the maps were then converted into CLASS format for further analysis. The main-beam efficiency of HARP is 0.64, and all line intensities were converted into a T_MB scale. The average rms in the CO map is 2.0 K, on the T_MB scale. The beam size of HARP at 345 GHz is 14′′.

We followed the procedure outlined by Parsons et al. (2018) to correct continuum data for line contamination, using the HARP ¹² CO(3–2) integrated intensity map as a negative fake map to be combined with the SCUBA-2 data in the map maker (see also Drabek et al. 2012). For this we made use of a script kindly provided by Dr. Parsons.

Fig. 1 Herschel HiGAL 160 μm map (colour scale) with fields observed in 2012 (1–6; green) and 2015 (7–12; blue). The location and approximate size of star clusters are marked in yellow (Massi et al. 2015). The linear size of the region shown is approximately 22 × 17 pc2.

3 Results

The ¹²CO(3–2) map used to correct the 850 μm continuum emission is shown in Fig. 2.

Following the procedure outlined in Sect. 2, we obtained CO-corrected 850 μm maps for each of the 12 observed fields. The individual fields at each wavelength were then combined into the single map shown in Fig. 3 using the task ‘wcsmosaic’.

3.1 Core identification

The 850 μm map has thebest combination of sensitivity and angular resolution (see Sect. 2), and we used this to identify and extract the cores. Because we wanted to identify cores in this complex, we needed to isolate the most compact structures, which are more likely to collapse and form a single star or a close multiple system, depending on the level of fragmentation, and are therefore more closely connected with the IMF. To facilitate the identification, we first removed the extended 850 μm emission in which the cores are embedded. We therefore performed a multi-scale decomposition of the image, using a median filter (see Belloche et al. 2011). We used seven levels of decomposition on S∕N maps, using only those with scales ≲ 64′′ in the image reconstruction, thus filtering out any structure more extended than that. This procedure results in the image shown in Fig. 4. We note that the spatial resolution at 850 μm is ~ 0.1 pc, so cores are at best barely resolved. Nevertheless, hereafter we will designate the compact structures found as ‘cores’.

We first constructed a data cube using pyfits, inserting empty two-dimensional planes before and after the actual continuum data, as suggested in Mookerjea et al. (2004), and modifying the header accordingly; on this data cube we ran Gaussclumps (Stutzki & Guesten 1990) as implemented in GREG⁴. The algorithm decomposition of the filtered S/N map was then visually inspected to identify and remove spurious sources, confirmed as such by a comparison with the Herschel Hi-GAL maps. We also excluded cores outside of the region observed in CO with HARP, and at the noisy edges of the continuum mosaic. A minimum size threshold was set at 13′′ (~ beam at 850 μm). The final catalogue contains 1221 cores.

The distributions of fluxes and sizes of the cores in this catalogue are shown in Figs. 5 and 6, respectively. In the latter figure, we have also made a distinction between cores in regions with and without intense FUV-flux, as defined in Sect. 3.2.

Fig. 2 12CO(3–2) integrated emission ∫ TMB dV; − 20 to + 20 km s−1). Contours values are 80(80)470 Kkm s−1 (low(step)high).

3.2 ‘Active’ and ‘quiescent’ regions

To investigate how the enhanced FUV radiation fields from newly born massive stars affect the molecular environment and the ongoing star formation activity, and to avoid effects due to possible dust opacity differences in intense and moderate FUV environments (see Sect. 3.5), we divided the mapped region into two sub-regions, hereafter the ‘active’ and the ‘quiescent’ regions. The active region includes the parts of the cloud complex more exposed to a high FUV-photon flux, such as G353.2+0.9. On the other hand, the quiescent region includes molecular clouds farther from high-mass stars. Thus, we basically separated the complex into its eastern and western halves (see Fig. 7). The eastern half contains the cluster Pismis 24 (along with the HII region G353.2+0.9), the cluster AH03J1725–34.4 (along with the HII region G353.1+0.6), and the cluster B (along with the HII region G353.2+0.7). This represents the active region. The western half is bound to be less affected by FUV radiation and has been labelled as quiescent region. We note that we are only interested in finding out if a global signature of different levels in the FUV is apparent in the core populations in the two areas that may be further investigated in future works. Thus we do not need very accurate selection criteria to define the two regions at this stage.

3.3 Spectral energy distributions

Because we are interested in the CMF, we need to transform the fluxes of the cores into masses. Making use of the Hi-GAL images we constructed the core SEDs by complementing our SCUBA-2 850 μm and 450 μm flux densities with those from the 160 and 70 μm frames taken at the central pixels of the cores identified in the 850 μm map. The Herschel- and the 450 μm maps were first convolved tothe resolution of the 850 μm map and had their extended emission subtracted in the same way as done for the 850 μm frame. We fitted a greybody to the SEDs for those cores with at least three data points and a 450 μm flux density larger than the CO-corrected 850 μm one, varying the dust temperature between 5 and 80 K in steps of 0.5 K, and the surface density in the range 10⁻³−10^1.5g cm⁻², in steps of 0.1 in logarithmic scale. The dust opacity κ_ν was assumed to be 0.9 cm² g⁻¹ at 1.3 mm (Ossenkopf & Henning 1994), and the index β equal to 1.8, consistent with the dust models used. This yields κ_ν = 1.85 cm² g⁻¹ at 870 μm, the value adopted for the ATLASGAL survey (Schuller et al. 2009). A more detailed analysis on the adopted dust opacity is reported in Sect. 3.5 We were able to derive a dust temperature for 686 cores; their locations are shown in Fig. 7. The active region contains 411 cores and the quiescent region 275, out of the 686 considered. The distribution of their temperatures is displayed in Fig. 8. Some examples of SED-fits are shown in the appendix.

An important difference between the cores populating the two regions can be seen in Fig. 8. The temperature of cores in the quiescent region peaks at ~ 25 K, whereas that of cores in the active region peaks at ~35 K. Cores in the quiescent region are thus significantly colder than those in the active region. Formal errors on the temperature from the greybody fits are nearly always less than 1 K. This clearly demonstrates the effect of FUV photons on the gas. This result is contrary to that found by Bobotsis & Fich (2019), who do not find evidence of HII regions affecting the surrounding cores in their study of 38 complexes of ionised and molecular gas. However, their target regions are mostly farther away than 2 kpc and do not contain very massive ionising stars.

3.4 Core masses

The dust temperatures (T_d) of the best-fit models are used to calculate the mass of the cores, from the integrated 850 μm flux densities from Gaussclumps by means of: $$M=\gamma D^2{S}_{\nu }/[{\kappa }_{\nu }{B}_{\nu }(T_{\text{d}})],$$(1)

where D is the distance to the complex, B_ν(T_d) is the Planck function evaluated at temperature T_d, S_ν is the total flux density and κ_ν is the dust opacity, both at 850 μm; γ is the gas-to-dust ratio (assumed to be 100). In Table 2 we list the properties of all the bona fide cores. The CMF in the NGC 6357 region was constructed through a kernel density estimate, using an Epanechnikov kernel with a bandwidth estimated using the direct plug-in method (Feigelson & Babu 2013). Figure 9 shows the results of this procedure (green dashed line). The slope of the high-mass end of the CMF (M > 5 M_⊙; the completeness limit – see Sect. 3.6) is − 2.1 ± 0.2, derived by directly fitting a power-law distribution to the data points using PyMC (Patil et al. 2010).

The total gas mass of dense cores is ~1000 and ~1500 M_⊙ in the quiescent and active regions, respectively. Assuming a molecular mass of ~ 1.8 × 10⁵ M_⊙ for the whole NGC 6357 complex (Cappa et al. 2011, re-scaled to 1.7 kpc), that implies a fraction of gas in dense cores of ~ 1.4%. Such low values are typical in the Galaxy (Csengeri et al. 2016), reflecting a low SFE (e.g. Bobotsis & Fich 2019).

Fig. 3 SCUBA-2 observations of the dust emission. Left: 850 μm mosaic; the emission has been corrected for the integrated 12CO(3–2) line emission. Right: 450 μm mosaic. Units are mJy beam−1.

Fig. 4 850 μm map after filtering out emission from extended structures as explained in the text. Units are mJy beam−1. The black lines indicate the region mapped in 12CO(3-2) (see Fig. 2).

3.5 Dust opacity uncertainties

As mentioned above (Sect. 3.3), we adopted a dust opacity of 0.9 cm² g⁻¹ at 1.3 mm and β = 1.8, leading to 1.93 cm² g⁻¹ at 850 μm, based on the theoretical work of Ossenkopf & Henning (1994). Their tabulated values are obtained for dense protostellar cores (~ 10⁵−10⁸ cm⁻³) that have remained stable for at least 10⁵ yr. We have also assumed that the dust opacity is constant throughout the sample of cores, so it is essential that we should check whether this is a good approximation before deriving global properties such as mass distributions and mass-size relationships. Although Ossenkopf & Henning (1994) note that uncertainties on the tabulated values should not add up to more than a factor of 2, one has to take into account that most of the cores in NGC 6357 have been exposed to an intense FUV field, a scenario that was cautioned but not modelled by Ossenkopf & Henning (1994). In addition, from Tables 2 and 3 in their paper it is clear that dust opacitydepends on the gas density and that its core-to-core variations will be much smaller if ice mantles have grown on the grains. Unfortunately, one can reasonably expect that part of these ice mantles either have been removed or did not grow at all in the presence of a FUV field. We note that systematic errors of even a factor of 2–3 in the dust opacity are acceptable, in that they do not affect global properties of the cores such as the shape of their mass distribution. This is not true for large core-to-core variations in the dust opacity, which in principle may heavily affect the CMF slope.

To test the effects of dust opacity uncertainties in the region, one can then assume the worst case scenario of grains with no ice mantles at all. We derived the dust opacity at 850 μm as a function of density from Tables 2 and 3 of Ossenkopf & Henning (1994) by a linear fit in the log-log space. For the 452 cores with a deconvolved size of at least half the telescope beam and using the masses derived in Sect. 3.4, we derive a range of average densities ~ 5 × 10⁴ − 10⁵ cm⁻³. Many cores therefore may have densities less than 10⁵ cm⁻³; the fit allows us to extrapolate down to the lowest values. This may not give the correct dust opacity in those cases, but it should be good enough for our test, taking into account the limited range spanned by the core densities.

Starting from the masses obtained from Eq. (1) with the adopted dust opacity (1.93 cm² g⁻¹ at 850 μm), we computed the density of each of the 452 cores and recalculated their mass using the dust opacity for their density, as given by the fitted function. We then iterated this procedure for each core until convergence was reached. The masses obtained by assuming a constant opacity (of 1.93 cm² g⁻¹) and the ones obtained with a density-dependent dust opacity are compared in Fig. 10a. For most cores, the masses obtained with a density-dependent dust opacity are within a ~ 30% of the masses obtained by adopting a constant opacity. Only three very massive cores exhibit significant differences (with the density-dependent opacity always yielding lower masses).

We performed another test by constructing and comparing the mass distributions, which are shown in Fig. 10b. These are mostly consistent with each other within ~ 1 σ (for a Poisson statistics) and their difference is almost negligible in the log-log space, provided the two most massive bins are excluded. According to the mass-size relationships derived in Sect. 4.2, it is easy to see that the most massive cores are also the less dense, so presumably they are more affected by opacity variations. In fact, a linear fit to the distributions in Fig. 10b (for M ≥ 0.4 M_⊙) yields a slope of − 0.76 ± 0.06 for constant opacity and − 0.71 ± 0.07 for density-dependent opacity, in other words consistent with each other within 1 σ. So the shape of the core mass distribution is only moderately affected by the choice of the dust opacity. A further simple check can be performed to see if the sample masses depend on the size. We derived mass-size relationships by a linear fit in the log-log space, namely: $$\log (M/M_{\odot })=p\times \log (D/\text{arcsec}).$$(2)

We obtained p = 2.35 ± 0.15 in the constant opacity case, and p = 2.49 ± 0.12 in the density-dependent opacity case. Again, the difference is at a ~ 1 σ level.

We note that these results were obtained by using the average core densities. However, we expect that significantly higher densities arereached only in a small fraction of a core’s total mass. Furthermore, considering that the density-dependent dust opacity is computed for the worst case of grains with no ice mantles, the effects of dust opacity uncertainties on the derived masses may be even lower. A constant dust opacity is then a good approximation, yielding core-to-core mass uncertainties of < 30% in most cases.

Fig. 5 Distribution of fluxes for the 1221 cores identified in the 850 μm mosaic.

Fig. 6 Distribution of core sizes for all 1221 cores (uncorrected for beam; black solid line); the red curve indicates the active region and the blue one the more quiescent part.

Fig. 7 All 686 cores with a valid SED fit superimposed on the mosaic at 850 μm. Darker symbols indicate lower temperatures. The white line indicates the dividing line between active (to the east) and quiescent (to the west) regions.

Fig. 8 Temperature distribution for all 686 cores (solid black line), cores in the active region only (411; solid red line), and cores inthe quiescent region only (275; solid blue line).

Fig. 9 Core mass function in NGC 6357 (dashed green curve for all cores) obtained as density distribution of sources after smoothing with a Epanechnikov kernel. The red curve indicates the active region, and the blue curve indicates the more quiescent part. The vertical line marks the 5 M⊙ completeness limit (see Sect. 3.6). The small coloured bars show the masses of the single sources for each sample, following the colour coding of the distributions.

Fig. 10 Core masses and core mass distribution. (a) Core masses derived using a density-dependent opacity (grains with no ice mantles) vs. masses derived using a constant dust opacity of 1.93 cm 2 g−1 (at 850 μm). The full line marks equal masses, whereas the dashed line marks density-dependent opacity masses equal to 70% of the constant density masses. (b) Core mass distributions for masses obtained with a constant dust opacity of 1.93 cm 2 g−1 (at 850 μm; full line) and for masses obtained with a density-dependent opacity (dashed line). The error bars have been computed assuming a Poisson statistics.

Fig. 11 Completeness of the cores in our sample, as a function of mass. Starting from the completeness as a function of flux density for point sources, we have used a representative source size of 22′′ and the observed temperature distribution (Fig. 8) to derive the distribution in mass shown here (see text). The dashed line indicates the mass where the completeness reaches its maximum.

Fig. 12 Continuum emission at 850 μm towards G353.2+0.9, overlaid with the locations of the gas cores found by Giannetti et al. (2012; red circles of size equal to that of a clump) and of the dense cores retrieved in this work (black ellipses). The gas clumps are labelled according to Giannetti et al. (2012). Units are mJy beam−1. Contour levels are 50, 150, and 300 mJy beam−1.

3.6 Completeness

To evaluate the completeness of our data we followed these steps: Firstly, we inserted a set of 20 artificial cores in the CO-corrected 850 μm map with a FWHM of 14′′ (the SCUBA-2 PSF at 850 μm) and a constant peak flux density. Secondly, the images with the synthetic cores were filtered as done for the original data. Thirdly, Gaussclumps was run on this processed image and a catalogue of cores was produced, and finally this catalogue was cross-matched with the position of the artificial cores. If one or more sources were detected within a beam, we compared their fluxes and sizes, which must be within a factor of 2.5 and 3 of the input values, respectively,to decide whether the injected source was detected. This was done for peak flux densities of 5, 10, 20, 40, 60, 80 and 100 mJy beam⁻¹, and the procedure was repeated 100 times for each of these values. In this way we obtained the completeness as a function of peak flux density. A completeness of 90% is reached for 60 mJy beam⁻¹. This corresponds to a point-like source with a mass of M_compl,point ~ 0.9 M_⊙ for a dust temperature of 20 K and M_compl,point ~ 0.5 M_⊙ for a dust temperature of 30 K. One way to estimate the mass completeness limits of extended sources would be by assuming a Gaussian spatial distribution of the emission, and then use the distribution over sizes to see how the completeness of point sources would change; one can expect that $$M_{\text{compl}}~M_{\text{compl},\text{point}}\times (D^2+B^2)/B^2,$$(3)

where B and D are the beam- and deconvolved source sizes, respectively.

We used an alternative method. For each peak flux density value used for the artifical core experiments (i.e. 5, 10, 20, 40, 60, 80 and 100 mJy beam⁻¹) we computed a cumulative mass distribution from Eq. (1) according to the following procedure. Once the peak flux density is fixed, we considered a source size of ~ 22′′, halfway between unresolved sources (~ 13′′) and the largest (not beam-deconvolved) size values of ~ 30′′, and derived the mass distribution obtained from the total flux and the observed distribution of dust temperatures (Fig. 8). Thus, each peak flux density is associated with a range of masses that reflects the empirical temperature distribution and whose completeness was assumed to be that estimated for this peak flux density from the artificial core experiments described earlier. The final mass completeness curve was then computed by averaging the completeness levels derived for each mass value from the adopted set of peak flux densities, and is shown in Fig. 11.

The completeness as a function of mass reaches its maximum around 5 M_⊙. The maximum is not 100% due to confusion effects that prevent the recovery of artificial cores that end up in bright and crowded regions. We note that our estimates of the mass completeness level assume that the efficiency in retrieving a core mostly depends on its peak flux density, irrespective of its size. This appears reasonable. However, one has to consider that a large fraction of the flux from weak extended cores is likely to be lost in the noise. So the retrieved total flux density from weak extended sources is bound to be affected by larger errors.

4 Discussion

4.1 Comparison with previous works

The CMF for all cores was derived as discussed in Sect. 3.4 and shown in Fig. 9. To assess the goodness of both the core retrieval algorithm and the masses obtained, we compare our results to those of previous works. Giannetti et al. (2012) mapped the area towards G353.2+0.9 in several molecular mm-lines and used data from ATLASGAL (Schuller et al. 2009), with a beam of ~ 19.′′2 to find the corresponding dust emission in the same region. They identified 14 clumps from line emission, only three of which with a clear compact counterpart in the continuum emission (their clumps C1, C2, and E). Figure 12 displays the continuum emission at 850 μm, overlaid with the positions of the gas clumps found by Giannetti et al. (2012) and with those of our dust cores.

Clearly, the increase in spatial resolution has resulted in being able to decompose the continuum emission in a much larger number of compact structures. The goodness of the decomposition operated by Gaussclumps can be confirmed by a visual examination of the figure. Clump C1 corresponds to our core ID 1, clump C2 to our core ID 4, and clump E to our cores ID 2 and ID 20 (see Table 2). We note that Giannetti et al. (2012) assume a distance of 2.56 kpc, so their masses have to be scaled to our adopted distance of 1.7 kpc.

As a first test, we degraded our (line-corrected) 850 μm map to the ATLASGAL resolution. Then we performed a rough aperture photometry over the areas of the continuum cores found by Giannetti et al. (2012, which they label 1 to 5; see their Table 3). By assuming the same dust temperature (i.e. 30 K), we retrieve masses that are a factor of 0.5−1.2 of theirs, confirming that the calibration at 850 μm is consistent with that of the ATLASGAL data.

As for the gas clumps A to P of Giannetti et al. (2012), we found that the masses of our cores are a factor of 0.3−0.5 of those of the corresponding gas clumps. The gas clumps are colder than our dense cores and using the gas temperature would drive the masses of our cores into accordance with those derived from the gas. But it should be noted that using the gas temperature for the dust cores 1 to 5 of Giannetti et al. (2012) would increase their masses well above those of the corresponding gas clumps. This indicates that the difference in mass between our cores and the gas clumps is not due to the adopted temperatures, but to the higher spatial resolution and sensitivity of the 850 μm map, which allow us to retrieve denser and more compact structures embedded in the gas. This is also confirmed by the fact that a few gas clumps are only associated with faint compact continuum sources or faint diffuse emission.

Our results can also be compared with those of Russeil et al. (2019). These authors used Herschel/PACS and SPIRE data to select a robust sample of 155 dense cores in NGC 6357 with the algorithm getsources (Men’shchikov et al. 2012). The ranges of size and mass that they find are roughly consistent with ours. However, their mass distribution is peaked around 30–50 M_⊙, with a median mass of 22 M_⊙ and a maximum mass of 386 M_⊙, while our CMF exhibits a maximum mass of 176 M_⊙ and a peak at ~1 M_⊙. Thus we have achieved a much lower mass completeness limit than Russeil et al. (2019), which is likely to be caused by not only the differentalgorithms used to decompose the continuum emission but also our preliminary filtering out the more extended emission. We further performed a closer comparison selecting their three most massive sources (their MDC 1, 2, and 5). As shown in Fig. 13, we have retrieved those sources (our cores ID 13, 7, and 14, respectively), along with a number of fainter compact structures around them. The masses of our corresponding cores are a factor of 4− 5 lower than theirs, which cannot be due to the temperatures used or the difference in the adopted dust mass opacity (their prescriptions would yield ~ 1.24 cm² g⁻¹ at 850 μm). This clearly indicates that the cores in our sample are systematically less massive than theirs. As mentioned above, this difference is likely to arise not only due to the different algorithms used to decompose the emission, but also to our preliminary filtering out lower spatial frequency emission.

4.2 Pre-stellar and protostellar cores

Pre-stellar cores are defined as gravitationally bound cores that do not host protostars yet. Separating protostellar and pre-stellar cores is critical, in that the CMF of ‘protostellar’ cores may be affected by mass erosion due to increased turbulence and outflows, and the development of internal temperature gradients, and is less suitable to compare with stellar IMFs. Furthermore, dust opacity is sensitive to the physical environment and may vary by up to a factor of ~ 2 in different physical environments, as discussed in Sect. 3.5. It is then essential that the selected sample contains a class of objects as homogeneous as possible, with a limited range in density.

Identifying protostellar cores requires probes of protostellar activity, such as compact emission at 70 μm. Thus, we cross-checked our 850 μm data and Herschel Hi-GAL 70 μm maps for spatial coincidences. All sub-mm cores associated with 70 μm emission have then been labelled as protostellar, whereas all remaining cores have been labelled as ‘starless’. It has been shown that 70 μm emission is a sensitive probe of protostellar activity (see e.g. Koenig et al. 2012). The distribution of the peak flux densities of the 70 μm emission associated with the sub-millimetre cores, has a maximum at ~ 4 mJy beam⁻¹. This represents a rough completeness limit. Assuming that the infrared emission comes from a Gaussian spatial distribution ~ 20′′ in size and a FWHM~9′′ for the Herschel beam at 70 μm, this yields a flux of ~24 mJy. Then by using Eq. (2) of Dunham et al. (2008) one can estimate a detection limit for the luminosity of the central protostar L_int ~ 0.3 L_⊙. One problem with this approach is that in the active region, dust is heated from the outside as well, through FUV radiation, thus increasing the core emission at 70 μm possibly mimicking the effects of a central protostar. In fact, we found 129 ‘70 μm-bright’ sources out of 411 cores in the active regions, and only ten 70 μm-bright sources out of 275 cores in the quiescent region. All 70 μm-bright cores in the active regions are associated with one of the two HII regions G353.2+0.9 and G353.1+0.6. This indicates a possible bias due to far-infrared emission from warm dust in the outer parts of cores located inside or near to the HII regions.

We then applied a further selection filter to the 70 μm-dark cores, based on the spatial coincidence of cores and red WISE sources. After retrieving all point sources from the WISE Source Catalog⁵ projected towards NGC 6357, we selected those with colours of Class I sources according to Koenig et al. (2012). All contaminants (extragalactic sources, PAH emission, shock emission) were discarded following the criteria of Koenig et al. (2012). We also included sources detected in only two bands with d log(λF_λ)∕d logλ > −0.3 and sources detected at only 22 μm. In addition, we required photometric errors <0.3 mag in all the relevant bands. Then we cross-correlated this cleaned list of infrared sources and the one of sub-mm cores. All 70 μm-dark cores whose ellipse overlaps the positional uncertainty ellipse of a WISE source were re-classified as protostellar; this concerned 31 cores in the active region and 41 in the quiescent region. In the end, we found 251 starless and 160 protostellar cores in the active region and 224 starless and 51 protostellar cores in the quiescent region. In principle we could also use the WISE sources to verify the nature of the 70 μm-bright cores in the areas of the active region closest to the HII regions and in the PDR. Unfortunately this is impossible, due to the saturation of the WISE bands; therefore in those areas all 70 μm-bright cores remain tentatively classified as protostellar.

In Fig. 14 we show the temperature distributions for the various (sub)samples. To avoid effects due to incompleteness, only cores with M > 5 M_⊙ have been selected. Figure 14a shows the distribution for all cores, and emphasizes what we saw in Fig. 8 for cores of all masses. The temperature distributions for 70 μm-dark and 70 μm-bright plus WISE-associated cores in the active and quiescent regions are shown in Figs. 14b, c. The temperatures in the quiescent region for both classes are similar, peaking between 15–20 K with an average T ~ 19 K. On the other hand, in the active region the distribution of 70 μm-bright cores is clearly shifted to higher temperatures with respect to that of 70 μm-dark cores (from an average T ~ 22 K to an average T ~ 28 K).

The sensitivity limit for protostellar cores associated with WISE sources can be estimated as follows. We will assume here that the central protostellar source has an infrared flux increasing with wavelength, with a spectral index d log(λF_λ)∕d logλ ~ 2. This is typical of embedded Class I sources. As the flux is higher in the WISE band at 22 μm, we need to derive the WISE sensitivity to point sources in this band. The 5σ sensitivity in the 22 μm band quoted in the explanatory supplement⁶ is 5.4 mJy (i.e. the magnitude at 22 μm ([22]) ~ 8). As the flux at 22 μm and the spectral index are known, we can compute the flux that such a source would exhibit in the other WISE bands. Following Kryukova et al. (2012), who relate infared spectra to central bolometric luminosities, F₂₂ ~ 5.4 mJy and spectral index 2 correspond to a protostar of ~ 3.5 L_⊙. However, WISE is much more sensitive in the other bands, thus central sources with the same spectral index can be undetected at 22 μm but detected in the lower wavelength bands, increasing the actual sensitivity. On the other hand, due to saturation in the WISE bands there is in fact no way to discriminate between starless and protostellar cores in the areas towards the two HII regions, that is to say the most critical ones to understand the effects of FUV feedback. Even using the MSX Point Source Catalogue, which has a worse spatial resolution and sensitivity but saturation at higher infrared fluxes, yields only 20 matches near to the HII regions, still indicating sensitivity problems. All cores associated with MSX point sources have temperatures in the range 26–63 K, with only two cores above 5 M_⊙ (T = 26 and 34.5 K). Moreover, the histogram plotting the number of protostellar sources versus [22] displays peaks at [22] ~ 1.5 mag and ~2.5 mag in the active and quiescent regions, respectively, indicating that the actual completeness limit is much worse than inferred from the nominal sensitivity limit.

A rough analysis of the CMF shows that little is changed in the statistical properties of the population of starless cores when excluding those associated with WISE sources. Binning the CMF in intervals of 0.4 in logarithmic masses and fitting the high-mass end (log[M∕M_⊙] > 0.8) with a N ~ M^α function, one obtains α = −1.79 ± 0.09 for 70 μm-dark cores and α = −1.71 ± 0.1 for 70 μm-dark cores cleaned of WISE-associated sources in the quiescent regions. In the active regions, one obtains α = −1.72 ± 0.09 for 70 μm-dark cores and α = −1.62 ± 0.12 for 70 μm-dark cores cleaned of WISE-associated sources. The difference is not statistically significant, probably due to the fact that only a few sources are moved between samples if their mass is above the completeness limit.

The discussion above highlights the difficulties inherent in obtaining a census of real starless cores. Although 70 μm emission is in principle a very sensitive probe of protostellar emission, nevertheless we found many cores associated with red WISE sources and no emission at 70 μm. In addition, we have no way of checking thenature of 70 μm-bright cores around the two HII regions, which are the sources most affected by the FUV feedback. On the other hand, we have shown that the association with red WISE sources is not a sensitive tool in the whole region, because of saturation. As a result, thereis currently no way of unambiguously separating starless and protostellar cores in large areas towards NGC 6357.

For the sake of simplicity, we will refer to the sample of 70 μm-dark cores cleaned of those WISE-associated as the ‘starless’ sample, despite its uncertain degree of contamination from protostellar cores. It is composed of 224 sources in the quiescent region and 251 sources in the active region. On the other hand, there are only 51 candidate protostellar cores in the quiescent part and 160 possible protostellar cores in the active region (31 cores that are associated with red WISE sources outside the HII regions plus 129 70 μm-bright cores).

Not all starless cores are also pre-stellar. Some are just transient, short-lived structures that will dissipate rather than form stars. To check how much our sample of starless cores would also be representative of pre-stellar cores we followed the analysis of André et al. (2014). In Fig. 15, core masses and sizes can be compared to the half-mass loci of Bonnor-Ebert critical mass for T = 30 K and T = 20 K. It is clear that cores with masses larger than the completeness limits derived from the sub-mm emission, cannot be thermally supported against gravity. Even though turbulence and magnetic fields can oppose gravitational collapse, it is unlikely that they will prevent it. So, most of our starless cores, if really starless, would also be pre-stellar in nature, provided that only masses above the completeness limit are taken into account.

Mass and size of all cores in the two regions are plotted in Fig. 15a. We note that only cores with deconvolved sizes greater than half the beam size have been used. A linear fit (in log–log space) to the data according to Eq (2) yields p = 2.3 ± 0.2 both for the active and the quiescent region. This is consistent with log(M∕M_⊙) ~ 2 × log(D∕arcsec), typical of core populations with constant column density. Similarly, a linear fit to the mass versus deconvolved size of starless cores (shown in Fig. 15b) yields p = 2.4 ± 0.3 in the active region and p = 2.3 ± 0.2 in the quiescent region. The mass-size relations obtained from the linear fit to protostellar cores (shown in Fig. 15c) are quite close to those: p = 2.0 ± 0.3 in the active region and p = 2.4 ± 0.4 in the quiescent region. The differences in the power-law indices are clearly not significant. In addition, one has to be aware that the fit results may be biased due to size-dependent incompleteness (see Sect. 3.6). Nevertheless, the mass-size relationships derived in various star-forming regions or from various samples of cores or clumps displays power-law indices in the range 1.7−2.2 and somewhat depending on the algorithm used, as discussed by Massi et al. (2019), hence similar to the value found in NGC 6357.

Fig. 13 Continuum emission at 850 μm towards the fields of cores MDC 1, 2, and 5 (labelled in figure) of Russeil et al. (2019). Units are mJy beam−1. The locations of the cores identified in this work are overlaid (small open squares).

Fig. 14 Distribution of temperatures for cores in the active (full line) and quiescent (dashed line) regions. Only cores with M > 5 M⊙ (i.e. above the mass completeness limit) have been included. (a) All cores. (b) Cores with no emission at 70 μm and no red WISE sources associated. (c) Cores associated with 70 μm emission and/or red WISE sources.

Fig. 15 Massvs. deconvolved size for (a) all cores; (b) 70 μm dark cores cleaned from WISE-associated cores; (c) 70 μm bright cores plus cores associated with WISE sources. Cores with a deconvolved size smaller than half the beam size have been discarded. Small filled squares mark cores in the active region, whereas open green triangles mark cores in the quiescent region. The critical half-mass loci of Bonnor-Ebert spheres are drawn for T = 20 K and T = 30 K (black dash-dotted lines). The dotted red line marks the completeness limit for T = 20 K dust (from Eq. (3)). Linear (log-log) fits are indicated by a solid (cores in the quiescent region) and a dashed (cores in the active region) line.

4.3 Nature of the cores in the PDRs

Although we were not able to reliably separate protostellar and starless cores in the molecular gas adjacent to the HII regions, nevertheless their nature can be assessed if we assume that those cores are embedded in PDRs. As we said, all 70 μm-bright cores in the active region lie near either G353.2+0.9 or G353.1+0.6. Thus, we applied to this sample the models developed by Gorti & Hollenbach (2002). These authors investigated the evolution of FUV-illuminated clumps (we note that they designate cores as ‘clumps’) in PDRs by examining both the case of turbulent clumps in a vacuum and that of pre-existing clumps in a confining interclump medium (ICM). They found that clump evolution in a PDR is driven by η_c0 (the ratio of the initial clump column density to that of the FUV-heated clump surface, ~ 10²¹ cm⁻²), ν (the ratio of sound speed at the surface to that in the innermost clump regions), and the turn-on time of the heating flux (compared to the initial sound crossing time in the clump).

We computed the deconvolved size for each of the 160 70 μm-bright cores in the active region and discarded all those smaller than half the beam size. We then determined volume and column density for each of the remaining 102 sources. These were further divided into two subsamples, based on their location either towards G353.2+0.9 or G353.1+0.6. The FUV field (G₀) impinging on each core, in Habing units of 1.6 × 10⁻³ erg cm⁻² s⁻¹ was derivedfrom Eq. (14) of Giannetti et al. (2012) by scaling to a distance of 1.7 kpc and an inclination of the ionising sources to the core direction of 45°. We adopted a total FUV stellar luminosity of 1.2 × 10⁶ L_⊙ for G353.2+0.9 (Giannetti et al. 2012) and scaled the FUV field derived by Massi et al. (1997) to the estimated star-core distance for G353.1+0.6. We did not account for extinction of the FUV field and shadowing for the outermost cores, we just assumed that the clumpiness of the gas allows the radiation to penetrate deeply in the molecular clouds bordering the HII regions. Furthermore, the relatively high adopted inclination of 45° may partly compensate for extinction.

From G₀ and the average density we derived the core surface temperature using the PDR model of Kaufman et al. (1999). In turn, we derived ν for each core from the surface and far-IR temperatures. It should be noted that ν measures the intensity of the FUV field. The dimensionless column density parameter η, namely the ratio of the current core column density to that of the FUV-heated core surface (~ 10²¹ cm⁻²), was also computed from the mass and the deconvolved size of each core. The results are shown in Fig. 16, adapted from diagnostic diagrams from Gorti & Hollenbach (2002; see e. g. their Figs. 2 and 7).

Gorti & Hollenbach (2002) consider two classes of cores populating a cloud PDR. One class is represented by turbulent cores in a vacuum that form quickly in the PDR in a time less than that needed by the IF to cross the cloud PDR. It should be noted that these authors do not account for gravity in their analytical treatment. Assuming a PDR thickness of ~ 1 pc and an IF speed of ~1 km s⁻¹, the lifetime of these cores would be <1 Myr. According to Gorti & Hollenbach (2002), these cores (if massive enough) are impulsively heated by the FUV field and would rapidly attain a constant column density irrespective of their initial column density. This critical column density only depends on ν, the non-thermal to thermal pressure ratio (α), and the magnetic to thermal pressure ratio (β); dot-dashed lines in Fig. 16 mark the loci of critical density for α = β = 1 and α = β = 2, see Eqs. (6)–(7) of Gorti & Hollenbach (2002).

The second class these authors consider are pre-existing cores. Pre-existing cores much smaller than the thickness of the cloud PDR and embedded in an ICM would enter gradually the PDR and would experience a steadily increasing FUV field (and ν). The solid lines in Fig. 16 mark the evolution of the column density parameter η against the parameter ν, which will steadily increase as the core moves deeper into the PDR (Eqs. (24)–(25) of Gorti & Hollenbach 2002). The core experiences photoevaporation and the fraction of mass lost during the evolution is indicated by dotted lines in the figure. The dashed line marks the gas column density that would be completely photoevaporated given a FUV field ν (Eq. (26) ofGorti & Hollenbach 2002).

The cores in the PDRs around G353.2+0.9 or G353.1+0.6 show the same distribution and are clearly more consistent with pre-existing cores embedded in an ICM than turbulent cores. Although the value G₀ has not been corrected for extinction, we note thatthe highest ν values always derive from the highest values of G₀. Therefore, these cores represent the ones nearest to the ionising sources, which are the least affected by extinction. In the density range spanned by the rest of the sample, the core surface temperature is not much affected by G₀ and even a decrease of an order of magnitude would result in a decrease in temperature by a factor of 2–3. Therefore, log (ν) would decrease at most by ~0.24.

The cores appear to cluster around the critical density for pre-existing gas structures or lie below it in Fig. 16. This suggests that when the increasing ν (i.e. the IF moves closer and closer to the core) causes the eroded gas column density to decrease below the critical density, a core will be quickly photoevaporated. However, this is inconsistent with the cores with the highest column density displaying the highest values of ν. In fact, this would imply that the cores with the highest column density are located in regions where G₀ is the highest, that is to say close to the IFs. But there is no reason why equally high column densities should not be found closer to the bulk molecular gas (exhibiting low values of ν). Then, eitherturbulence and/or gas compression becomes important close to the IFs and new cores are formed or, more likely, lack of spatial resolution artificially merges smaller cores increasing the measured column density (we note that the IF is seen edge-on towards both HII regions).

This analysis suggests that the PDRs preceding the IFs penetrate into the molecular clouds through a population of pre-existing molecular cores embedded in a lower density ICM. Some of these cores are then likely to experience an RDI with a subsequent gravitational collapse (see Bisbas et al. 2011, and references therein). This is consistent with the elephant-trunk structures, some associated with IR sources, observed in the region (Massi et al. 2015). Decataldo et al. (2019) numerically modelled the photoevaporation of Jeans-unstable molecular cores, hence including gravity, finding that Jeans-unstable cores will remain Jeans-unstable after the RDI, although losing part of their mass. This indicates that the most massive cores in the cloud PDRs towards NGC 6357 should collapse and host forming stars, rather than be just photoevaporated. In other words, accelerated or even triggered star formation is expected around the HII regions.

Fig. 16 η-ν (see text) plot for the cores associated with the PDRs around the HII regions G353.2+0.9 (filled red squares) and G353.1+0.6 (open blue squares). The dash dotted lines mark the critical densities for turbulent cores for different values of the non-thermal to thermal pressure ratio (α) and the magnetic to thermal pressure ratio (β), and a frozen magnetic field. The solid lines mark the evolution of the column density ratio η against the parameter ν (which increases with increasing G0) for pre-existing cores embedded in an interclump medium. The dotted lines indicate the fraction of mass lost by such cores due to photoevaporation when exposed to a steadily increasing G0. The dashed line marks the critical density for such cores, i. e. the column density of a core that would be completely photoevaporated by a FUV field ν.

4.4 The CMF and comparison with the IMF

We derived the CMF for the samples of cores in the active and quiescent region and estimated the slope of the high-mass end of the CMF, α = −d log(dN∕dM)∕d log(M), using only cores with M ≥ 5 M_⊙. To minimise statistical biases, which affect simple linear fits in the log-log space (see Maschberger & Kroupa 2009), we used the python task PyMC and the maximum likelihood estimator (ML) discussed by Maschberger & Kroupa (2009). The results are listed in Table 3.

Three trends can be pointed out from Table 3: (i) in each region (quiescent and active), all samples and estimates are consistent with one another within 1σ, except for the values for pre-stellar cores in the active area yielded by the ML; (ii) all estimates in the quiescent area are consistent with a Salpeter IMF (α = 2.35) within 1σ; (iii) all estimated CMFs in the active area are flatter than in the quiescent area at least at a 2σ significance level. Any effect of the FUV field on the dust opacity, discussed in Sect. 3.5, should not significantly affect α in the pre-stellar samples (except possible global shifts in mass, different for each region), which do not include cores embedded in the PDRs surrounding G353.2+0.9 and G353.1+0.6. We also note that the CMFs in the quiescent area displays slopes α that are consistent with those of the IMFs derived by Massi et al. (2015) for Pismis 24, adopting the smaller extinction range. When the cluster members are selected in the extinction range A_V ~ 3.2−7.8 mag, by assuming the same distance as in this work, Massi et al. (2015) find slopes in the range α ~ 1.9−2.4 for stellar masses between ~3−6 and ~60−80 M_⊙, depending on the estimator used.

The flatter CMF in the active region would be another signature of the effects of an intense FUV field, already discussed for the cores in the PDRs in Sect. 4.3. The FUV field would photodissociate the smaller cores more quickly, resulting in a top-heavy CMF. However, the effect would be more noticeable in the low-mass end of the CMF, which cannot be probed by our observations. This does not necessarily imply that a flatter CMF would in turn result in a flatter IMF as the FUV field may affect the SFE of the cores as well.

We have assumed isothermal cores. The development of temperature gradients may impact on core mass determination in two ways. Protostellar cores are heated from the inside, thus using a single dust temperature, a mean value weighted by the coldest outermost layers, can lead to overestimating a core mass. More massive cores should be more affected than less massive cores, so a residual protostellar contamination of our starless sample in the quiescent (lower FUV-flux) region may cause the actual CMF to be slightly steeper than found. On the other hand, cores in the active (high FUV-flux) region are heated from the outside, so a single dust temperature is a mean value dominated by the warmer outermost layers, resulting in an underestimate of the mass. A prototellar contamination of the sample would result in cores with both outer and inner heating and would probably just mitigate the effect. Thus the actual CMF in the UV-active region might be even flatter than found. Temperature gradients in cores are therefore likely to enhance the difference between CMFs rather than smooth it.

4.5 Cloud-cloud collision scenario

Fukui et al. (2018a; erratum in Fukui et al. 2018b) propose that the miniburst of star formation in NGC 6334 and NGC 6357 originated in a cloud-cloud collision scenario. It is therefore interesting to check whether the CMFs that we have derived show any signature of such a process, provided this is still at work. This scenario has been theoretically investigated through numerical simulations by a number of authors (e.g. Takahira et al. 2014, 2018; Balfour et al. 2015; Shima et al. 2018; Wu et al. 2018; Fukui et al. 2021) generally finding that cloud-cloud collisions increase the SFE and boost the birth of high-mass stars. Takahira et al. (2018) find a CMF with α = 1.6 in the high-mass end for a low collision speed (5 km s⁻¹), which is consistent with our determinations in the active region. However, when the collision speed is increased, a break will develop on the CMF (300 M_⊙ at v = 10 km s⁻¹ shifting to lower masses for higher collision velocities) above which the CMF steepens. Shima et al. (2018) derived an IMF steeper than a Salpeter one in the high-mass end (at collision velocities of 10 and 20 km s⁻¹). The inclusion of photoionising feedback by these authors increases the SFE and results in a slightly flatter IMF, although still steeper than a Salpeter one.

Fukui et al. (2018a) find a velocity separation between clouds of ~ 12 km s⁻¹) towards NGC 6357, so a comparison between simulations and the results summarised in Table 3 shows that the CMF in the active region is consistent with the cloud-cloud collision scenario simulated by Takahira et al. (2018), but not with the simulations of Shima et al. (2018). Nevertheless, the latter agrees with the high-end IMF derived by Massi et al. (2015) in Pismis 24 using the wider extinction range, which is actually steeper than a Salpeter one. In addition, Takahira et al. (2018) predict much higher fractions of gas mass in dense cores than we estimated in NGC 6357 (1.4%).

On the other hand, Fukui et al. (2021) found a top-heavy CMF for a collision speed v = 20 km s⁻¹, which appears to be flat up to 6–10 M_⊙ at ≥ 0.6 Myr from the beginning of the collision. However, the CMF steepens for M > 6−10 M_⊙. Our derived CMF does not exhibit such a break, but because of the completeness limit at ~ 5 M_⊙ we cannot be conclusive. Fukui et al. (2021) also propose that a signature of a cloud-cloud collision can be found in the distribution of core-core separation. We constructed this distribution by using the edges of the minimum spanning tree, as in Fukui et al. (2021). Unfortunately, our spatial resolution does not allow us to resolve their peak at ~ 0.1 pc, if any. Our distribution is nearly flat from 0 to ~50 pc (the complex radius), with a low peak between 10 and 20 pc that may resemble the second peak yielded by the simulations of Fukui et al. (2021). Ultimately, we cannot draw any clear conclusion on a cloud-cloud collision scenario in NGC 6357 based only on the comparison between observed and theoretical CMFs.

5 Summary and conclusions

We have used JCMT SCUBA-2 observations at 450 and 850 μm, complemented with Herschel (Hi-Gal) maps at 70 and 160 μm, to study theproperties of the dense core population in the galactic HII regions-young star clusters-star forming complex NGC 6357. In particular, we aimed to assess the effect of intense FUV radiation on dense cores.

We mapped the region in the ¹²CO(3–2) line as well, to correct the emission at 850 μm from line contamination in order to construct a map of pure (dust) continuum emission. We then used the algorithm Gaussclumps on the map at 850 μm to identify the dense cores after removing the extended emission.

We retrieved 1221 cores. For 686 of these we could derive a dust temperature by fitting greybodies to their SEDs obtained from submm (SCUBA-2) and FIR (Herschel) fluxes. The observed (beam-convolved) sizes lie in the range ~ 13′′ (the beam size at 850 μm) to ~ 30′′ (~ 0.1 to 0.2 pc); 452 of these have a deconvolved size of at least half the beam size. We divided the mapped area into two sub-fields, designated as active and quiescent regions, which are roughly the eastern and the western halves of the NGC 6357 complex, respectively. The active region contains the three HII regions G353.2+0.9, G353.2+0.6, and G353.2+0.7 and their associated star clusters. The cores in this region (411) are expected to be affected by an intense FUV field. The quiescent region contains 275 cores farther from the ionising sources, hence less affected by FUV photons. We therefore used the quiescent region as a control sample to study the effects of the FUV radiation on dense cores in the active region.

Our main results are the following:

• 1.

  We found different core temperature distributions in the two regions, peaking at ~25 K (quiescent region) and ~35 K (active region), which we attribute to the effects of the FUV radiation.

• 2.

  We further subdivided our sample into starless and protostellar cores, by exploiting the association with protostellar tracers such as WISE point sources and emission at 70 μm. We found 51 candidate protostellar and 224 starless cores in the quiescent region, and 160 candidate protostellar and 251 starless cores in the active region. Unfortunately, the sample of the active region is biased towards protostellar cores due to the effects of the external FUV field which heats the dust increasing its emission at 70 μm and saturation in the WISE bands towards the HII regions.

• 3.

  We derived the core masses by assuming a plausible dust opacity (we also discussed possible effects of FUV radiation). We estimated a mass completeness limit of ~5 M_⊙. We found mass-size relations of log(M∕M_⊙) ~ a × log(D∕arcsec), with a in the range 2.0−2.4, consistent with those from other regions or large samples of clumps. The starless cores above the mass completeness limit are likely to be gravitationally bound, hence pre-stellar in nature.

• 4.

  The estimated fraction of molecular gas in dense cores is 1.4% in NGC 6357. This is consistent with values from most galactic star formation regions.

• 5.

  We showed that the properties of the cores nearer to the HII regions are consistent with pre-existing cores gradually being engulfed in a PDR and photoevaporating.

• 6.

  We constructed the CMFs for M > 5 M_⊙, finding a Salpeter-like CMF in the quiescent region (α ~ 2.1−2.4) and a significantly flatter (at a 2 σ level) one in the active region (~1.7−1.9 with α = −d log(dN∕dM)∕d log(M)). The difference becomes even more significant for pre-stellar cores when using the maximum likelihood estimator (α ~ 2 vs. α ~ 1.2). We attribute this to the effects of the FUV radiation as well.

• 7.

  We compared the CMFs with those predicted by simulations of cloud-cloud collisions, finding no conclusive evidence for cloud-cloud collisions giving rise to the cores, rather than them being pre-existing.

In conclusion, we found differences in the global properties of the cores in the active region and those in the quiescent region. In particular we found a statistically significant difference in the slope of the CMF between the two regions. We attribute these differences to the influence of the FUV radiation.

Acknowledgements

The James Clerk Maxwell Telescope was operated by the Joint Astronomy Centre until 2015, on behalf of the Science and Technology Facilities Council of the United Kingdom, the National Research Council of Canada and the Netherlands Organisation for Scientific Research. Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for Innovation. The James Clerk Maxwell Telescope has been operated by the East Asian Observatory since 2015, on behalf of The National Astronomical Observatory of Japan, Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, the National Astronomical Observatories of China and the Chinese Academy of Sciences (Grant No. XDB09000000), with additional funding support from the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada. We thank Harriet Parsons for the Perl script used to correct the 850 μm data with the CO-emission, and Eugenio Schisano for providing us with the Herschel 160 μm map used inFig. 1. This research has made use of NASA’s Astrophysics Data System Bibliographic Services (ADS) and of the SIMBAD database, operated at CDS, Strasbourg, France. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. A.G. and F.M. were partly supported by INAF through the grant Fondi Mainstream ’Heritage of the current revolution in star formation: the Star-forming filamentary Structures in our Galaxy’.

Appendix A Additional figures

Here we show examples of SED fits to cores in the active and quiescent regions (see Sect. 4.2).

Fig. A.1 Example of SED fits to two cores in the active region. Top panel: all flux densities can be fitted with a single greybody; there is no excess at 70 μm, and this core is recognised as a pre-stellar core. The core in the lower panel has an excess flux density at 70 μm (the point has not been excluded from the fit, hence the large χ2). This could be due either to the core being heated from the outside by the FUV radiation field, or to a protostar inside the core (see Sect. 4.2).

Fig. A.2 Same as Fig. A.1, but for two cores in the quiescent region.

References

All Tables

All Figures

	Fig. 1 Herschel HiGAL 160 μm map (colour scale) with fields observed in 2012 (1–6; green) and 2015 (7–12; blue). The location and approximate size of star clusters are marked in yellow (Massi et al. 2015). The linear size of the region shown is approximately 22 × 17 pc2.
In the text

	Fig. 2 12CO(3–2) integrated emission ∫ TMB dV; − 20 to + 20 km s−1). Contours values are 80(80)470 Kkm s−1 (low(step)high).
In the text

	Fig. 3 SCUBA-2 observations of the dust emission. Left: 850 μm mosaic; the emission has been corrected for the integrated 12CO(3–2) line emission. Right: 450 μm mosaic. Units are mJy beam−1.
In the text

	Fig. 4 850 μm map after filtering out emission from extended structures as explained in the text. Units are mJy beam−1. The black lines indicate the region mapped in 12CO(3-2) (see Fig. 2).
In the text

	Fig. 5 Distribution of fluxes for the 1221 cores identified in the 850 μm mosaic.
In the text

	Fig. 6 Distribution of core sizes for all 1221 cores (uncorrected for beam; black solid line); the red curve indicates the active region and the blue one the more quiescent part.
In the text

	Fig. 7 All 686 cores with a valid SED fit superimposed on the mosaic at 850 μm. Darker symbols indicate lower temperatures. The white line indicates the dividing line between active (to the east) and quiescent (to the west) regions.
In the text

	Fig. 8 Temperature distribution for all 686 cores (solid black line), cores in the active region only (411; solid red line), and cores inthe quiescent region only (275; solid blue line).
In the text

Fig. 9 Core mass function in NGC 6357 (dashed green curve for all cores) obtained as density distribution of sources after smoothing with a Epanechnikov kernel. The red curve indicates the active region, and the blue curve indicates the more quiescent part. The vertical line marks the 5 M⊙ completeness limit (see Sect. 3.6). The small coloured bars show the masses of the single sources for each sample, following the colour coding of the distributions.

In the text

Fig. 10 Core masses and core mass distribution. (a) Core masses derived using a density-dependent opacity (grains with no ice mantles) vs. masses derived using a constant dust opacity of 1.93 cm 2 g−1 (at 850 μm). The full line marks equal masses, whereas the dashed line marks density-dependent opacity masses equal to 70% of the constant density masses. (b) Core mass distributions for masses obtained with a constant dust opacity of 1.93 cm 2 g−1 (at 850 μm; full line) and for masses obtained with a density-dependent opacity (dashed line). The error bars have been computed assuming a Poisson statistics.

In the text

	Fig. 11 Completeness of the cores in our sample, as a function of mass. Starting from the completeness as a function of flux density for point sources, we have used a representative source size of 22′′ and the observed temperature distribution (Fig. 8) to derive the distribution in mass shown here (see text). The dashed line indicates the mass where the completeness reaches its maximum.
In the text

	Fig. 12 Continuum emission at 850 μm towards G353.2+0.9, overlaid with the locations of the gas cores found by Giannetti et al. (2012; red circles of size equal to that of a clump) and of the dense cores retrieved in this work (black ellipses). The gas clumps are labelled according to Giannetti et al. (2012). Units are mJy beam−1. Contour levels are 50, 150, and 300 mJy beam−1.
In the text

	Fig. 13 Continuum emission at 850 μm towards the fields of cores MDC 1, 2, and 5 (labelled in figure) of Russeil et al. (2019). Units are mJy beam−1. The locations of the cores identified in this work are overlaid (small open squares).
In the text

	Fig. 14 Distribution of temperatures for cores in the active (full line) and quiescent (dashed line) regions. Only cores with M > 5 M⊙ (i.e. above the mass completeness limit) have been included. (a) All cores. (b) Cores with no emission at 70 μm and no red WISE sources associated. (c) Cores associated with 70 μm emission and/or red WISE sources.
In the text

Fig. 15 Massvs. deconvolved size for (a) all cores; (b) 70 μm dark cores cleaned from WISE-associated cores; (c) 70 μm bright cores plus cores associated with WISE sources. Cores with a deconvolved size smaller than half the beam size have been discarded. Small filled squares mark cores in the active region, whereas open green triangles mark cores in the quiescent region. The critical half-mass loci of Bonnor-Ebert spheres are drawn for T = 20 K and T = 30 K (black dash-dotted lines). The dotted red line marks the completeness limit for T = 20 K dust (from Eq. (3)). Linear (log-log) fits are indicated by a solid (cores in the quiescent region) and a dashed (cores in the active region) line.

In the text

Fig. 16 η-ν (see text) plot for the cores associated with the PDRs around the HII regions G353.2+0.9 (filled red squares) and G353.1+0.6 (open blue squares). The dash dotted lines mark the critical densities for turbulent cores for different values of the non-thermal to thermal pressure ratio (α) and the magnetic to thermal pressure ratio (β), and a frozen magnetic field. The solid lines mark the evolution of the column density ratio η against the parameter ν (which increases with increasing G0) for pre-existing cores embedded in an interclump medium. The dotted lines indicate the fraction of mass lost by such cores due to photoevaporation when exposed to a steadily increasing G0. The dashed line marks the critical density for such cores, i. e. the column density of a core that would be completely photoevaporated by a FUV field ν.

In the text

Fig. A.1 Example of SED fits to two cores in the active region. Top panel: all flux densities can be fitted with a single greybody; there is no excess at 70 μm, and this core is recognised as a pre-stellar core. The core in the lower panel has an excess flux density at 70 μm (the point has not been excluded from the fit, hence the large χ2). This could be due either to the core being heated from the outside by the FUV radiation field, or to a protostar inside the core (see Sect. 4.2).

In the text

	Fig. A.2 Same as Fig. A.1, but for two cores in the quiescent region.
In the text