The dynamic centres of infrared-dark clouds and the formation of cores

ABSTRACT

High-mass stars have an enormous influence on the evolution of the interstellar medium in galaxies, so it is important that we understand how they form. We examine the central clumps within a sample of seven infrared-dark clouds (IRDCs) with a range of masses and morphologies. We use 1-pc-scale observations from the Northern Extended Millimeter Array (NOEMA) and the IRAM 30m telescope to trace dense cores with 2.8-mm continuum, and gas kinematics in C¹⁸O, HCO⁺, HNC, and N₂H⁺ (J = 1–0). We supplement our continuum sample with six IRDCs observed at 2.9 mm with the Atacama Large Millimeter/submillimeter Array (ALMA), and examine the relationships between core- and clump-scale properties. We have developed a fully automated multiple-velocity component hyperfine line-fitting code called mwydyn which we employ to trace the dense gas kinematics in N₂H⁺ (1–0), revealing highly complex and dynamic clump interiors. We find that parsec-scale clump mass is the most important factor driving the evolution; more massive clumps are able to concentrate more mass into their most massive cores – with a log-normally distributed efficiency of around 9 per cent – in addition to containing the most dynamic gas. Distributions of linewidths within the most massive cores are similar to the ambient gas, suggesting that they are not dynamically decoupled, but are similarly chaotic. A number of studies have previously suggested that clumps are globally collapsing; in such a scenario, the observed kinematics of clump centres would be the direct result of gravity-driven mass inflows that become ever more complex as the clumps evolve, which in turn leads to the chaotic mass growth of their core populations.

1. INTRODUCTION

The stellar populations of the many billions of galaxies in the Universe are dominated, in absolute number, by low-mass stars (⁠|$m_* \lesssim 2\ \mathrm{M}_\odot$|⁠), while the much rarer high-mass stars (⁠|$m_*\ \gt\ 8\ \mathrm{M}_\odot$|⁠) have an enormous influence on the interstellar medium (ISM) and chemical evolution of galaxies. The relative abundance of low- through to high-mass stars as they enter the main sequence is given by the stellar initial mass function (IMF), and our understanding of the evolution of galaxies depends critically upon our understanding of how the IMF arises.

For low-mass stars, observations of nearby star-forming regions have revealed that the gravitational fragmentation of filaments with a (super)critical mass-per-unit-length appears to play a crucial role in determining the masses of 0.1-pc-scale cores (André et al. 2014). These cores are thought to be the precursors of individual stellar systems, and the distribution of core masses – the core mass function (CMF) – is found in these regions to closely resemble the IMF, leading to speculation that the latter is inherited directly from the former (e.g. André et al. 2010, 2014; Polychroni et al. 2013; Könyves et al. 2020; Ladjelate et al. 2020). However, for high-mass stars, the picture is more complicated (see Motte, Bontemps & Louvet 2018b, for a recent review). High-mass star-forming regions are inherently more difficult to observe as a consequence of both their short lifetimes (∼0.5–2 Myr, Battersby, Bally & Svoboda 2017; Sabatini et al. 2021), and their relatively low rate of occurrence. The relatively small number of observations of high-mass star-forming regions that are able to robustly measure the CMF finds them to deviate significantly from the IMF, with top-heavy distributions (e.g. Motte et al. 2018a; Kong 2019), and fragmentation studies (e.g. Bontemps et al. 2010; Sanhueza et al. 2017) routinely demonstrate that Jeans-like fragmentation is unable to produce cores that could be the progenitors of high-mass stars, indicating a different formation pathway.

One of the most massive protostellar cores that has been observed in the Galaxy, SDC335-MM1 (⁠|$\sim\!{500}\,$| M_⊙, Peretto et al. 2013), is located at the centre of a so-called hub-filament system (HFS, Myers 2009): a focal point of converging filaments of dense gas. A growing sample of HFSs has since been studied in the literature, which often find massive cores at the centre of the systems, velocity gradients along the filaments which may be interpreted as longitudinal gas flows that bring material to the centre of the system at rates of ∼10⁻⁴–10⁻³ M_⊙ yr⁻¹ (e.g. Kirk et al. 2013; Peretto et al. 2014; Lu et al. 2018; Williams et al. 2018; Chen et al. 2019; Anderson et al. 2021). Whether the gas flows arise as a result of global hierarchical gravitational collapse (e.g. Vázquez-Semadeni et al. 2019) or supernova feedback-powered inertia (e.g. Padoan et al. 2020) is an ongoing matter of debate, but such large flow rates are clearly sufficient to transport a significant quantity of material into the centre of the region over a few Myr, and may be a crucial aspect of high-mass star formation.

HFSs have most commonly been identified within infrared-dark clouds (IRDCs) in mid-IR survey data (Peretto & Fuller 2009) which have superior resolution offered compared with submillimetre or millimetre survey data; e.g. the 8-|$\mu$|m Galactic Legacy Infrared Midplane Survey Extraordinaire (GLIMPSE) survey data (Benjamin et al. 2003) have a resolution of 2 arcsec compared with 18 arcsec of the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL; Schuller et al. 2009). IRDCs also make excellent targets for studies of the earliest phases of star formation due to their low level of emission at mid- and far-IR wavelengths. The relative level of 8-|$\mu$|m brightness has been shown to be a particularly effective tracer of the evolutionary state of molecular structures (e.g. Battersby et al. 2011; Rigby et al. 2021), being dark at the early stages, and bright when star formation is advanced. This was quantified in a new metric, the infrared-bright fraction (f_IRB, Rigby et al. 2021; see also Watkins et al., in preparation), which measures the fraction of pixels within a structure that are brighter than the local background. Rigby et al. (2021) used this quantity to demonstrate that parsec-scale clumps within the Galactic Star Formation with NIKA2 (Adam et al. 2018) Galactic Plane Survey (GASTON-GPS) tend to grow in mass at early times, supporting the gravitational collapse scenario for high-mass star formation as in Peretto et al. (2013).

The f_IRB metric was used by Peretto et al. (2022) in conjunction with another new metric – the filament convergence parameter, f_C, that quantifies the HFS-like nature of a clump by giving systems at the centre of radially converging filaments a higher value – to examine the evolution of HFSs in the same sample of GASTON clumps. The value of f_C tends to increase with f_IRB, suggesting that either HFSs are a late stage of clump evolution, or that IR-dark HFSs have short lifetimes compared with IR-bright HFSs, or both. A similar conclusion was reached by Kumar et al. (2020), who found that high-mass star-forming regions are systematically associated with HFSs. The latter two studies are, crucially, based on single-dish data, and are limited to angular resolutions of ≳10 arcsec, which prevent core scales to be resolved for typical distances of a few kpc for high-mass star-forming regions.

At higher angular resolution, Anderson et al. (2021) compared the properties of 0.1-pc-scale cores within a sample of six IR-dark HFSs with the cores within a larger sample of clumps mapped with the Atacame Large Millimeter/submillimeter Array (ALMA), and found that IR-dark clumps tended to have the highest efficiencies for formation of the most massive cores (MMCs), i.e. a greater fraction of the total clump mass was contained within the single MMC, than for IR-bright ones. Similarly, Traficante et al. (2023) found that the surface density of clumps and the MMCs increase with evolutionary stage. Studies from the ALMA-IMF large programme (Motte et al. 2022), Nony et al. (2023) are finding a time dependence in the evolution of CMFs in W43 (Pouteau et al. 2023), and Nony et al. (2023) found that the protostellar CMF in W43 is more top-heavy than the Salpeter-like pre-stellar CMF, indicating that the MMCs accrete a larger fraction of mass than less massive ones. Taken together, these results suggest that the MMCs grow most rapidly at earlier times, while the rest of the core population catches up at later times as the clumps continue to gain material from their wider environment. However, we lack a detailed picture of how the dense gas in the centres of hub-filaments facilitate this growth of the MMCs, and if and how these conditions differ from those found within non-HFS clumps, and across the range of masses. In this paper, we attempt to address this by examining the relationship between core masses and dense gas kinematics traced by N₂H⁺ across a sample of IR-dark clumps with a range of properties.

The fundamental rotational transition of the molecular ion diazenylium, N₂H⁺, at 93 GHz was first detected by Turner (1974), and later identified by Green, Montgomery & Thaddeus (1974), and confirmed by Thaddeus & Turner (1975) and has since been widely utilized as a tracer of high column density molecular gas. It is usually optically thin, and thus is ideal for studying the kinematics of dense gas linking them to the formation of dense cores in this study. Its rotational transitions – including the J = 1–0 transition at 93.174 GHz that we have observed in this study – have at least seven hyperfine components (Caselli, Myers & Thaddeus 1995), though they often appear as a triplet in the ISM where linewidths are typically supersonic. The high-frequency F₁, F = 0, 1 → 1, 2 component (sometimes called the ‘isolated component‘) is located at a frequency equivalent to a Doppler-shifted velocity of −8 km s⁻¹, and a low-frequency triplet whose component centroids extend to ∼+7 km s⁻¹ with respect to a central triplet that defines the rest frequency.

Due to this hyperfine structure, maps of the first and second moment of the full emission line cannot be interpreted as simply as is the case for molecular emission lines that have no hyperfine splitting, such as CO. Performing fits of models to the hyperfine components in order to recover centroid velocities, linewidths and excitation temperatures has routinely been done within the literature with observations of N₂H⁺; however, the high-spatial resolutions that are now achievable with facilities like the Northern Extended Millimeter Array (NOEMA) and ALMA have meant that multiple separate velocity components can be distinguishable within the spectra, which complicates fitting considerably. Some studies have achieved this by ignoring the majority of the hyperfine structure (e.g. Henshaw et al. 2014; Hacar et al. 2018; Barnes et al. 2021), focusing their efforts on the isolated component, which may then be fit using a combination of a small number of Gaussian profiles. This is only achievable when the isolated components are clearly identifiable, which is often not the case in regions of very high column density and linewidth such as those within this study. Furthermore, the isolated component is a factor of a few weaker than the brightest component, except in optically thick regions, and so provides a reduced effective signal-to-noise ratio (SNR) compared with fitting the full spectrum. To overcome these problems, we have developed a fully automated multiple velocity component hyperfine line-fitting code to analyse N₂H⁺ cubes, called mwydyn (the Welsh word for worm, pronounced ‘muy-din’, in recognition of the wiggly nature of the spectra).

In this paper, we examine the centres of a sample of IRDCs with existing NIKA (Monfardini et al. 2010) and NIKA2 (GASTON-GPS) observations covering a range of properties, using high-resolution observations of N₂H⁺ (1–0) alongside several other key molecules and continuum, in order to examine the kinematics surrounding the MMCs, and understand how the gas motions relate to the assembly of material within them. In Section 2, we described the sample selection and observations. In Section 3, we describe the analysis of the data, including a description of a new fully automated fitting code for lines with hyperfine structure, and extraction of core- and clump-scale properties. Our main results are described in Section 4, and discussed in Section 5, and we present our conclusions in Section 6.

2. TARGETS AND OBSERVATIONS

2.1 Sample selection

We compiled a sample of seven IRDCs from the catalogue of Peretto & Fuller (2009) to examine with NOEMA, which are a subsample of the study of Peretto et al. (2023). The IRDCs were selected to cover a range of masses and morphologies, and to lie at comparable distances with 3.5 < d/kpc < 5.0 thus limiting the impact of a varying spatial resolution element. We also have overlapping 1.2-mm observations from our NIKA (Rigby et al. 2018) and NIKA2 (Rigby et al. 2021) observing programmes for these seven IRDCs. Systemic velocities for each of the IRDCs in our sample had been determined by observations of N₂H⁺ (1–0) in Peretto et al. (2023) using EMIR on the IRAM 30m telescope. The systemic velocities were then used to determine heliocentric distances using the batch input version (v2.4.1) of the Reid et al. (2016)¹ Bayesian Distance Calculator, which adopts the latest values of fundamental Galactic parameters from Reid et al. (2019), including a Sun–Galactic Centre distance of 8.15 kpc. We make the following modifications to the configuration files: (i) we insist that, since our targets are IRDCs, only the near kinematic distance solution is considered, by setting P(far) to zero; (ii) since we do not know a priori that the IRDCs reside within spiral arms, we disable the influence of the prior relating to the spiral arm model upon the distance probability density functions. The result is that our distances are primarily given by the near kinematic distance, but incorporating information for any maser sources with measured parallaxes. Details of the observations are given in Table 1, including parsec-scale properties of the central clumps that are characterized in Section 3.2.

2.2 Observations

The column density peaks of the seven IRDCs – which we refer to hereafter as the central clumps – were observed with single NOEMA pointings between 2019 May 21 and August 30 using the PolyFiX wideband correlator in Band 1 (3 mm), providing ∼31 GHz instantaneous bandwidth split over two sidebands each consisting of two polarizations. The local oscillator was tuned to a frequency of 99.5 GHz, providing spectral windows covering approximately 88–96 and 103–111 GHz. We selected several high-resolution spectral windows with 62.5 kHz channel widths in order to provide 0.2 km s⁻¹ channels for our main target lines which were the (1–0) transitions of HCO⁺, HNC, N₂H⁺, and C¹⁸O.

A total of nine 15-m dishes were used in the nine-dimensional (9D) configuration, which used stations W12, W09, W05, E10, E04, N13, N09, N05, and N02 with projected baselines ranging from 15 to 176 m. The receiver bandpass was calibrated using the quasars 1749+096, 1830−210, 3C454.3, 3C345, and 3C273. Phase and amplitude calibration was performed using observations of the quasars 1829–106, 1830–210, and 1908–201. The star MWC 349 was used as the primary flux calibrator, which is accurate to typically less than 10 per cent² in the 3-mm band. Short-spacing observations for our spectral lines were obtained using the IRAM 30m telescope with dedicated observing runs in 2019 November, and 2020 February, May, July, and October. For the short-spacing observations we used EMIR with a set-up optimized for our main target lines, with the local oscillator tuned to 99.5 GHz, and the FTS back-end in 50-kHz mode to best match our NOEMA observations.

The 2.8 and 3.2 mm continuum images, and integrated intensity images from the main targeted lines are shown in Fig. 1, along with Spitzer/GLIMPSE 8 |$\mu$|m (Benjamin et al. 2003), Herschel/Hi-GAL 70 |$\mu$|m (Molinari et al. 2016), and IRAM 30-m/NIKA and GASTON 1.2 mm (Rigby et al. 2018, 2021) continuum images showing the wider environment.

2.3 Data reduction

Data reduction was performed using software from the gildas suite.³ The raw data were calibrated using the clic software, resulting in calibrated uv-tables, which were then imaged and cleaned using mapping. For the continuum sidebands at 2.8 and 3.2 mm, the procedure involved, first, producing a cleaned cube in order to identify and filter any spectral lines within the band that may result in overestimated continuum flux densities, before imaging the cube with filtered channels, and cleaning.

Pseudo-visibilities for the single-dish data were combined with the NOEMA uv-tables for the four main target lines before imaging. To clean the images, we used robust weighting with a ‘robust’ parameter of 0.5 to achieve a compromise between point-source sensitivity and resolution. We did not use any cleaning support in order to prevent the introduction of unwanted bias into the resulting emission maps. Two continuum images per source were also produced following the same procedure, with the exception that there are no short-spacing observations. The half power primary beamwidth, which defines the field of view, is 54.1 arcsec for the N₂H⁺ (1–0) and 3.2-mm continuum observations, for which the synthesized beam is typically 6.6 × 3.2 arcsec across our targets. For the 2.8-mm continuum observations, the average beam size is 5.6 × 2.7 arcsec, with a field of view of 47.0 arcsec. The individual beam sizes and rms values for each target are listed in Table 1.

2.4 Supplementary high-resolution continuum data

We also supplement our NOEMA continuum observations with 2.9-mm continuum observations of six IRDCs mapped previously with ALMA, as described in Anderson et al. (2021): SDC326.476, SDC335.579, SDC338.315, SDC339.608, SDC340.969, and SDC345.258. These supplementary observations have a greater spatial extent, and have a higher angular resolution (with beams of typically 3 × 2 arcsec) than the observations presented in this study, having fully mapped the clouds, and so we first processed the ALMA observations to be more directly comparable to our single-pointing NOEMA data: The ALMA data were spatially smoothed to the average beam size for our NOEMA observations. The observations were then resampled onto a pixel grid of the same pixel size and extent (i.e. cropping to the NOEMA primary beamwidth) of our NOEMA maps, replicating the same field of view as our 2.8-mm observations, with the equivalent of the single-pointing centres being placed at the peak intensity of the Herschel Infrared Galactic Plane Survey (Hi-GAL) 500-|$\mu$|m observation (Molinari et al. 2016) of the same cloud. This latter choice was made to determine the pointing centres in a similar way than was done for our NOEMA sources, which were based on the peak intensities from our 1.15-mm NIKA and NIKA2 data. The characteristics of the resulting processed data, including beam sizes, can be found in Table 1.

2.5 Ancillary data

We make use of 8-|$\mu$|m continuum imaging data from the Spitzer GLIMPSE (Churchwell et al. 2009), which have an angular resolution of 2 arcsec. We also make use of 70-, 160-, and 250-|$\mu$|m imaging from the Hi-GAL (Molinari et al. 2016), and column density maps generated from these images from Peretto et al. (2016). The 70-, 160-, and 250-|$\mu$|m data have angular resolutions of ∼8, 12, and 18 arcsec, respectively, and the column density maps have the same resolution of the 250-|$\mu$|m images. Finally, we made use of 1.15-mm NIKA2 (Adam et al. 2018) imaging from the GASTON-GPS (Rigby et al. 2021), as well as the NIKA maps of SDC18.888 and SDC24.489 from Rigby et al. (2018). Both the NIKA (Monfardini et al. 2010) and NIKA2 data have an angular resolution of 11 arcsec.

3. ANALYSIS

3.1 Core population

The population of compact sources, within each field in the 2.8-mm continuum data was extracted using astrodendro,⁴ a Python-based implementation of the dendrogram (Rosolowsky et al. 2008) segmentation method. Sources were extracted using a minimum valid flux density level of 3 times the rms noise value, with a minimum significance for structures of 2 times the rms noise value, meaning that the minimum peak intensity of 5 times the rms noise value. For each map, the rms noise value was first determined using emission-free regions of the cleaned image. Finally, sources were required to have an area equal to at least half of the beam area. The latter choice was made as opposed to requiring a full beam area because, by construction, the dendrogram technique always underestimates compact source sizes as a result of clipping the wings of the source profile after convolution with the telescope beam, an effect that is particularly pronounced for point sources near the detection limit, or in crowded regions. To alleviate this problem, aperture corrections were applied to the core fluxes that are a function of the source area and peak SNR, which are detailed in Appendix  A. We also find that the inclusion of this minimum source size assists with breaking down oversized sources during the construction of the dendrogram, and recovering compact sources more accurately in crowded regions. The 2.8-mm data were favoured over the 3.2-mm data for source detection due to both their slightly higher angular resolution, and the greater intensity of dust continuum emission, at the cost of a slightly reduced field of view.

The source extraction procedure produces a mask for each of the sources it identifies, and we used these masks to calculate the properties of interest: integrated flux density, peak flux density, SNR, and radius. We are particularly interested in the dendrogram ‘leaves’, which are smallest structures within which no further substructure is discernible. Hereafter, we will refer to these objects as ‘cores’, although we note that the beam size gives a spatial resolution of 0.1 pc at the mean IRDC distance, and so these cores will almost certainly contain further unresolved substructures within them.

Core masses were calculated using:

where S_ν is the flux density integrated over the source area, d is the distance to the source, B_ν(T_d) is the Planck function evaluated at frequency ν and dust temperature T_d, and the dust opacity is given by: |$\kappa _\nu = 0.1 \, \mathrm{cm}^2 \, \mathrm{g}^{-1} (\nu /\mathrm{THz})^{\beta }$|⁠. This value incorporates a gas-to-dust mass ratio of 100, and is the typical value adopted by e.g. Marsh, Whitworth & Lomax (2015), Peretto et al. (2016), and Anderson et al. (2021). For consistency with the latter study in particular, we adopt a value for the dust emissivity spectral index, β = 2.

To determine the dust temperatures, we follow the same approach as Anderson et al. (2021). Firstly, for all cores we adopt the flux density-weighted colour temperature from the maps of Peretto et al. (2016) in which the temperature is estimated from the Herschel/Hi-GAL 160-/250-|$\mu$|m colour, again assuming β = 2. This technique is limited to the resolution of the 250-|$\mu$|m data, which is 18 arcsec, and so we may overestimate the temperature of cold cores compared with the ambient dust temperature on larger scales. For cores containing any embedded star formation, we estimate the core temperature from the 70-|$\mu$|m luminosity of any associated compact sources (i.e. protostars) within the catalogue of Molinari et al. (2016). The bolometric luminosity of compact sources is known to be strongly correlated with 70-|$\mu$|m luminosity (Dunham et al. 2008; Ragan et al. 2012), and so we use the relationship of Elia et al. (2017):

to calculate bolometric luminosities, where S₇₀ is the integrated flux density of the associated 70-|$\mu$|m compact source. For optically thin dust heated by an internal star, the temperature at a particular radius can then be determined as (Terebey, Chandler & Andre 1993)

where q = 2/(4 + β), β = 2, and T₀ = 25 K is the dust temperature at a reference radius r₀ = 0.032 pc for a star of luminosity L₀ = 520 L_⊙. For the radius of our cores, we adopt the ‘effective’ radius of a circle with the source area (i.e. the sum of the area of all pixels in the leaf), |$r = R_\mathrm{eff} = \sqrt{A/\pi }$|⁠, deconvolved by the NOEMA beam.

In this way, we determined the masses of a total of 47 cores, with 21 and 26 cores from the NOEMA and ALMA data sets, respectively, and we present their properties in Table 2. The cores range in mass between ∼2 and 1300 M_⊙, where the latter value belongs to SDC335.579 MM1, a well-known high-mass core (Peretto et al. 2013). With 10 per cent and 5 per cent uncertainties on the integrated flux densities for the NOEMA and ALMA data, respectively, a 0.5 K uncertainty on colour temperatures (as recommended by Peretto et al. 2016), and errors on the distance determinations that are typically ∼1 kpc, the uncertainties are dominated by a factor of 2 uncertainty on κ₀ (Ossenkopf & Henning 1994). The core masses are therefore thought to be accurate to within a factor of ∼2.

3.2 Clump-scale properties

We calculated the masses of the central clumps within the 13 IRDCs from the Peretto et al. (2016) column density maps⁵ in order to provide context of the wider environment for our observations. The mass was measured within a 1-pc-diameter aperture centred on the positions given in Table 1 – which is approximately the diameter of the NOEMA primary beam full width at half maximum (FWHM) for our pointings – since it is difficult to describe another definition that works well in the varying backgrounds across the sample. This is the approximate size-scale of what are typically called ‘clumps’ (e.g. Ellsworth-Bowers et al. 2015; Urquhart et al. 2018; Elia et al. 2021), though we will hereafter refer to this particular measurement as the ‘1-pc clump mass’ (M_1pc), as a reminder of our fixed-aperture calculation. Because this technique adopts a fixed aperture size, the measurement is equivalent to the average surface density. 1-pc clump masses are calculated according to

where the surface area element dA = d²dΩ, for the source distance d and pixel solid angle dΩ. In this case, the column densities used have been background-subtracted by first subtracting the minimum value of the column density within the aperture from each pixel: |$N_{\mathrm{H}_2}^\prime = N_{\mathrm{H}_2}-N_{\mathrm{H}_2\mathrm{, bg}}$|⁠. We adopt a value of 2.8 for the mean molecular weight per hydrogen molecule, |$\mu _\mathrm{H_2}$|⁠, which is the result of assuming mass fractions of 0.71 for hydrogen, 0.27 for helium, and 0.02 for metals. The recovered 1-pc clump masses range from 330 to 4360 M_⊙, corresponding to mean surface densities of Σ_pc ∼0.1–1.2 g cm⁻² across the aperture. For context, this range of densities is encompassed by the density range of 0.05–1.0 g cm⁻² of high-mass star-forming clumps found in ATLASGAL (Urquhart et al. 2014), indicating that all of our sources are capable of forming high-mass stars.

One measure of the dynamic status of a clump is through the virial parameter, which measures the ratio of gravitational to the kinetic energy. A clump is in virial equilibrium when the gravitational energy is equal to twice the total kinetic energy, 2E_k = E_G. Following the formulation of Bertoldi & McKee (1992), we calculate the virial parameter:

where σ is the total linewidth (including thermal and non-thermal contributions), R is the radius, and M is the mass. The choice of a fixed radius in our methodology would result, in cases where the parsec-mass is contained within an area that is smaller than the fixed aperture, in an overestimated virial parameter. Therefore, we calculated an adjusted radius and mass for the virial parameter determination only. To do this, we first determined the lowest closed contour in the N₂H⁺ (1–0) integrated intensity. We adopted the radius of a circle with the same area A enclosed by the contour, |$R_\mathrm{ctr} = \sqrt{A/\pi }$|⁠. We then calculated the background-subtracted column density within this same contour from the Herschel-derived column density maps as in equation (5), and term this M_ctr, which we use as our mass measurement (and note that this measurement more closely resembles the usual ‘clump mass’). Finally, the linewidths σ_ctr were then obtained by fitting the N₂H⁺ (1–0) spectrum, averaged over the region covered by the same contour as for the mass determination. We performed the fit using the gildas: class model, which is described in more detail in Section 3.3, though we fit only a single component here. We list the values derived for the virial parameter determination in Table 3.

A clump in virial equilibrium has α_vir = 1, though we note that equipartition of kinetic and gravitational energy occurs at α_vir = 2, and thus any clump with α_vir ≤ 2 is considered to be gravitationally bound. We note that this formulation of the virial parameter incorporates a factor of order unity which accounts for a non-uniform and non-spherical mass distribution, though the equation is derived under the assumption of a uniform density sphere. This means that for a source with a radial density profile of ρ(r) ∝ r⁻², equipartition of kinetic and gravitational energy occurs at α_vir = 3.3. For the seven NOEMA clumps in our sample for which we have N₂H⁺ (1–0) linewidths, virial parameters are in the range of 0.3–1.5, and so they are all considered to be gravitationally bound in keeping with the Peretto et al. (2023) measurements of the same sources.

The bolometric luminosity of each clump is calculated using equation (2), and using the total integrated flux density of all 70-|$\mu$|m compact sources from Molinari et al. (2016) that lie within the 1-pc aperture. Bolometric luminosities calculated in this way vary from ∼15 to 50 000 L_⊙ across the sample, with all but two of the thirteen sources having L_bol > 1000 L_⊙ – the limit at which massive young stellar objects and H ii regions are associated with ATLASGAL clumps (Urquhart et al. 2014), and roughly corresponding to an embedded B3 (∼6 M_⊙) or earlier type star.

We also supplement these properties with two new quantifications of each clump’s evolutionary status. Firstly, for all 13 sources we calculate the infrared-bright fraction f_IRB, which determines the fraction of pixels within the clump at 8 |$\mu$|m (from Spitzer/GLIMPSE imaging, Churchwell et al. 2009) that are brighter than the local background within a 4.8-arcmin-wide box, and which has shown to be a good tracer of relative evolution (see Rigby et al. 2021, Watkins et al., in preparation). Sources evolve from f_IRB = 0 at the earliest stages of star formation to f_IRB = 1 at the latest stages, although we note that absolute time taken to evolve through this sequence is probably a function of mean density (via the free-fall time) and therefore mass.

Secondly, we calculate the filament convergence parameter f_C as presented by Peretto et al. (2022) for the seven NOEMA sources that have been observed by NIKA (Rigby et al. 2018) or NIKA2 (Rigby et al. 2021). The convergence parameter quantifies the level of local filament convergence associated with each pixel by identifying the filaments within a field of a given radius and then quantifying how close they are to being radially aligned with that pixel. In this case, the field radius is set to 39 arcsec, corresponding to 1 pc at a distance of 5.2 kpc, which is the median distance for clumps in the GASTON field. The filaments are identified using the second derivative method of e.g. Orkisz et al. (2019), which determines a map of topological curvature, and identifies pixels which have eigenvalues below a threshold of −3 times the local standard deviation as being associated with filaments. These filaments are then reduced to single pixel-wide ‘skeletons’ using scikit-learn’s skeletonize routine. Peretto et al. (2022) formulated the convergence parameter – which is calculated for a pixel with coordinates (x, y) – as

where N_fil and N_pix are the number of filaments, and the number of filament pixels within the search radius, respectively, Δθ is the angle between the radial direction from position (x, y) to pixel i, and the filament direction at pixel i. C_n is a normalization constant which ensures that a pixel at the convergence point of six parsec-long radially distributed filaments has a value of f_C = 1.

For the five sources falling within the GASTON-GPS field, we take the maximum value within the NOEMA fields of view from the convergence parameter map of Peretto et al. (2022), and for the two sources covered by NIKA, we calculate new convergence parameter maps in exactly the same way from the 1.2-mm maps of Rigby et al. (2018). This parameter quantifies the HFS-like nature of a clump by identifying nearby filaments, and providing a high value f_C → 1 for sources which are located at the converging point of several nearby filaments, or a low value of f_C = 0 for clumps that do not have any local filaments pointing towards themselves. We find that the NIKA and NIKA2 data are perfect for determining f_C because they represent a combination of sampling the long-wavelength part of a dust spectral energy distribution (SED) that is most sensitive to dust column density, and relatively high-angular resolution (∼11 arcsec), and so we are unable to determine an equivalent value for the ALMA sample of IRDCs. In contrast to the NIKA and NIKA2 data, the closest match in angular resolution from Herschel/Hi-GAL (Molinari et al. 2016) would be the 160- or 250-|$\mu$|m data at 12- and 18-arcsec resolution, respectively, but the continuum emission at these wavelengths is much more sensitive to local variations in temperature, and so are less suitable for tracing dust column density. Conversely, the 500-|$\mu$|m data are more suitable for tracing dust column density, but have relatively poor angular resolution, at 36 arcsec.

3.3 Automated multicomponent N₂H⁺ line-fitting: mwydyn

Our N₂H⁺ (1–0) observations provide our primary means of assessing the kinematics of the dense gas within our clumps. As discussed earlier, it is necessary to fit the full hyperfine structure for these spectra in order to determine the intrinsic velocity dispersion and velocity centroids, and in many cases we can discern multiple discrete components within the spectra. We have therefore devised a fully automated fitting algorithm – called mwydyn – to decompose these spectra with up to three distinct velocity components that will enable us to carry out our kinematic analysis. In principle, mwydyn is extensible to any molecule with hyperfine structure, and has been developed with testing on ALMA data (Anderson et al., in preparation) as well as the NOEMA data used in this study.

The procedure follows that employed by gildas: class,⁶ which fits an individual N₂H⁺(1–0) spectral component with a model containing four free parameters, but is in principle extensible to any molecule with a comparable hyperfine structure (e.g. HCN and NH₃). The method assumes that each component of the hyperfine multiplet shares the same excitation temperature and linewidth, and that the opacity varies as a function of frequency with a Gaussian profile. The total opacity is the sum of the opacity of the N hyperfine component profiles, which can be expressed as

where τ_total is the sum of the opacities at the individual component line-centres, r_i is the fractional intensity of the ith component (the sum of which is normalized to unity), δv_i is velocity offset of the ith component relative to the velocity of the reference component at v_cen, and Δv is the FWHM linewidth common to all components. The total line profile is then given by

Analytically, T₀ is proportional to τ_total, scaled by a factor that encapsulates the line amplitude.

Running mwydyn on our N₂H⁺ (1–0) cubes results in the fitting of between one and three velocity components to each spectrum, with each fit being described by the four parameters: τ_total, v_cen, Δv, and T₀. We detail the mwydyn algorithm in Appendix B1, but in summary, the algorithm runs by:

• Determining an appropriate noise map.

• Cycling through each individual spectrum (i.e. pixel-by-pixel) whose peak intensity exceeds an SNR threshold. During this process, initial guesses are first generated for a 1-component fit, which is then performed. Next, 2- and 3-component fits are attempted based on the result of the 1-component fit. Finally, the 1-, 2-, and 3- component models are compared, and the best-fit is selected, provided that higher component fits exceed a threshold of improvement over the simpler models.

• Comparing the fits to each spectrum with the fits of neighbouring spectra in order to determine if better solutions have been found locally.

• Write out data products.

Fig. 2 illustrates some of the results from running mwydyn on the N₂H⁺ (1–0) cube for SDC18.888, including a sample of spectra with their fitted components and combined model overlaid, and the output integrated intensity map alongside the residual image. We can see that mwydyn produces a model that matches the data very well, with deviation in the integrated residuals on the order of ∼4 per cent at most. In the case of SDC18.888, emission is detected in every single pixel as a result of the pointing sampling only the centre of the IRDC, and so there is a model fit for every position in the image. The region of the largest integrated residuals is located just to the south-west of the brightest 3.2-mm continuum emission source, indicating features in the spectral lines that were not perfectly modelled by mwydyn, but it should be noted that this region is different to the location of the worst fits as quantified by the Bayesian Information Criterion (BIC) map. The component map also demonstrates that the largest number of components is not necessarily associated with the densest gas as traced by the continuum imaging.

In Fig. 3, we show a 3D illustration of the fitting results to the full data cube for SDC18.888. In colour, the centroid of each component is plotted in (ℓ, b, v_cen) coordinates, where the colour denotes the fitted FWHM of the component, and with a point opacity that is normalized by the integrated intensity of the component. Each surface illustrates a projection of the number of components along the three different axes. It is immediately obvious that the gas in the region is structured and highly complex. We note that high velocity dispersions between the multiple N₂H⁺ (1–0) emission components detected in each spectrum does not necessarily indicate the presence of structures at different physical separations, i.e. coherent structures identified in position-position-velocity (PPV) space do not always map on to coherent structures in 3D space, and that the complexity in PPV space arises naturally in an inhomogeneous turbulent flow (Clarke et al. 2018). We display similar figures for the other six IRDCs in Appendix B2, and we interpret the results in Section 3.4.

mwydyn is written in Python, and is fully parallelized. Run-times on these data were approximately 5 seconds per spectrum per CPU, on a computer cluster dating from 2016. The code is publicly available on GitHub.⁷

3.4 Dense gas tracers

Our observational set-up included a number of molecular species that probe different densities and conditions within our targets: C¹⁸O, HCO⁺, HNC, and N₂H⁺(1–0), and in this section we present an overview of the general picture that they provide. In Fig. 1, we show the integrated intensity (moment 0) maps for these four molecular tracers alongside the continuum images. It is clear that C¹⁸O (1–0) is tracing almost exclusively the material outside of the cores, which have a range of densities that barely overlaps a value of 10⁴  cm⁻³ where we expected CO freeze-out on to dust grains to be significant (Bergin & Tafalla 2007). The map exhibits little correspondence to the highest column densities traced by the 2.8- and 3.2-mm continuum, though the cores are in fact visible in the C¹⁸O (1–0) cubes, but the emission is relatively weak compared with the larger scale emission. We consider this molecule to be a fairly accurate tracer of the clump envelopes.

At the opposite end of the density scale, N₂H⁺ becomes detectable at moderate densities of ≳10⁴  cm⁻³ (Priestley, Clark & Whitworth 2023), and is generally optically thin. The N₂H⁺ (1–0) maps in Fig. 1 are a much better match to the continuum images, but while there clearly is emission that is co-spatial with continuum emission, it is also much more widespread. N₂H⁺ (1–0) is, therefore, tracing both the ambient clump material (at densities above the CO freeze-out) as well as the cores, and is valuable for tracing the transition from clump to core. It is important to recall here that the NOEMA data of all four of the molecular lines imaged here have been complemented with IRAM 30m observations to provide short-spacing information that allows the large-scale emission to be recovered. There are no such complementary observations for the continuum images, and so the differences in the morphology of the emission are caused by a combination of both chemical and observational (i.e. spatial filtering) effects.

HCO⁺ and HNC show similar behaviour and are thought to trace similar conditions (e.g. Barnes et al. 2020; Tafalla, Usero & Hacar 2021), and the maps are very similar, though the HNC maps also resemble a mix of both HCO⁺ and N₂H⁺ (1–0) in several sources (e.g. SDC23.367 and SDC24.433). The critical density for both molecules is similar, with n_crit ∼ 10⁵  cm⁻³ at 10 K, but when accounting for more realistic excitation conditions and photon trapping, their effective excitation density is more like ∼10³  cm⁻³ (Shirley 2015). They are not useful as tracers of column density as a result; however, towards dense cores optically thick tracers like these can be useful as tracers of dynamic processes. For example, HCO⁺ (1–0) often exhibits self-absorbed and blue-asymmetric spectra (e.g. Myers et al. 1996; De Vries & Myers 2005; Fuller, Williams & Sridharan 2005), which are a characteristic signpost of infall under gravitational collapse, and it is indeed useful in the sample for these purposes.

In Fig. 4, we show spectra for each of the targets. In all cases, we show the spectrum at the position of the peak of the 3.2-mm continuum emission (the centre of the MMC in all cases), as well as a secondary spectrum which shows the maximum pixel value across all channels in the field. In the maximum spectra, all targets except SDC24.489 show velocities with emission features that do not correspond to the target clump, indicating that there are other molecular gas structures along the line of sight. This is unsurprising given the intermediate density range probed by C¹⁸O (1–0), and the location of the targets in the inner Galaxy where the spiral structure of the Galaxy is relatively crowded in terms of velocity (e.g. Rigby et al. 2016). We expect CO molecules to become depleted within the dense and cold molecular gas of cores, which may result in a drop in the optical depth that would allow the core systemic velocity to be distinguishable a peak. In general, the peak C¹⁸O (1–0) emission at the position of the cores are well-modelled with Gaussian fits, but we do see such dips in the peak spectra for SDC18.888, SDC24.118, and SDC24.489, which may arise from self-absorption in optically thick gas, coupled with a temperature gradient along the line of sight. SDC24.118 is the clearest example of a complex doubly peaked profile, while SDC24.618 shows the strongest degree of asymmetry.

The picture revealed by the HCO⁺ and HNC (1–0) spectra is more complicated, though generally the two molecules share similar spectral features. We see blue-asymmetric spectra, possibly indicating infall motions, most clearly in SDC23.367, SDC24.489, and SDC24.618, while SDC18.888 shows blue-asymmetric spectra in the most IR-dark positions of the map, but not at the position of the continuum peak. In fact, the HCO⁺ and HNC (1–0) spectra for SDC24.489 show the kind of self-absorbed and blue-asymmetric profiles that are archetypical infall spectra, and it is interesting to note that this is our most IR-dark clump, with the highest f_MMC, and the second strongest HFS as measured by f_C. Large wings in the spectra indicate the presence of outflows in most of these sources, and SDC24.489 is the outlier in this case too, with far less prominent outflow wings in the peak spectrum, though the redshifted side of the spectrum does appear to show some outflow-like features. SDC24.618 is the only target in the sample in which the HNC (1–0) emission differs substantially from HCO⁺ (1–0), and appears to be more closely related to the N₂H⁺ emission.

Our densest gas tracer, N₂H⁺ (1–0) has the most complicated spectra due to its hyperfine structure, and these show wide variation across the sample. The peak spectra in every case require at least two separate velocity components in the mwydyn models indicating that the density structure is complex, even for SDC24.489 whose spectra are the most simple. SDC23.367 and SDC24.630 appear to show signs of self-absorption, which mwydyn tends to model as three separate velocity components, with redshifted and blueshifted components either side of the central weaker component. The three sets of hyperfine lines tend towards a similar brightness temperature when the gas is optically thick, as is seen for these sources. In the case of SDC18.888, the N₂H⁺ (1–0) spectra also seem to be showing outflow wings that match the broad wings seen in HCO⁺ and HNC (1–0).

4. RESULTS

4.1 Core-mass fractions

In this section, we examine several ratios that characterize the relationships between the core populations and their clump-scale environments. The fraction of the 1-pc clump mass that resides within the MMC is

and similarly, the core formation efficiency (CFE) is the fraction of the 1-pc clump mass within its core population:

Finally, we utilize a metric that measures the relative dominance of the MMC over the rest of the core population:

We list these quantities, along with the number of cores per clump, in Table 4. In Fig. 5, we plot the core masses, f_MMC, and CFE – as functions of 1-pc clump mass, f_IRB – indicating the evolutionary status of the clumps, and f_C – indicating the strength of HFS morphology. Both f_MMC and CFE distributions are approximately log-normal, with a mean value of |$\log _{10}(f_\mathrm{MMC}) = -1.03 \approx 9~{{\ \rm per\ cent}}$| and a standard deviation of 0.37 dex, while the CFE distribution has a mean value of |$\log _{10}(CFE) = -0.80 \approx 15.7~{{\ \rm per\ cent}}$| and with a standard deviation of 0.26 dex. We find a strong correlation between the 1-pc clump mass and MMC mass and total mass in cores, and we note that the quoted values for the correlations between the logarithms of these quantities are even stronger, indicating a power-law relationship. We detail the correlation coefficients for these pairs of variables, and for all of the relationships examined in this and subsequent sections in Fig. 10, indicating those that are statistically significant.

In terms of core formation efficiency, we do not see any correlation between either f_MMC or CFE and 1-pc clump mass. In Fig. 5e, we note a lack of f_MMC values greater than 10 per cent for sources with 1-pc clump masses of less than ∼700 M_⊙, which might suggest that high CFE values are not obtainable if the local surface density is too low, though the number statistics here are very small. We see no correlation between the 1-pc clump mass, and the total number of cores per clump.

Examination of the CFE as a function of evolution, as traced by f_IRB, weakly suggests that the largest values are present at earliest times – a trend that is seen both in the f_MMC and CFE. The three sources that go against this trend – SDC326.476, SDC338.315, and SDC340.969 – are from the sample of IR-dark HFSs of Anderson et al. (2021). We find a weak correlation between f_IRB and the number of cores per source, but the trend is not statistically significant given the sample size. Since f_IRB traces only the relative evolution for clumps, and is expected to increase at different rates for clumps of different masses (with denser clumps having shorter free-fall times), it is possible that L_bol/M_1pc may provide a more suitable tracer of evolution, and this quantity is used extensively in the literature for this purpose (e.g. Molinari et al. 2008; Urquhart et al. 2014; Elia et al. 2017). However, we see no correlation between L_bol/M_1pc and any of the core mass-related properties, or the number of cores per clump reported in this section. The L_bol/M_1pc is indeed correlated with T_core within the sample, though the quantities are not always independent, since for cores that are associated with 70-|$\mu$|m compact sources, their temperatures are determined from their 70-|$\mu$|m-derived bolometric luminosities in equation (3).

We see no link between either the filament convergence parameter f_C or the virial parameter α_vir of the sources, and the associated masses or formation efficiencies of their constituent core populations. However, four out of the total of thirteen IRDCs examined are dominated by a single core, which constitutes more than 75 per cent of the total mass in cores within the clump. We find that this quantity, q_MMC shares a negative monotonic relationship with f_C, illustrated in Fig. 6, though this is not statistically significant unless the outlier SDC24.489 is excluded. With a standard dendrogram extraction, such a result could arise as a consequence of cores in a more crowded environment being separated at a higher contour level, artificially lowering their mass. However, we have mitigated this effect with aperture corrections that account for this splitting. We compare the mean densities of the cores to f_C, assuming the cores are spherical:

where the mean molecular weight per hydrogen molecule μ = 2.8. The densities range from |$\sim\!{6} \times 10^4$| to |$\sim\!{4} \times 10^7\, \mathrm{cm}^{-3}$| and are strongly correlated with the convergence parameter, when considering the density of the MMCs, and all cores. This may be an indication of fragmentation into higher density cores in the most strongly convergent systems.

We note that the various mass ratios used in this section (CFE, f_MMC, and most massive core-mass fraction) are distance-independent quantities, meaning that if the dust opacities also share a single value within clumps, that the uncertainties are very small. We present the quantities in Table 4.

4.2 N₂H⁺ (1–0) kinematics revealed by mwydyn

In Section 3.3, we described our method for performing a multiple-velocity component fit to each spectrum of each N₂H⁺ (1–0) cube in our sample. In this section, we discuss the results of this analysis, and how the gas kinematics relates to the positions of the cores seen in the continuum images.

4.2.1 Interpretation of fitted parameters

Before we examine the results of our mwydyn fitting for the seven NOEMA IRDCs, it is important clarify some aspects of the resulting parameters.

One of the key operating concepts of mwydyn is that the emission within each N₂H⁺ (1–0) spectrum can be described as a linear (non-radiatively interacting) combination of 1–3 discrete parcels of gas with the various assumptions of the local thermodynamic equilibrium (LTE)-based fitting procedure (detailed in Section 3.3), such as a single centroid velocity, velocity dispersion, and temperature for each component. In reality, the gas distribution of these regions is continuous across 3D space, and with spatially varying densities over many orders of magnitude which are potentially radiatively coupled. Consequently, mwydyn does not produce unique solutions for each model spectrum, but rather gives a combination of parameters which best fit the data, given the limitations of the model, and the parameter space that it is allowed to explore. We built in an approach where mwydyn will prefer to model a spectrum in the simplest way, with the fewest components, and only add additional velocity components if the fit is significantly improved by them.

mwydyn produces reliable results for spectra which have Gaussian component profiles that are not heavily blended, and this has been tested against synthetic spectra whose parameters can be accurately reproduced. However, real emission spectra are not always so well-behaved, and there are certain features that are present in our N₂H⁺ (1–0) cubes that are not captured by mwydyn. Such complicating features include line-of-sight temperature fluctuations, the presence of outflows, which reveal themselves by extremely broad wings in the profiles, and self-absorption, whose characteristic ‘m’-shaped spectral profiles mwydyn tends to model as two or three distinct velocity components. In all of these cases, the spectral profiles will deviate from Gaussian distributions, and mwydyn will compensate for these by adding extremely optically thick components, that are unlikely to realistically represent the gas properties, but that provide a better statistical fit.

Due to the construction of our observations, which target the central parsec of the IRDCs in question, single-component fits tend to be more common in low-mass clumps, and towards the edges of the ∼1-pc fields (see Fig. C1). 2- and 3-component fits dominate the centres of almost all of our targets, which is an indication that the dense gas kinematics there are more complex, and these often consist of one or more extremely optically thick component (with τ_total ≲ 30); we do not believe these to be accurate measurements of the optical depth of such components, and we caution against interpreting the values of T₀ and τ_total in particular, too literally. Rather, the distribution of these parameters indicate that the gas is showing significant departures from single-temperate LTE conditions with Gaussian profiles. Since the same process is applied to all sources equally, we place the most emphasis on exploring the relative differences between the fitted parameter distributions between, and within, our targets, and it is in this manner that mwydyn is a powerful tool.

4.2.2 Kinematic complexity

The mwydyn fits reveal a generally high-level of what we will refer to here as ‘kinematic complexity’, though ‘dynamical activity’ might be an equally valid phrase. Our single-pointing observations towards these clump centres are not complete mappings of the IRDCs, as is evident in Fig. 1, but are probing the central 1 pc, and consequently we are sampling what are probably the most extreme conditions within those clouds. In Fig. 7, we show distributions of N₂H⁺ (1–0) integrated intensity per pixel for spectra that were fit with 1-, 2-, and 3-component mwydyn fits. This illustrates, firstly, that the fraction of pixels in SDC18.888 that are fit by single components is small relative to the 2- and 3-component fits, and secondly that the higher multiple component fits are found towards the spectra with the greatest integrated intensity. We also find that for the spectra with the greatest integrated intensities in Fig. 7, these tend to be more often composed of 2-component fits as opposed to 3-component fits; this is where the blending of the spectral components has coupled with high-velocity wings that are so extreme that the 3-component fits do not managed to make a substantial improvement, and mwydyn prefers the simpler models. The spectrum at position ‘b’ in Fig. 2 is a good example of this kind of spectrum.

We observe the same trend in all of the sources. In Fig. 9, we show that the mean number of mwydyn-fitted components per spectrum (N_comp) increases as a function of the 1-pc clump mass, and the 2D distributions of this can also be seen in row (f) of Fig. C1. In general, we capture more of the quiescent outer edges of the sources with lower 1-pc clump masses, such as SDC24.118, SDC24.489, SDC24.618, and SDC24.630, for which a border of 1-component fits. For the most massive clumps, SDC18.888 and SDC24.433, we barely probe any quiescent regions, which we have essentially cropped out, resulting in a very small fraction of 1-component fits. Anderson et al. (in preparation) find this same trend of quiescent outskirts and complex interiors in a kinematic follow-up to the six fully mapped IR-dark HFSs of Anderson et al. (2021), with mwydyn fits to the N₂H⁺ (1–0) data.

4.2.3 Distributions of linewidths

To first order, and under the assumption that N₂H⁺ (1–0) is generally optically thin, the combination of the linewidths and the integrated intensities of the fitted components from mwydyn encode the kinetic energy within the gas along a given line of sight. In Fig. 8, we compare the distribution of fitted linewidths from mwydyn for pixels associated with the MMC in each IRDC, and for pixels not associated with any core. We show these distributions using Gaussian kernel-density estimators (KDEs) to estimate the probability density functions. The shaded distributions have been weighted by integrated intensity (a proxy for mass), and we indicate the mean and weighted mean values for each of the distributions with solid and dashed vertical lines, respectively.

In several of the sources, SDC24.630, SDC24.118, and SDC24.618, the distributions that correspond to the MMCs have a greater fraction of spectra in the high-linewidth region of the distribution than the non-core distributions, as indicated by the relative positions of the (weighted and non-weighted) mean values. In some cases, such as SDC24.489 and SDC24.630, and 23.367, the MMC distributions are more sharply peaked at an intermediate linewidth, while the non-core distributions of linewidths are more widely spread, and we note that these are the three most IR-dark clumps with the lowest values f_IRB. For SDC24.433, the MMC and non-core linewidth distributions are very similar, and in the case of SC23.367, the distributions are quite different, with a bimodal distribution for the non-core linewidths, and a skewed but singly peaked distribution for the MMC. It is evident these distributions are systematically skewed to higher FWHMs in the weighted case than in the non-weighted case, indicating that higher column density spectra tend to be more dynamic.

For each source, we test that the FWHM distributions for MMC and non-core components, and that the weighted and un-weighted components are significantly different using a two-sample Anderson–Darling test. The greatest similarity between any of the distributions is between the unweighted MMC and non-core samples of for SDC24.118, and even in this case, the p-value recovered is 0.005, and in all cases, p ≪ 0.05, indicating that the null hypothesis that the two samples are drawn from the same underlying distribution can be rejected. However, the Anderson–Darling tests, with such large samples, become very stringent, and even very slight differences between distributions may result in very low p-values.

There seems to be no single simple relationship that describes the behaviour in MMC-associated and non-core-associated pixels throughout the sample, and their interaction with clump-scale properties such as mass, evolution, and morphology. For example, the two most massive clumps, SDC18.888 and SDC24.433, have significantly higher weighted mean values in the MMCs than in the non-core distributions, but there is no clear mass-related trend amongst the rest of the sample. The position–position–velocity view of the clumps, seen in Figs 3 and B1 show a high level of kinematic complexity. These 3D distributions exhibit many features, such as layers, gradients, sinusoidal oscillations within substructures, cavities, and loops. Such kinematic complexity is not captured by the distributions of linewidths, and so it is not surprising that the linewidth distributions alone are insufficient to elucidate the relationship between linewidths and clump-scale properties.

For four out of seven of these sources, the MMC distributions look similar in shape to the non-core distributions. SDC24.433, SDC24.489, and SDC23.367 have qualitatively differently shaped distributions for the MMC spectra compared with the non-core spectra, but they too do not seem to share any particular set of characteristics, with a diversity of masses, f_IRB, and f_C values (and this extends to the core-related statistics, f_MMC, CFE, and q_MMC). It is possible to come up with reasons for why these sources might be different; SDC24.433 is the most advanced in terms of evolution, with the highest f_IRB and L_bol values, and so is more likely to contain strong stellar feedback that may results in higher linewidths. SDC24.489 is also an outlier in that it is a very strong HFS, and very IR-dark, with an anomalously high q_MMC value (see Fig. 6). However, SDC23.367 does not appear to be remarkable in any particular way. Our sample is probably too small to identify the underlying trends that may be responsible for these similarities and differences. Despite these differences, we suggest that the distributions are not radically different, with similar mean and median values, and spanning a similar range, suggesting that the overall kinematics within the MMCs are not hugely different from the clump interiors, and that the MMCs are not kinematically decoupled from their parent clump at the scales probed by our observations.

4.2.4 Global statistics

To quantify the level of kinematic complexity, a number of further statistics, weighted by the integrated intensity, were calculated. As an initial step, we use an agglomerative clustering algorithm (from sklearn.cluster), to link together all components that are separated by <3 arcsec in the two spatial axes, and 0.5 km s⁻¹ in the spectral axis. We reduce our data to components that reside in the largest cluster identified in this way, in order to exclude any contamination from structures that are probably not spatially associated, but which fall along the same line of sight.

These statistics were calculated for N_comp total number of components within each target. With between 3918 and 4537 spectra per target, we note that the high sensitivity of our observations mean that we have detected and modelled at least one N₂H⁺ emission component in more than 85 per cent of spectra for every source. The total number of modelled components, N_comp, ranges from 6341 to 10 694 across the sample.

We first calculated the weighted mean centroid velocity over all components in each target:

where w_i is the weight of the component, for which we adopted the associated integrated intensity. Similarly, we calculated the weighted mean linewidth:

where |$\sigma _{\mathrm{line}, i} = \mathrm{FWHM}_i/\sqrt{8 \ln {2}}$|⁠), from our mwydyn results. We next calculated the weighted dispersion between the individual velocity components with respect to the weighted mean centroid velocity:

Finally, we calculated the weighted mean total velocity dispersion – a quantity which encapsulates both the spread in centroid velocities about the global (weighted) mean, and the component linewidths:

In Fig. 9, we compare four quantities that summarize the level of ‘kinematic complexity’ in each IRDC: global weighted mean values for the centroid dispersion, linewidths, total velocity dispersion, and number of components with the four quantities that describe the evolutionary state of the host clumps: 1-pc clump mass, f_IRB, f_C, and α_vir. In each case, we have calculated the Pearson correlation coefficients and p-values to test for linear correlations. We further calculated these same properties for only those spectra associated with the MMCs, for which we display the weighted mean values in Fig. 9 too. Correlation coefficients are also presented in Fig. 10. These weighted mean values are generally larger in the MMC spectra compared with the full sample, but this is not always the case.

Considering the statistics for all components (empty markers), we find that the general trend is for three of the four quantities, |$\bar{\sigma }_\mathrm{comp}$|⁠, |$\bar{\sigma }_\mathrm{total}$|⁠, and |$\bar{N}_\mathrm{comp}$| to increase with 1-pc clump mass, while |$\bar{\sigma }_\mathrm{comp}$| and |$\bar{N}_\mathrm{comp}$| are also correlated strongly with convergence parameter. None of the quantities are correlated with f_IRB. There are no correlations between the four mwydyn-based component statistics and the virial parameters within our sample. When considering the difference between the MMC-only statistics and the global statistics, it is interesting that the most of the correlation between 1-pc scale and kinematic properties disappear, while the correlations between M_1pc and |$\bar{\sigma }_\mathrm{line}$| and strengthens marginally.

We have also performed the same calculations on a pixel-wise bases in order to produce 2D maps of these quantities, which we include in Appendix  C, and display in Fig. C1. We discuss our interpretation of these findings in Section 5.

4.3 Relationships between clump-scale, core-scale, and kinematic properties

Thus far, we have examined the relationships between clump- and core-scale properties, and between clump-scale and kinematic properties. Here, we explore the relationships between all of the properties. Since it would be prohibitively tedious to examine individual figures concerning every pair of parameters, in Fig. 10 we present the Pearson and Spearman rank-order correlation coefficients (ρ_P and ρ_S, respectively) between all pairings of the parameters presented thus far. The Pearson coefficients measure the strength of the linear relationship between two variables, while the Spearman’s rank measures the strength of monotonic relationship that is not necessarily linear. We highlight those relationships which are statistically significant, given the p-values of the corresponding tests (p ≤ 0.05), as well as those that are close to that limit (0.05 ≤ p ≤ 0.10) which may be of interest in future studies with larger samples. We note that some of these relationships have stronger and more statistically significant correlations in logarithmic space (resulting in larger ρ_P values) – e.g. log₁₀(M_1pc) versus log₁₀(M_MMC) – but we do not explore these here for the sake of simplicity.

The strongest relationships tend to be between properties of the same type; e.g. core properties such as M_MMC and M_total are most strongly correlated with the other core properties such as f_MMC and CFE. With the exception of the clump-scale properties, this is neither surprising nor interesting because these variables are not independent. However, this illustration allows the identification of new relationships that may be of further interest. For example, sources greater weighted mean linewidths tend to be associated with a greater number of cores.

5. DISCUSSION

5.1 The evolution of clumps

Peretto et al. (2022) identified ∼2000 filaments within the ∼2 deg² GASTON-GPS field, and examined the relationship between their orientation and the physical properties of the ∼1400 clumps within the field. They found that clumps with a higher value of the filament convergence parameter f_C tend to be: (i) more massive and (ii) more IR-bright (as measured by f_IRB) than clumps with a lower f_C value. These suggest that HFSs are either a late-stage configuration in clump evolution, or that the evolution of HFSs is initially rapid, such that examples if IR-dark HFSs are relatively rare. The same survey data were also used to show that clumps accrete mass during the early stages of evolution (Rigby et al. 2021). These results provide an important context for the discussion of the results presented in this study.

On the 1-pc (i.e. ‘clump’) scale, the properties of our sources broadly agree with this picture of clump evolution, with 1-pc clump mass being correlated with convergence parameter, and negatively correlated with virial parameter (albeit both correlations have low statistical significance, with Pearson’s p = 0.19 and 0.30, respectively). This is unsurprising given that five of the seven sources are within the GASTON-GPS field, but this would not be guaranteed to be the case given the dramatically smaller sample size. We do not see a correlation between clump mass and our tracer of evolution, f_IRB, but given the range in clump masses (and therefore free-fall time-scales), f_IRB is not expected to show a strong correlation.

In Section 4.1, we found that both the total mass in cores, and the mass of the MMCs are strongly correlated with the 1-pc clump mass. Anderson et al. (2021) found a similar correlation between MMC mass and total mass in a sample of 35 clumps at various evolutionary stages, and found a much stronger correlation when limiting their sample to six IR-dark HFSs within the sample. Traficante et al. (2023) also report the same correlation in sample of 13 high-mass dense (Σ > 1.0 g cm⁻²) clumps at various stages of evolution. On the other hand, Morii et al. (2023) did not find a correlation between MMC mass and clump mass in the ASHES sample of 3970-|$\mu$|m-dark clumps, which are likely to be at an earlier stage in evolution than the aforementioned studies, though they do find that MMC mass correlates with clump surface density; due to our methodology, this latter result is consistent with the results presented here. When considering the various selection criteria in terms of evolutionary status of the targets, these results may all be consistent in a clump-fed scenario of star formation, where the cores’ growth is promoted by the continuing accretion of material from the wider environment (e.g. Peretto et al. 2020; Rigby et al. 2021), and the correlation between MMC mass and 1-pc clump mass is weak at the earliest stages (i.e. 70-|$\mu$|m-dark stages and 8-|$\mu$|m-dark) and strengthens over time. Indeed, Fig. 5 provides hints that the CFE is relatively high for the least evolved clumps in our sample with low values of f_IRB (i.e. 70-|$\mu$|m-bright but 8-|$\mu$|m-dark), though the sample size is too small to really tell.

We also found that the fraction of the total mass in cores contained within the single MMC decreases as a function of f_C (see Fig. 6). This quantity, q_MMC, may indicate the shape of the CMF with a high value of q_MMC indicating a top-heavy mass function within the clump centres, and vice versa (the quantity is clearly not, however, a robust measure of the complete pre-stellar CMF within each source, as both the spatial extent and resolution are prohibitive in these observations). When coupled with the suggestion that f_C increases over time for clumps, this indicates that mass accretion from the wider environment is initially concentrated on the MMC, with the surrounding cores receiving a greater fraction of infalling material at later stages. This does not necessarily suggest that the MMC stops accreting altogether, though this may well be the case since its early evolution is likely to be the most rapid, and therefore most likely to result in strong stellar feedback (indeed there are no 70-|$\mu$|m-dark pre-stellar MMCs in our sample). Rather, the low q_MMC values in the most hub-filamentary clumps simply suggest that infall or filamentary accretion rates towards the non-MMC cores are higher relative to the MMC at later stages, or this may indicate fragmentation of the MMC at later stages. SDC24.489 is an outlier in our sample for this, whose relatively high value of f_C given its early evolution (with f_IRB = 0) appears to have resulted in an unusually high value of q_MMC – a point which we will return to.

The distributions of N₂H⁺ (1–0) centroid dispersions and, to a lesser extent, linewidths in Fig. 9 also show a tendency to reach higher values for clumps which are more massive and have a higher f_C value. These result in higher total dispersion values, though this trend is driven primarily by the inter-component dispersion per spectrum compared with linewidths. The general kinematic complexity of these regions increases with clump mass and f_C too, indicated by the increasing mean values of the number of mwydyn-fitted components per spectrum. These properties all negatively correlate with virial parameter, too, though the trends are not statistically significant. The relationship with virial parameter is complicated: If the most strongly bound clumps (with low virial parameters) are accompanied by infall signatures in a globally collapsing cloud, the infall itself should increase the linewidth, counteracting the effect on the virial parameter (e.g. Larson 1981; Kauffmann, Pillai & Goldsmith 2013). The weaker correlations are, therefore, expected. The fact that clumps with greatest mass have the lowest virial parameters despite their greater kinematic complexity indicates that gravity continues to dominate the overall balance of energy.

Overall, our measurements would seem to suggest that parsec-scale clump mass (or, equivalently, surface density) is the most important property in determining the dynamic properties of the dense gas in their interiors, and in determining the distribution of core masses. Hub-filament morphology, as measured by f_C, is a secondary factor, and probably one that is also driven by clump mass, with filaments converging more strongly towards more massive clumps over time (Peretto et al. 2022). We do not see any evidence for HFS having properties that separate them as a distinct sample from the non-HFS clumps in our sample, with no visible transition into different behaviour that occurs at f_C = 0.2. The MMCs are located in the most massive clumps, where filament convergence is high, and with a high velocity dispersion between different velocity components of dense gas that requiring a higher level of complexity in mwydyn models. This conclusion is in agreement with Kumar et al. (2020), who suggested that all high-mass stars form in HFSs.

5.2 The dynamic centres of clumps

The gas kinematics, as probed by our four main molecular tracers C¹⁸O, HCO⁺, HNC, and N₂H⁺ (1–0) described in Section 3.4, consistently portray highly dynamic clump interiors. Despite our sample being chosen to cover the earliest IR-dark phases of clump evolution, and an absence in many cases of embedded objects at 8 |$\mu$|m, they all contain compact 70-|$\mu$|m sources, and display evidence for outflows, primarily visible in the HCO⁺ (1–0) spectra, but also often in HNC (1–0) and even, in some cases, N₂H⁺ (1–0). Infall signatures are present in at least three of the seven NOEMA sources, and these tend to be concentrated towards the most 8-|$\mu$|m-dark regions, and are clearest in SDC24.489, our most IRDC and second strongest HFS (with f_IRB  = 0 and f_C  = 0.39). However, an absence of optically thick, self absorbed, and blue-asymmetric spectra does not guarantee the absence of infall. The classical Myers et al. (1996) infall profile assumes spherical symmetry in the density and temperature gradients, while infall occurring under more realistic and chaotic ISM conditions may result in a wide variety of spectral profiles, including redshifted ones, depending viewing angle (Smith et al. 2012), so it is quite possible that there is ongoing infall in the other clumps too. This result is broadly in keeping with Jackson et al. (2019), who found that blue-asymmetric profiles were present in the majority of HCO⁺ (1–0) profiles within a sample of ∼1000 high-mass clumps from the MALT90 survey, and that such profiles were more common at early stages of clump evolution.

The IRDCs in our sample were selected from the sample of Peretto et al. (2023), whose analysis of the velocity dispersion radial profiles on scales of a few tens of parsecs down to a few tenths of a parsec by showed that clumps located within these IRDC centres are dynamically decoupled from their parent molecular clouds, i.e. the radially averaged velocity dispersion as traced by N₂H⁺(1–0) remains constant within the central few parsecs (i.e. the region probed by our observations). The proposed interpretation is that clumps are globally collapsing, while the rest of the molecular gas is not. In that context, the complex kinematics that is being seen in our N₂H⁺ (1–0) data in particular (as well as the other lines) has to be the result of intertwined gravity-driven mass inflows whose density structures vary according to the local conditions.

The dynamic conditions of the dense gas are most clearly visible in the results of our mwydyn fits to the N₂H⁺ (1–0) spectra presented in Section 4.2. In Fig. 8, we compared the distributions of linewidths in the MMCs and in non-core spectra, and they are qualitatively fairly similar, indicating that the MMC kinematics are not radically different from the surrounding gas. The |$\bar{\sigma }_\mathrm{total}$| points in Fig. 9 representing the condensed mwydyn-based total linewidths for the MMC spectra are almost universally located at only marginally higher values than the overall distributions, indicating that the level of kinematic complexity is somewhat elevated in the MMCs. Overall, we suggest that this represents a level of kinematic similarity between the ambient clump material and the gas material, that may be expected if the cores are continually accreting material from the wider environment, as suggested by many other studies (e.g. Peretto et al. 2020; Rigby et al. 2021), with both environments being highly complex, with supersonic linewidths. We temper this conclusion with the caveat that N₂H⁺ (1–0) may not be accurately tracing the kinematics of the most dense gas associated with the cores. Observations of higher density-tracing species such as N₂D⁺ (1–0) may reveal a more accurate picture of the kinematics at core scales. Conversely, the kinematics of the diffuse gas in the clump and core environments (associated with outflows, for example) will not be represented in these distributions.

As mentioned in Section 4.2, the structures revealed in position–position–velocity in Figs 3 and B1 are incredibly diverse. One of the most striking features is the pervasive apparently sinusoidal oscillations that are common across our sample, and which are similar to structures that have been identified elsewhere in the literature. Hacar et al. (2013) identified similar structures (‘fibres’) within C¹⁸O (1–0) observations at similar spatial scales to this study in the nearby filamentary complex L1495/B213 (see also Tafalla & Hacar 2015), and in N₂H⁺ (1–0) in the Orion Integral Shaped Filament (Hacar et al. 2018), and Henshaw et al. (2020) have identified such features across a wide range of scales, from 0.1 pc to kpc within the Milky Way and nearby galaxies. The possible explanation for these features is varied, with Henshaw et al. (2014) finding that they can be reasonably modelled as the result of infall or outflow. These kinds of structures are clearly visible in our data within Figs 3 and B1 in the grey-scale 2D position–velocity projections on the back of each axis, as well as within the coloured data points in the centre. Manipulation of the 3D data suggests that, at least in our IRDC-centre observations, these structures may arise through the projection of concave and convex undulations within sheet-like layer; the oscillatory behaviour visible in 2D projects appear to be representing undulating velocity motions in 3D space that do not have a preferred spatial axis, and it is intriguing to note the lack of apparent filamentary structures (or fibres) on these scales in 3D. The N₂H⁺ (1–0) fibres found in Orion (Hacar et al. 2018) have lengths of ∼0.15 pc and widths of ∼0.03 pc, and so if similar structures were present in our targets, we might reasonably expect to detect (but not resolve) their presence. The close packing of the fibres may cause a level of confusion in data at lower physical resolution, such as ours, so the sheet-like structures in our observations might be a manifestation of this.

Four of the seven NOEMA sources surpass the threshold for being classified as HFSs with f_C > 0.2 (Peretto et al. 2022), and so it is not clear whether we should expect to observe filaments within the central 0.5 pc of such systems. SDC24.618 is our lowest mass clump, with the second-lowest filament convergence parameter (f_C = 0.12) and this is our best candidate for a filamentary or fibrous system. While the general structure in this clump also follows the sheet-like interpretation of the apparent fibres that the eye is drawn to in Fig. B1, where they arise mainly as the projection of convex and concave structures in PPV space, there are a small number of isolated and velocity-coherent structures with large aspect ratios that could be considered to be filaments. If it is generally true that more massive clumps tend to evolve towards a hub-filamentary configuration, then in the very centres of these systems (which our observations target) we might expect the filaments to merge, eradicating the evidence for pre-existing fibres that make those filaments.

5.3 The formation of cores in clump centres

The values of f_MMC and CFE are consistent with being drawn from log-normal distributions centred on 9 and 16 per cent, respectively. The presence of these log-normal distributions indicate that the individual values may be determined by a series of random processes, and we attribute this to the chaotic clump interiors discussed in the previous section, which are generated by conversion of gravitational potential energy to kinetic energy during the global collapse of the clump. In this interpretation, the clumps with extreme values of f_MMC and CFE are simply sampled from the same underlying distribution. The fact that we do not see any clear relationship in Fig. C1 between the locations of extrema in the velocity gradients in our mwydyn fits and the positions, or extrema in the total linewidth of the identified cores, also supports this.

We see tentative evidence for two other aspects in Fig. 5 panels (e) and (f). The first of these aspects is a tendency for the most IR-dark clumps to have the higher values of CFE, which is consistent with Anderson et al. (2021) who also found that clumps at different evolutionary stages (as traced by 8-|$\mu$|m brightness) tended to show different values of f_MMC, with IR-dark sources having high values, and IR-bright sources having lower values. However, the six IR-dark HFSs of that study are also represented in our sample, and so the results are not independent, and our sample size prevents us from drawing a statistically significant result on this point. If the f_MMC and CFE distributions are time variable, then we need a dramatically larger sample to identify this.

Secondly, we point out that none of the clumps with M < 700 M_⊙ – corresponding to a surface density of around Σ_1pc < 0.2 g cm⁻² – achieve |$f_\mathrm{MMC}\ \gt\ 10{{\ \rm per\,cent}}$|⁠. Given the sample size, this is not hugely surprising, but the probability of drawing four consecutive samples that are below 10 per cent is somewhat unlikely, with odds of approximately 1/16. Various studies have proposed thresholds in gas density, above which star formation becomes more efficient. For example, Lada, Lombardi & Alves (2010) proposed a threshold of 116 M_⊙ pc⁻², above which the surface density of star formation correlates linearly with gas surface density, and Heiderman et al. (2010) proposed a similar threshold of 129 M_⊙ pc⁻² (both values correspond to ∼0.025 g cm⁻²). On clump scales, Traficante et al. (2020) found that a sample of 70-|$\mu$|m-dark clumps that become increasingly dominated by non-thermal motions at values of surface density greater than 0.1 g cm⁻². It is possible that what we see here is evidence of a similar threshold, though we note that there are still several clumps with |$f_\mathrm{MMC}\ \lt\ 10{{\ \rm per\,cent}}$| that exceed the threshold. If such a threshold is present, it is not convincing from these data alone, and we require a much larger sample to reach a robust conclusion.

5.4 The role of IR-dark HFSs

Peretto et al. (2022) showed that, on average, clumps hosting HFSs are both more IR-bright, and more massive than non-HFS clumps. Furthermore, Rigby et al. (2021) showed evidence for the accretion of material onto clumps throughout their early- to mid-stages of evolution. In this study, we have also shown that more massive clumps are associated with more massive MMCs (and core populations overall), and greater kinematic complexity in dense gas. All of these results would seem to suggest that high-convergence HFSs are the result of the evolution of clumps, as opposed to a coincidental set of initial conditions that give rise to the growth of the densest and most massive clumps. What, then, is the role of IR-dark hub-systems, such as those of Anderson et al. (2021)? We have identified one source within our sample, SDC24.489, that bucks the trend in Fig. 6, with an extremely large value of q_MMC given its f_C value, compared with the rest of the sample. It is tempting to suggest that, as an extremely strong and IR-dark HFS, SDC24.489 is something of a special source.

However, there are some observational effects that we can expect are having an impact on our interpretation of HFSs, and the use of the convergence parameter f_C. Firstly, much of the focus on HFSs in recent years has depended on targets being identified from within Spitzer 8-|$\mu$|m data as IR-dark features (e.g. Peretto & Fuller 2009; Peretto et al. 2013), which may result in a selection bias, because we cannot locate IR-bright HFSs at similar resolution in blind surveys. In these data, the extinction features allow a view of the gas column density at a resolution of 2 arcsec, but only where the location and status of the cloud are, in some sense, ‘lucky’. Sources which are located at greater distances will be more difficult to see against the diffuse background, as will those which are more evolved (with warmer gas), and we do not have any way of surveying total column density across large areas of the Galaxy at comparable resolution for sources which do not satisfy these selection criteria. Even with the Herschel 250-|$\mu$|m data, the resolution is a factor of nine worse than Spitzer, limiting the visibility of filaments to sources closer than 5 kpc in Kumar et al. (2020). If the 50m AtLAST telescope (Klaassen et al. 2020) is ever built, and equipped with a Band 10 (350-|$\mu$|m) continuum receiver, only then will we have access to large-area and unbiased surveys of column density at these resolutions. For now, the best available data for tracking column density across a wide range of evolutionary stages of clumps appear to be the GASTON-GPS data, but even here the resolution is a factor of six worse than Spitzer. The combination of greater resolution and higher sensitivity may partly explain why evolutionary trends in clump properties had been identified in GASTON-GPS (Rigby et al. 2021) and especially with respect to filament convergence (Peretto et al. 2022) that were not present in ATLASGAL (e.g. Urquhart et al. 2018), though source extraction techniques also differ substantially.

Secondly, the viewing angle surely has an impact on how we study HFSs. If the filaments in HFSs are isotropically arranged in three-dimensions, then it should be expected that some viewing angles should be more favourable than others for their identification. For instance, in a HFS composed of three filaments converging on a central clump with opening angles that are as widely separated as possible, we should expect from some viewing angles to see only two of these three filaments, and we would regard it as a weaker HFS than if we could see all three filaments (and record a lower f_C value). In the extreme case where HFSs may be completely planar (with some indications of this seen in Treviño-Morales et al. 2019 and Anderson et al., in preparation), we may only see their true extent – and measure the highest f_C value – if our viewing angle is such that they are seen face-on. In this case, some fraction of IRDCs that appear as filamentary may turn out to be HFSs at a high angle of inclination. In this case, SDC24.489 may not be an outlier in Fig. 6 at all, but merely have the most fortuitous viewing angle.

To unravel whether SDC24.489 really is an intrinsically extreme source we require an expanded sample that has not been selected from a sample of IRDCs.

6. SUMMARY AND CONCLUSIONS

We have examined the central ∼1 pc of a sample of seven IRDCs with spectral line observations in the 3-mm band using NOEMA and the IRAM 30m telescope, and an expanded sample including a further six IRDCs observed by ALMA in 3-mm continuum. We have developed a Python program called mwydyn, which can decompose N₂H⁺ (1–0) spectra (and, in principle, any molecular transition with hyperfine structure) into up to three distinct velocity components in a fully automated manner, and is a powerful tool for kinematic analysis.

We explore the relationship between the properties of the core populations, and the kinematics traced primarily by N₂H⁺(1–0), and several key properties of their host 1-pc-scale environments which may reveal clues about their formation; the clump mass measured within a 1-pc aperture (equivalent to the local surface density), the infrared-bright fraction (f_IRB) that traces evolution, the filament convergence parameter (f_C) that quantifies the local morphology, and the virial parameter (α_vir), tracing the balance of kinetic and gravitational energy. For the ALMA sources there are no GASTON-GPS-like single-dish data, and so we could only compute f_C values for the sources observed with NOEMA.

Our main conclusions are as follows:

• We identified 47 cores within the central ∼1 pc of our 13 target clumps at a spatial resolution of ∼0.08 pc. The cores range in mass between ∼5 and 1300 M_⊙, with beam-deconvolved radii between ∼0.01 and 0.2 pc.

• The mass of the MMCs, and the total mass in cores are both strongly correlated with the 1-pc clump mass, and occur with an approximately log-normally distributed formation efficiency over the full mass range, with |$f_\mathrm{MMC}\approx 9~{{\ \rm per\ cent}} \pm 0.35$| dex, and CFE |$\approx 16~{{\ \rm per\ cent}} \pm 0.25$| dex, respectively.

• We report tentative evidence that clumps are most efficient at concentrating mass into their MMCs during the earliest and most IR-dark phases.

• We find that the fraction of the total mass in cores that resides within the single MMC is negatively correlated with the filament convergence parameter, implying that HFSs tend to be associated with less dominant MMCs than in less convergent systems.

• Our analysis in C¹⁸O, HCO⁺, HNC, and N₂H⁺ (1–0) consistently reveal a picture of highly dynamic clump centres, which may give rise to the log-normal distributions of f_MMC and CFE.

• The weighted mean values of measures of kinematic complexity of the dense gas traced by N₂H⁺ (1–0), including inter-component velocity dispersion, total linewidth, and number of velocity components, increase with 1-pc clump mass and f_C.

• The distributions of fitted N₂H⁺ (1–0) linewidths in spectra associated with the MMCs are qualitatively similar to the distributions for spectra associated with ambient parsec-scale material. This suggests that the MMCs inherit their kinematic properties – at least in terms of dense gas – from the wider (chaotic) clump-scale environment, and are not kinematically decoupled on ∼0.1-pc scales.

• Filamentary substructure is not found to be a general characteristic of our sample, and we instead find that apparently filament-like oscillatory structures seen in position–velocity projections arise as a consequence of projecting concave and convex structures that appear in a sheet-like morphology in 3D (PPV)-space into 2D (position-velocity)-space. This is most likely a result of observational limitations, though filament merging in the centres of the clumps may also eradicate their signatures.

Our data support a picture of clump evolution in which hub-filamentary morphology is a by-product of the evolution itself, as opposed to being a result of particular initial conditions. The clump’s gravitational potential draws in surrounding filamentary structures over time which, in turn, promote the growth of the MMCs at the convergence point. The MMCs grow most rapidly at early times, with an initially top-heavy CMF. The rate of growth of the neighbouring core masses increases relative to the MMCs at later times, either as a result of accelerating infall, or MMC accretion rates slowing due to protostellar feedback, resulting in an apparent steepening of the CMF. Stellar feedback must ultimately halt the collapse on to the hub centre, with ionizing radiation from young massive stars dispersing the diffuse ambient clump material, leaving the exposed spines of the dense filaments radially arranged around the newly formed central cluster.

We temper these results with the caveat that our sample size is small, with only between seven and thirteen sources (for the kinematic and core-based results, respectively), clump evolution is a multidimensional problem. A much larger sample will be returned to in the future, based on ALMA observations that are ongoing at the time of writing. What we have presented here is a powerful suite of tools – especially mwydyn – that will become increasingly useful as sample sizes expand, enabling the community to untangle many facets of clump evolution and star cluster formation.

ACKNOWLEDGEMENTS

We would like to thank the anonymous referee, whose careful reading and review helped to improve the quality of this manuscript. AJR would like to thank Charlène Lefèvre for assistance with preparing the NOEMA observations, and calibrating, imaging, and cleaning the data. AJR would like to thank the School of Physics and Astronomy at the University of Leeds for a postdoctoral fellowship. AJR and NP acknowledge the support of the Science and Technology Facilities Council consolidated grant number ST/N000706/1 and ST/S00033X/1. GAF also acknowledges support from the University of Cologne and its Global Faculty Program. EJW gratefully acknowledges funding from the German Research Foundation (DFG) in the form of an Emmy Noether Research Group (grant number KR4598/2-1, PI: Kreckel). This work is based on observations carried out under project numbers W18AK and S19AG with the IRAM NOEMA Interferometer and 217–18 and 097–19 with the 30m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). The research leading to these results has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under grant agreement number 730562 [RadioNet].

Following software were used in this research: astropy (The Astropy Collaboration 2022), carta (Comrie et al. 2021), lmfit (Newville et al. 2014), matplotlib (Hunter 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), and sklearn (Pedregosa et al. 2011).

DATA AVAILABILITY

A combined version of Tables 1, 2, 3, and 4, is available in machine-readable format as online supplementary material to this article. The NOEMA continuum and spectral data are publicly available on Zenodo, in a repository accessible at https://doi.org/10.5281/zenodo.10455794. The Spitzer and Herschel data are publicly available. The NIKA maps of SDC18.888 and SDC24.489 were published in Rigby et al. (2018). In keeping with the IRAM Large Program Policy,⁸ the GASTON-GPS data will be made publicly available no later than 2024 October, and a data release paper is being prepared (Zhou et al., in preparation).

Footnotes

References

APPENDIX A: APERTURE CORRECTIONS

The source extraction performed in Section 3.1 used the astrodendro implementation of the dendrogram algorithm (Rosolowsky et al. 2008). The dendrogram algorithm extracts structures based on the topology of the image as a function of contour level. There are several aspects of the methodology which should be considered when using the integrated flux densities recovered by the algorithm, and these are related to two key aspects:

• Firstly, sources are only identified within contours which have a flux density specified by the min_value parameter which, in our case, is equal to a value of three times the rms noise value for each source. Any flux located outside of this contour will not be included in the integrated measurement.

• Secondly, sources that are not perfectly isolated may have their boundaries defined at a higher contour level, and consequently any flux located outside of this contour is not included in the measurement.

In Fig. A1, we present several case studies to illustrate why it is important to account for these effects. The image contains a total of six sources, numbered 1–6 according to the dendrogram extraction that was performed. These models present an idealized noiseless image of sources viewed by a telescope with a 3.8-pixel-FWHM Gaussian beam (similar to our observations). Sources 5 and 6 are point-like, with a peak intensity of 20 and 6, respectively, in arbitrary units. Sources 3 and 4 are extended with respect to the beam, and with peak intensities of 20 and 6, respectively. Sources 1 and 2 are extended, with peak intensities of 20 and 6, respectively, and are partially blended. The bottom three panels show profiles taken through the image at the positions of the dashed white lines, and contain shaded areas to demonstrate the intensity that is accounted for by our integrated flux density measurements. The underreporting of integrated flux densities is worst for point sources at low SNR, and the behaviour can become unpredictable in crowded environments where intensity profiles become blended; in the example here, source 2 acquires a fraction of source 1’s integrated flux density.

This is perfectly appropriate for isolated sources – although it will still underreport the flux as a consequence of point (i) above – but, in regions with either complex background emission, or for blended sources, this may incorporate unrelated flux. By contrast the clipped measurement will more severely underreport the integrated flux of isolated sources, but will do a better job of recovering the appropriate integrated flux of compact sources that are lying in top of a region of extended emission.

The 1D profiles of the sources in Fig. A1 show the fraction of the total integrated flux (in 2D) recovered by the 'clipped' and 'bijected' measurements (as per the scheme of Rosolowsky et al. 2008, printed on the top and bottom, respectively). We define the aperture correction factor, f_ap to be the reciprocal of these values that will, upon multiplication, restore the total flux. The models with the figure demonstrate that the integrated flux densities are more reliably recovered for sources that have a higher peak SNR, and which have larger areas. A further point of note is that, by construction sources 5 and 6 represent pure point sources, with a size equal to the observation beam. The area within the contour defined by min_value contains 21 pixels, compared with the nominal 36 pixels in the beam. If we had set the min_npix parameter to be equal to the number of pixels in the beam, we would not have recovered this source at all, and this motivates our adoption of the half-beam-area for min_npix in Section 3.1.

This far, our models have considered the simplistic case of a Gaussian telescope beam, with no noise. For our observations, we calculate a f_ap that is appropriate for each source, by measuring the equivalent fractions of the total flux in the integrated fluxes within the appropriate synthesized beam. The point spread function of each individual observation is used, and we establish the reference flux within an elliptical aperture that gives the maximum flux up to a radius of two beam radii. For the NOEMA-observed sources, the reference aperture is typically 1.8 beam radii, while for the ALMA-observed ones, this is 2.0 beam radii, with the difference arising from the different observational set-ups. In Fig. A2, we show f_ap as a function of SNR for sources with different areas ranging from 0.5 to 2 beam-areas for SDC18.888.

APPENDIX B: mwydyn

B1 Description of algorithm

The procedure runs as follows:

• First, the data cube is converted from units of surface brightness (Jy beam⁻¹) into main beam brightness temperature (K).

• A noise map is generated by calculating the rms of the first and last 25 channels of each spectrum, which is used to guide the algorithm to only fit the spectra which have a peak SNR of greater than 10.

• All of the valid pixels are cycled through one-by-one. First, a single N₂H⁺(1–0) component fit is performed on the spectrum. The fitting is performed using the Levenberg–Marquardt algorithm-based lmfit Python package (Newville et al. 2014), which requires initial guesses of the four parameters. lmfit offers several advantages over the more widely used scipy.optimize, including providing support for constraints on the parameters.

  For p₁, we simply take the amplitude of the brightest channel as the initial guess, and this channel also provides the initial guess for the centroid velocity p₂. For the linewidth p₃ the initial guess is 0.5 km s⁻¹, which roughly corresponds to the FWHM linewidth corresponding to the isothermal sound speed at 10 K, and for the total opacity p₄, we adopt a value of 0.2 on the basis that N₂H⁺ (1–0) is generally optically thin.

• Next, the algorithm attempts to fit second and third velocity components to the spectrum, comprising of 8- and 12-parameter models, respectively. The initial guesses for the first four parameters (i.e. the parameters for the first component) are made in the same way as in the previous step. For the second set of parameters, a duplicate of the first component initial guesses is used, but with a new centroid velocity guess. This new centroid velocity guess is estimated by comparing the velocity range of detected emission (i.e. the first and last channels which have intensities exceeding an SNR of 3.5), with the velocity range expected for the single-component fit (i.e. ∼15 km s⁻¹, assuming that only its hyperfine multiplets were detected). The value of 3.5 for the minimum SNR for ‘detected’ emission in this case was chosen such that a spurious noise spike would register as a false positive once in every ∼2000 channels, which greatly exceeds the 500 channels of the data cubes in this study.

  For the third component, the same initial guess values are adopted as for the second component.

• Each of the spectra are now cycled through again, in order to determine which of the single-, double-, or triple-component models best fit the data. This is primarily done by comparing the BIC quality of fit estimator for each of the fitted models. The BIC encapsulates the likelihood function, but includes a punishment term for the number of fit parameters such that increasing the number of fit parameters (i.e. adding additional velocity components) does not systematically result in an improved quality of fit. For each spectrum, the model resulting in the lowest BIC is selected as the preferred model, but for the higher order components, we impose that they must constitute an improvement in BIC (i.e. a lower value) of at least 20 over the lower order models.

• Thus far each spectrum has been treated independently of its neighbours, and while lmfit generally provides an excellent exploration of the parameter space, better solutions to similar spectra may have been found in the vicinity. Since the pixel size in our data cubes oversample the beam with ∼4–5 pixels across the major and minor axes, adjacent pixels are strongly correlated. For the final step we cycle through each spectrum in our initial fit, and identify any spectra within a search radius of 2 pixels which have better BIC, and attempt a refit using the corresponding fit parameters as the initial guesses. This process is repeated 10 times to allow any better fitting parameters percolate across the pixel grid, though on these data, the vast majority of spectra take on three or fewer refits. This approach is similar to that of Koch et al. (2021).

• With the fitting now complete, the results are written to a table containing a list of each component in each spectrum, along with its fit parameters, and .fits cubes of the model, and images of the ancillary data products such as BIC and component maps.

There are a number of parameters that control the operation of mwydyn. Firstly, the snrlim parameter sets the minimum peak SNR for a spectrum for mwydyn to attempt to perform a fit. We find that a value of 10 is optimum for this procedure, due to the relative height of the faintest hyperfine component: Below this level, the so-called isolated component begins to reach an SNR of around 3, which has negative implications for our initial guesses for the velocity centroid of the second component as detailed earlier in this section. Secondly, delbic controls the improvement in the BIC of which an additional component must provide in order to be accepted as a better solution. This value is set to 20 by default, which provides a good compromise between including minor new components to fit non-Gaussian line shapes, but still having the flexibility to fit a wide range of line profiles. The refrad parameter controls the number of pixels in the search radius for the spatial refinement stage, and the default value of 2 allows new solutions to be found reasonably rapidly without the increased computational costs of having a higher value. For cleaniter, the number of iterations of spatial refinement, a value of 10 provides enough leeway for new solutions to propagate across the image region. In the case of SDC18.888, this is more than enough, for which no spectrum is re-fitted more than 4 times in this stage.

B2 Fitting results

In Section 3.3, we described mwydyn, a fully automated multiple-velocity-component fitting algorithm for hyperfine spectra, and we displayed the resulting position–position–velocity diagram for SDC18.888. In Fig. B1, we show the same form of diagram with the other six sources in our NOEMA sample.

APPENDIX C: COMPONENT MAPS

In Section 4.2.4, we described several statistics, |$\bar{\sigma }_\mathrm{comp}$|⁠, |$\bar{\sigma }_\mathrm{line}$|⁠, and |$\bar{\sigma }_\mathrm{comp}$|⁠, that measure the overall level of kinematic complexity in each clump by considering the ensemble of mwydyn components as a single distribution. An alternative way to examine these properties is on a per-spectrum basis, which allows the construction of maps. In Fig. C1, we display these maps and overlay the 3.2-mm continuum contours to illustrate where the dust-traced dense gas resides. No consistent pattern between any of these properties and the location of the brightest continuum structures are seen, with the exception of the integrated intensity, which shows a reasonably strong coincidence with the continuum emission. Some of the cores (e.g. SDC18.888, SDC24.489, and SDC24.630) appear to be located in regions where the weighted mean centroid maps have a steep gradient, but this is not always the case. Some of the cores also seem to be located at local minima in the maps of mean linewidth (e.g. SDC23.367, SDC24.118, and SDC24.630), but SDC18.888 presents a counter-example. The MMC in SDC18.888 is associated with a local minimum in the map of the weighted centroid dispersion, but has a high value in the weighted mean linewidth, which suggests that the core resides at a convergence point of components which have large linewidths.