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

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 NOEMA and the IRAM 30-m telescope to trace dense cores with 2.8 mm continuum, and gas kinematics in C$^{18}$O, HCO$^+$, HNC, and N$_2$H$^+$ ($J$=1$-$0). We supplement our continuum sample with six IRDCs observed at 2.9 mm with 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$_2$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% - 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.


INTRODUCTION
The stellar populations of the many billions of galaxies in the Universe are dominated, in absolute number, by low-mass stars ( * ≲ 2 M ⊙ ), while the much rarer high-mass stars ( * > 8 M ⊙ ) have an enormous influence on the interstellar medium (ISM) and chemical evolution of galaxies.The relative abundance of lowthrough 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 ★ E-mail: a.j.rigby@leeds.ac.uk (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;Polychroni et al. 2013;André et al. 2014;Könyves et al. 2020;Ladjelate et al. 2020).However, for high-mass stars, the picture is more complicated (see Motte et al. 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-2Myr, Battersby et al. 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 find 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 (∼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 of ∼ 10 −4 -10 −3 M ⊙ yr −1 (e.g.Peretto et al. 2014;Kirk et al. 2013;Williams et al. 2018;Lu 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 IRDCs in mid-IR survey data (Peretto & Fuller 2009) which have superior resolution offered compared to sub-millimetre or millimetre survey data; for example, the 8 μm GLIMPSE survey data (Benjamin et al. 2003) have a resolution of 2 arcseconds compared to 18 arcseconds 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 μ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 (  IRB ; Rigby et al. 2021; see also Watkins et al. in prep.), 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  IRB metric was used by Peretto et al. (2022) in conjunction with another new metric -the filament convergence parameter,  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  C tends to increase with  IRB , suggesting that either HFSs are a late stage of clump evolution, or that infrared-dark HFSs have short lifetimes compared to 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 hub-filament systems.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 infrareddark HFSs to the cores within a larger sample of clumps mapped with ALMA, and found that infrared-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 most massive core, than for infrared-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 program (Motte et al. 2022), Nony et al. (2023) are finding a time dependence in the evolution of core mass functions (CMFs) in W43 (Pouteau et al. 2023), and Nony et al. (2023) found that the protostellar core mass function (CMF) in W43 is more top-heavy than the Salpeter-like pre-stellar CMF, indicating that the most massive cores 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 most massive cores, 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 2 H + across a sample of IR-dark clumps with a range of properties.
The fundamental rotational transition of the molecular ion diazenylium, N 2 H + , at 93 GHz was first detected by Turner (1974), and later identified by Green et al. (1974), and confirmed by Thaddeus & Turner (1975) and has since been widely utilised 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 transitionsincluding the =1-0 transition at 93.174 GHz that we have observed in this study -have at least seven hyperfine components (Caselli et al. 1995), though they often appear as a triplet in the interstellar medium (ISM) where linewidths are typically supersonic.The high-frequency  1 ,  = 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 −1 , and a low-frequency triplet whose component centroids extend to ∼ +7 km s −1 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 can not 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 2 H + ; however, the high-spatial resolutions that are now achievable with facilities like 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 compared to fitting the full spectrum.To overcome these problems, we have developed a fullyautomated multiple velocity component hyperfine line-fitting code to analyse N 2 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 2 H + (1-0) alongside several other key molecules and continuum, in order to examine the kinematics surrounding the most massive cores, 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 coreand clump-scale properties.Our main results are described in Section 4, and discussed in Section 5, and we present our conclusions in Section 6.

Sample selection
We compiled a sample of seven IRDCs from the catalogue of Peretto & Fuller (2009) to examine with NOEMA, which are a sub-sample 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 < /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 programs for these seven IRDCs.Systemic velocities for each of the IRDCs in our sample had been determined by observations of N 2 H + (1-0) in Peretto et al. (2023) using EMIR on the IRAM 30-m 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) 1 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 (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 characterised in Section 3.2.

Observations
The column density peaks of the seven IRDCs -which we refer to herafter as the central clumps -were observed with single NOEMA pointings between 21st May and 30th August 2019 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 GHz 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 −1 channels for our main target lines which were the (1-0) transitions of HCO + , HNC, N 2 H + and C 18 O.
A total of nine 15-m dishes were used in the 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 cent2 in the 3 mm band.Short-spacing observations for our spectral lines were obtained using the IRAM 30-m telescope with dedicated observing runs in November 2019, and February, May, July, and October 2020.For the short-spacing observations we used EMIR with a setup optimised for our main target lines, with the local oscillator tuned to 99.5 GHz, and the FTS backend 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 μm (Benjamin et al. 2003), Herschel/Hi-GAL 70 μm (Molinari et al. 2016), and IRAM 30-m/NIKA and GASTON 1.2 mm (Rigby et al. 2018(Rigby et al. , 2021) ) continuum images showing the wider environment.

Data reduction
Data reduction was performed using software from the gildas suite3 .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 mm 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 -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 arcseconds for the N 2 H + (1-0) and 3.2 mm continuum observations, for which the synthesised 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.

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 arcseconds) 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 Hi-GAL 500 μm observation (Molinari et al. 2016) of the same cloud.This latter choice was made Table 1.Summary of observations, with columns as follows: 1) Source designations identify them in the catalogue of Peretto & Fuller (2009); 2-3) coordinates of the NOEMA pointings; 4) systemic velocities from Peretto et al. (2023); 5-6) heliocentric distance and uncertainty; 7-8) major and minor axes of the synthesised beam, and rms sensitivity of the 2.8 mm continuum; 9-10) major and minor axes and rms sensitivity of the N 2 H + (1-0) observations; 11) total clump mass contained within a 1-pc-diameter aperture; 12) bolometric luminosity of compact sources within a 1-pc aperture; 13) infrared-bright fraction within the 1-pc aperture; 14) maximum filament convergence parameter within the 1-pc aperture.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 characterisitics of the resulting processed data, including beam sizes, can be found in Table 1.

Ancillary data
We make use of 8 μm continuum imaging data from the Spitzer Galactic Legacy Infrared Midplane Survey Extraordinaire (GLIMPSE; Churchwell et al. 2009), which have an angular resolution of 2 arcseconds.We also make use of 70, 160, and 250 m imaging from the Herschel infrared Galactic Plane Survey (Hi-GAL Molinari et al. 2016), and column density maps generated from these images from Peretto et al. (2016).The 70, 160, and 250 μm data have angular resolutions of ∼8, 12, and 18-arcsec, respectively, and the column density maps have the same resolution of the 250 μ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.489from Rigby et al. (2018).Both the NIKA (Monfardini et al. 2010) and NIKA2 data have an angular resolution of 11 arcsec.

Core population
The population of compact sources, within each field in the 2.8 mm continuum data was extracted using astrodendro4 , 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 signal-to-noise ratio, which are detailed in Appendix A. We also find that the inclusion of this minimum source size assists with breaking down over-sized 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, signal-tonoise ratio, 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   is the flux density integrated over the source area,  is the distance to the source,   ( d ) is the Planck function evaluated at frequency  and dust temperature  d , and the dust opacity is given by:   = 0.1 cm 2 g −1 (/THz)  .This value incorporates a gas-to-dust mass ratio of 100, and is the typical value adopted by e.g.Marsh et al. (2015), Peretto et al. (2016), andAnderson 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-μm colour, again assuming  = 2.This technique is limited to the resolution of the 250-μm data, which is 18 arcsec, and so we may overestimate the temperature of cold cores compared to the ambient dust temperature on larger scales.For cores containing any embedded star formation, we estimate the core temperature from the 70-μ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-μ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  70 is the integrated flux density of the associated 70-μm compact source.For optically thin dust heated by an internal star, the temperature at a particular radius can then be determined (Terebey et al. 1993): where  = 2/(4+ ),  = 2, and  0 = 25 K is the dust temperature at a reference radius  0 = 0.032 pc for a star of luminosity  0 = 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),  =  eff = √︁ /, deconvolved by the NOEMA beam.In this way, we determined the masses of a total of 47 cores, with 21 cores 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.579MM1, a well-known high-mass core (Peretto et al. 2013).With 10% and 5% 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  0 (Ossenkopf & Henning 1994).The core masses are therefore thought to be accurate to within a factor of ∼2.

Clump-scale properties
We calculated the masses of the central clumps within the thirteen IRDCs from the Peretto et al. (2016) column density maps5 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 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' ( 1pc ), as a reminder of our fixed-aperture calculation.Because this technique adopts a fixed aperture size, the measurement is equivalent to the number of blank pixels near the column density peaks due to saturation in either the 160 or 250 μm Herschel bands.These were filled by interpolation using the interpolate_replace_nans function from astropy.convolution with a single-pixel standard deviation Gaussian kernel.average surface density.1-pc clump masses are calculated according to: where the surface area element d =  2 dΩ, for the source distance  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: We adopt a value of 2.8 for the mean molecular weight per hydrogen molecule,  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-4360 M ⊙ , corresponding to mean surface densities of Σ pc ∼0.1-1.2 g cm −2 across the aperture.For context, this range of densities is encompassed by the density range of 0.05-1.0g cm −2 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, 2 k =  G .Following the formulation of Bertoldi & McKee (1992), we calculate the virial parameter: where  is the total linewidth (including thermal and non-thermal contributions),  is the radius, and  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 2 H + (1-0) integrated intensity.We adopted the radius of a circle with the same area  enclosed by the contour,  ctr = √︁ /.We then calculated the background-subtracted column density within this same contour from the Herschel-derived column density maps as in Eq. 5, and term this  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 2 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 () ∝  −2 , equipartition of kinetic and gravitational energy occurs at  vir = 3.3.For the seven NOEMA clumps in our sample for which we have N 2 H + (1-0) linewidths, virial parameters are in the range 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.
Table 2. Calculated properties of the cores identified in Section 3.1: designations ordered by integrated flux density for each target; integrated flux density at 2.8 mm and 2.9 mm for the NOEMA-and ALMA-observed sources, respectively; the aperture correction that has been applied to the integrated flux density; mass, with uncertainties corresponding to the 16th-84th percentiles of the Monte Carlo-sampled uncertainty distributions; beam-deconvolved equivalent radius; core temperature.The bolometric luminosity of each clump is calculated using Equation 2, and using the total integrated flux density of all 70-m compact sources from Molinari et al. (2016) that lie within the 1-pc aperture.Bolometric luminosities calculated in this way vary from ∼15-50000 L ⊙ across the sample, with all but two of the thirteen sources having  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.

Designation
We also supplement these properties with two new quantifications of each clump's evolutionary status.Firstly, for all thirteen sources we calculate the infrared-bright fraction  IRB , which determines the fraction of pixels within the clump at 8 μm (from Spitzer/GLIMPSE imaging; Churchwell et al. 2009) that are brighter than the local background within a 4.8-arcminute-wide box, and which has shown to be a good tracer of relative evolution (see Rigby et al. 2021, Watkins et al. in prep.).Sources evolve from  IRB = 0 at the earliest stages of star formation to  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  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 arcseconds, 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 (, ) -as: where  fil and  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 (, ) to pixel , and the filament direction at pixel .  is a normalisation constant which ensures that a pixel at the convergence point of six parsec-long radially distributed filaments has a value of  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  C → 1 for sources which are located at the converging point of several nearby filaments, or a low value of  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  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 μ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 μm data are more suitable for tracing dust column density, but have relatively poor angular resolution, at 36 arcsec.

Automated multi-component N 2 H + line-fitting: mwydyn
Our N 2 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 3 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 prep.) as well as the NOEMA data used in this study.
The procedure follows that employed by gildas: class6 , which fits an individual N 2 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, NH 3 ).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  hyperfine component profiles, which can be expressed as: where  total is the sum of the opacities at the individual component line-centres,   is the fractional intensity of the  th component (the sum of which is normalised to unity),   is velocity offset of the  th component relative to the velocity of the reference component at  cen , and Δ is the FWHM linewidth common to all components.The total line profile is then given by: Analytically,  0 is proportional to  total , scaled by a factor that encapsulates the line amplitude.Running mwydyn on our N 2 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 ,  cen , Δ, and  0 .We detail the mwydyn algorithm in Appendix B1, but in summary, the algorithm runs by: (i) Determining an appropriate noise map.(ii) Cycling through each individual spectrum (i.e.pixel-by-pixel) whose peak intensity exceeds a signal-to-noise ratio 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.(iii) Comparing the fits to each spectrum with the fits of neighbouring spectra in order to determine if better solutions have been found locally.
(iv) Write out data products.
Figure 2 illustrates some of the results from running mwydyn on the N 2 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% 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 are not necessarily associated with the densest gas as traced by the continuum imaging.
In Fig. 3, we show a 3-D illustration of the fitting results to the full data cube for SDC18.888.In colour, the centroid of each component is plotted in (ℓ, ,  cen ) coordinates, where the colour denotes the fitted FWHM of the component, and with a point opacity that is normalised 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 2 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 onto coherent structures in 3-D 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 parallelised.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 GitHub7 .

Dense gas tracers
Our observational setup included a number of molecular species that probe different densities and conditions within our targets: C 18 O, HCO + , HNC, and N 2 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 18 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 4 cm −3 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 18 O (1-0) cubes, but the emission is relatively weak compared to 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 2 H + becomes detectable at moderate densities of ≳ 10 4 cm −3 (Priestley et al. 2023), and is generally optically thin.The N 2 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 with continuum emission, it is also much more widespread.N 2 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 30-m 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 et al. 2021), and the maps are very similar, though the HNC maps also resemble a mix of both HCO + and N 2 H + (1-0) in several sources (e.g.SDC23.367 and SDC24.433).The critical density for both molecules is similar, with  crit ∼ 10 5 cm −3 at 10 K, but when accounting for more realistic excitation conditions and photon trapping, their effective excitation density is more like ∼10 3 cm −3 (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 et al. 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 most massive core 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.489show 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 18 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 18 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.618shows 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.888shows blue asymmetric spectra in the most infrared-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.489show 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 infrared-dark clump, with the highest  MMC , and the second strongest hub-filament system as measured by  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 red-shifted side of the spectrum does appear to show some outflowlike 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 2 H + emission.
Our densest gas tracer, N 2 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.489whose spectra are the most simple.SDC23.367 and SDC24.630appear  .118are shown in bold, while the thin spectra show the maximum value across all channels in the field.In the case of the C 18 O spectra, a single Gaussian fit has been overlaid, and for the N 2 H + (1-0) spectra, we show the fits obtained with mwydyn, where the parameter designation is given by the formalism in Section 3.3, where p1, p2, p3, and p4 refer to  0 ,  cen , Δ, and  total , respectively.The systemic velocities identified by the single-dish data are shown as dashed vertical lines.
to show signs of self-absorption, which mwydyn tends to model as three separate velocity components, with red-shifted and blue-shifted 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 2 H + (1-0) spectra also seem to be showing outflow wings that match the broad wings seen in HCO + and HNC (1-0).

Core mass fractions
In this Section, we examine several ratios that characterise the relationships between the core populations and their clump-scale environments.The fraction of the 1-pc clump mass that resides within the most massive core is: and similarly, the core formation efficiency (CFE) is the fraction of the 1-pc clump mass within its core population: Finally, we utilise 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 Tab. 4. In Fig. 5, we plot the core masses,  MMC , and CFE -as functions of 1-pc clump mass,  IRB -indicating the evolutionary status of the clumps, and  C -indicating the strength of HFS morphology.
Both  MMC and CFE distributions are approximately log-normal, with a mean value of log 10 (  MMC ) = −1.03≈ 9% and a standard deviation of 0.37 dex , while the CFE distribution has a mean value of log 10 () = −0.80≈ 15.7% 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  MMC or CFE and 1-pc clump mass.In Fig. 5e, we note a lack of  MMC values greater than 10 per cent for sources with -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 Figure 5.The top row shows 2.8 mm core masses (solid symbols are MMCs, empty symbols are the sum of the mass in all cores) as functions of: a) 1 pc mass (and its equivalent surface density on the secondary -axis); b) the infrared-bright fraction; c) filament convergence parameter; d) virial parameter.The bottom row shows core formation efficiencies as functions of the same four parameters (solid indicating  MMC , empty symbols for CFE).The colour coding shows the core populations for each IRDC, with IRDCs from the NOEMA and ALMA samples marked with circles and squares, respectively.The plots showing the filament convergence parameter have been shaded for values of  c > 0.2 which signifies the parameter space containing hub-filament systems (Peretto et al. 2022).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  IRB , weakly suggests that the largest values are present at earliest times -a trend that is seen both in the  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 hub-filament systems of Anderson et al. (2021).We find a weak correlation between  IRB and the number of cores per source, but the trend is not statistically significant given the sample size.Since  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  bol / 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  bol / 1pc and any of the core mass-related properties, or the number of cores per clump reported in this Section.The  bol / 1pc is indeed correlated with  core within the sample, though the quantities are not always independent, since for cores that are associated with 70 μm compact sources, their temperatures are determined from their 70 μm-derived bolometric luminosities in Eq. 3.
We see no link between either the filament convergence parameter  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 13 IRDCs examined are dominated by a single core, which constitutes more than 75% of the total mass in cores within the clump.We find that this quantity,  MMC shares a negative monotonic relationship with  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  C , assuming the cores are spherical: where the mean molecular weight per hydrogen molecule  = 2.8.The densities range from ∼ 6 × 10 4 cm −3 to ∼ 4 × 10 7 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,  MMC , 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.

N 2 H + (1-0) kinematics revealed by mwydyn
In Sec.3.3 we described our method for performing a multiplevelocity component fit to each spectrum of each N 2 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.

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 2 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 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 three-dimensional space, and with spatially varying densities over many orders of magnitude which are potentially radiatively coupled.Consequently, mwy-dyn 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 wellbehaved, and there are certain features that are present in our N 2 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 ∼1pc 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  0 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.

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 2 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 ( comp ) increases as a function of the 1-pc clump mass, and the 2-dimensional 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 1component fits.Anderson et al. (in prep.)find this same trend of quiescent outskirts and complex interiors in a kinematic follow up to the six fully-mapped infrared-dark HFSs of Anderson et al. (2021), with mwydyn fits to the N 2 H + (1-0) data.

Distributions of linewidths
To first order, and under the assumption that N 2 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 most massive core in each IRDC, and for pixels not associated with any core.We show these distributions using Gaussian kernel-density estimates (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.489and 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 infrareddark clumps with the lowest values  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 singlypeaked 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 -value recovered is 0.005, and in all cases,  ≪ 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 -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 noncore 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 sub-structures, 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.367have qualitatively differently-shaped distributions for the MMC spectra compared to the non-core spectra, but they too do not seem to share any particular set of characteristics, with a diversity of masses,  IRB , and  C values (and this extends to the core-related statistics,  MMC , CFE, and  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  IRB and  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 hub-filament system, and very infrared-dark, with an anomalously high  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.

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 less than 3 arcsec in the two spatial axes, and 0.5 km s −1 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  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 2 H + emission component in more than 85% of spectra for every source.The total number of modelled components,  comp , ranges from 6341 to 10694 across the sample.
We first calculated the weighted mean centroid velocity over all components in each target: where   is the weight of the component, for which we adopted the associated integrated intensity.Similarly, we calculated the weighted mean linewidth: where  line, = FWHM  / √ 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 dispersiona 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 summarise 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,  IRB ,  C , and  vir .In each case we have calculated the Pearson correlation coefficients and -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 to 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, σcomp , σtotal , and Ncomp to increase with 1-pc clump mass, while σcomp and Ncomp are also correlated strongly with convergence parameter.None of the quantities are correlated with  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  1pc and σline and strengthens marginally.We have also performed the same calculations on a pixel-wise bases in order to produce 2-dimensional maps of these quantities, which we include in Appendix C, and display in Fig. C1.We discuss our interpretation of these findings in Sec. 5.

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 -values of the corresponding tests ( ≤ 0.05), as well as those that are close to that limit (0.05 ≤  ≤ 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)for example log 10 ( 1pc ) vs. log 10 ( MMC ) -but we do not explore these here for the sake of simplicity.

The evolution of clumps
Peretto et al. ( 2022) identified ∼2000 filaments within the ∼ 2 deg 2 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  C tend to be: i) more massive and ii) more infrared-bright (as measured by  IRB ) than clumps with a lower  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  = 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,  IRB , but given the range in clump masses (and therefore free-fall timescales),  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 −2 ) 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 39 70 μmdark 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 clumpfed 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-μm-dark stages and 8-μ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  IRB (i.e.70-μm-bright but 8-μ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 most-massive core decreases as a function of  C (see Fig. 6).This quantity,  MMC , may indicate the shape of the core mass function with a high value of  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 core mass function within each source, as both the spatial extent and resolution are prohibitive in these observations).When coupled with the suggestion that  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-μm-dark pre-stellar MMCs in our sample).Rather, the low  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  C given its early evolution (with  IRB = 0) appears to have resulted in an unusually high value of  MMC -a point which we will return to.
The distributions of N 2 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  C value.These result in higher total dispersion values, though this trend is driven primarily by the inter-component dispersion per spectrum compared to linewidths.The general kinematic complexity of these regions increases with clump mass and  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 et al. 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 parsecscale 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  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  C = 0.2.The most massive cores 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.

The dynamic centres of clumps
The gas kinematics, as probed by our four main molecular tracers C 18 O, HCO + , HNC, and N 2 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 μm, they all contain compact 70 μ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 2 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 μm-dark regions, and are clearest in SDC24.489,our most infrared-dark cloud and second strongest hub-filament system (with  IRB = 0 and  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 red-shifted 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 2 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 2 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 2 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 points σ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 most massive cores.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 2 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 2 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 positionposition-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 18 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 2 H + (1-0) in the Orion Integral Shaped Filament (Hacar et al. 2018), andHenshaw 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 grayscale two-dimensional position-velocity projections on the back of each axis, as well as within the coloured data points in the centre.Manipulation of the three-dimensional 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 2-dimensional projects appear to be representing undulating velocity motions in 3-dimensional 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 3-D.The N 2 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 hub-filament systems with  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 lowestmass clump, with the second-lowest filament convergence parameter (  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.

The formation of cores in clump centres
The values of  MMC and CFE are consistent with being drawn from log-normal distributions centred on 9 and 16 percent, 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  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 IRdark 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 μm brightness) tended to show different values of  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  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  < 700 M ⊙ -corresponding to a surface density of around Σ 1pc < 0.2 g cm −2 -achieve  MMC > 10%.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 et al. (2010) proposed a threshold of 116 M ⊙ pc −2 , 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 −2 (both values correspond to ∼ 0.025 g cm −2 ).On clump scales, Traficante et al. (2020) found that a sample of 70-μm-dark clumps that become increasingly dominated by non-thermal motions at values of surface density greater than 0.1 g cm −2 .It is possible that what we see here is evidence of a similar threshold, though we note that there are still several clumps with  MMC < 10% 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 hub-filament systems Peretto et al. (2022) showed that, on average, clumps hosting HFSs are both more infrared-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 hub-filament systems 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 infrared-dark hubsystems, 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  MMC given its  C value, compared to the rest of the sample.It is tempting to suggest that, as an extremely strong and infrared-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 hub-filament systems, and the use of the convergence parameter  C .Firstly, much of the focus on HFSs in recent years has depended on targets being identified from within Spitzer 8-μm data as infrared-dark features (e.g.Peretto & Fuller 2009;Peretto et al. 2013), which may result in a selection bias, because we can not locate infrared-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 arcseconds, 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-μ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 50-m At-LAST telescope (Klaassen et al. 2020) is ever built, and equipped with a Band 10 (350-μ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 threedimensions, 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  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 prep), we may only see their true extent -and measure the highest  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.489really is an intrinsically extreme source we require an expanded sample that has not been selected from a sample of IRDCs.

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 30-m 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 2 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 2 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 (  IRB ) that traces evolution, the filament convergence parameter (  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  C values for the sources observed with NOEMA.
Our main conclusions are as follows: (i) 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.
(ii) The mass of the most massive cores, 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  MMC ≈ 9% ± 0.35 dex, and CFE ≈ 16% ± 0.25 dex, respectively.
(iii) We report tentative evidence that clumps are most efficient at concentrating mass into their most massive cores during the earliest and most infrared-dark phases.
(iv) We find that the fraction of the total mass in cores that resides within the single most massive core is negatively correlated with the filament convergence parameter, implying that HFSs tend to be associated with less dominant most-massive cores than in less convergent systems.
(v) Our analysis in C 18 O, HCO + , HNC, and N 2 H + (1-0) consistently reveal a picture of highly dynamic clump centres, which may give rise to the log-normal distributions of  MMC and CFE.
(vi) The weighted mean values of measures of kinematic complexity of the dense gas traced by N 2 H + (1-0), including intercomponent velocity dispersion, total linewidth, and number of velocity components, increase with 1-pc clump mass and  C .
(vii) The distributions of fitted N 2 H + (1-0) linewidths in spectra associated with the most massive cores are qualitatively similar to the distributions for spectra associated with ambient parsec-scale material.This suggests that the MMCs inherit their kinematic propertiesat least in terms of dense gas -from the wider (chaotic) clump-scale environment, and are not kinematically decoupled on ∼0.1-pc scales.
(viii) Filamentary substructure is not found to be a general characteristic of our sample, and we instead find that apparently filamentlike 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 3-D (PPV)-space into 2-D (position velocity; PV)-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 hubfilamentary 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 core mass function.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 core mass function.Stellar feedback must ultimately halt the collapse onto the hub centre, with ionising 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 multi-dimensional 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.

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: (i) 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.
(ii) 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 idealised 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 demonstrated the intensity that is accounted for by our integrated flux density measurements.The under-reporting of integrated flux densities is worst for point sources at low signal-tonoise ratio, 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 under-report 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 under-report 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 1-D profiles of the sources in Fig. A1 show the fraction of the total integrated flux (in 2-D) recovered by the clipped and bĳected measurements (printed on the top and bottom, respectively').We define the aperture correction factor,  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 signalto-noise ratio, 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 to 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  ap that is appropriate for each source, by measuring the equivalent fractions of the total flux in the integrated fluxes within the appropriate synthesised beam.The PSF 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 setups.In Fig. A2, we show  ap as a function of signal-to-noise ratio for sources with different areas ranging from 0.5 to 2 beam-areas for SDC18.888. is 0.5 km s −1 , which roughly corresponds to the FWHM linewidth corresponding to the isothermal sound speed at 10 K, and for the total opacity  4 , we adopt a value of 0.2 on the basis that N 2 H + (1-0) is generally optically thin.
(iv) 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 S/N of 3.5), with the velocity range expected for the single-component fit (i.e.∼15 km s −1 , assuming that only its hyperfine multiplets were detected.The value of 3.5 for the minimum S/N 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.
(v) 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 Bayesian Information Criterion (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.
(vi) 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 3 or fewer refits.This approach is similar to that of Koch et al. (2021).
(vii) 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 .fitscubes 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 signal-to-noise ratio 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 a signal-to-noise ratio 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 Sec.3.3 we described mwydyn, a fully-automated multiplevelocity-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 Sec.4.2.4,we described several statistics, σcomp , σline , and σ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,SDC24.630), but SDC18.888presents a counter-example.The most massive core 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.

Figure 1 .
Figure 1.Each row shows all seven IRDCs in different bands.From top to bottom, the rows show: 1) Spitzer/GLIMPSE 8-μm continuum (Churchwell et al. 2009); 2) Herschel/Hi-GAL 70-μm continuum (Molinari et al. 2016); 3) 1.2 mm continuum from NIKA for SDCs 18.888 and 24.489 (Rigby et al. 2018), otherwise from NIKA2 (Rigby et al. 2021); 4) NOEMA 107 GHz / 2.8 mm continuum emission; 5) NOEMA 93 GHz / 3.2 mm continuum emission; 6) NOEMA + IRAM 30-m C 18 O (1-0) integrated intensity; 7) NOEMA + IRAM 30-m HCO + (1-0) integrated intensity; 8) NOEMA + IRAM 30-m HNC (1-0) integrated intensity; 9) NOEMA + IRAM 30-m N 2 H + (1-0) integrated intensity.The beam size is shown in the lower-left corner of each image, and a scalebar representing a distance of 0.5 pc is show in the lower-right corner.All coordinates are given with reference to the central coordinates given in Table 1.The top three rows (showing single-dish data) have a field width a factor of 2 larger than for the bottom six rows to provide more environmental context, with the smaller fields of view displayed as dashed boxes.Contours are shown for the 2.8 mm continuum for each image at levels of 10 and 50 times the rms value.

Figure 2 :
Figure 2: An illustration of the resulting fits from mwydyn for N 2 H + (1-0) in SDC18.888NOEMA.The top-left panel shows the integrated intensity, and is marked with 3 crosses corresponding to spectra shown in the top-right panel.The observed spectra are shown as steps, and the best-fitting model is shown as the solid line, along with its constituent components in dashed, dash-dotted, and dotted lines with centroid velocities indicated by arrows.On the bottom row, the left panel shows the model integrated intensity, with the residual image (data−model) shown in the middle-left panel.The middle right panel is a map of the BIC, and the right panel shows the number of velocity components for each spectrum.The black contours in all images show the 3.2 mm continuum flux density, starting from 1 mJy beam −1 , with steps of 2 mJy beam −1 .

Figure 3 :
Figure3: A three-dimensional representation of the fittingresults from mwydyn for SDC18.888, in which the centroid velocity for each component is plotted against its Galactic coordinates.The points are coloured according to the FWHM linewidth of the component, with a transparency normalised by the integrated intensity of the line.Projections of the fitting results are also shown along each of the back surfaces as grayscale, in which the point colour is determined by the number of components along that axis.The lower surface shows the number of components fitted to each spectrum in grayscale, accordingly.

Figure 4 .
Figure 4. Spectra of C 18 O (1-0), HCO + (1-0), HNC (1-0) and N 2 H + (1-0) emission at the position of the brightest pixel of 2.8 mm continuum emission for SDC18.888,SDC23.367, and SDC24.118 are shown in bold, while the thin spectra show the maximum value across all channels in the field.In the case of the C 18 O spectra, a single Gaussian fit has been overlaid, and for the N 2 H + (1-0) spectra, we show the fits obtained with mwydyn, where the parameter designation is given by the formalism in Section 3.3, where p1, p2, p3, and p4 refer to  0 ,  cen , Δ, and  total , respectively.The systemic velocities identified by the single-dish data are shown as dashed vertical lines.

Figure 6 .
Figure6.Top: The fraction of the total mass in cores that resides within the most massive core for the NOEMA clumps as a function of filament convergence parameter.Bottom: the mean density of cores (empty circles) and the mean density of the MMCs, per source, as a function of filament convergence parameter.Clumps with  C > 0.2 are considered to be hubfilament systems.

Figure 8 .
Figure8.Gaussian kernel density estimators illustrating the distribution of the FWHM linewidths of components identified by mwydyn for each source, ordered by 1-pc clump mass.Each data point in the filled distributions is weighted by the integrated intensity of the component, while the distributions shown in outline only were not weighted.The lower distributions for each source represent components within spectra that fall within the mask of the most massive core in each IRDC, while the upper distributions represent components that are not associated with any core.Solid and dashed vertical lines indicate the mean and weighted mean values, respectively.The number of components contributing to each KDE is displayed in the upper left.

Figure 9 .
Figure 9. Summary statistics of dense gas kinematics (rows) as a function parsec-scale properties (columns) for each target.Top row: mean standard deviation of fitted centroid velocities; Second row: mean linewidth; Third row: mean total linewidth.Bottom row: mean number of fitted components per spectrum.These points are shown as functions of parsec-scale properties in each column, including: 1-pc mass, infrared-bright fraction  IRB , filament convergence parameter  C , and virial parameter  vir .In each case the empty circles indicate the values for the total distribution, while the solid circles show the values for the spectra associated with the most massive cores.For σcomp , σline , and σtotal , the error bars show the range covered between the 16th and 84th percentiles of the full distribution, while the error bars for Ncomp give the standard deviation.Square markers, where visible, denote the median values.In all cases, the reported statistics (means, percentiles, and standard deviations) are weighted by the integrated intensity.For each panel, the Pearson correlation coefficients  and -values are given for the MMCs (top) and all spectra (bottom).

Figure 10 .
Figure 10.Pearson correlation coefficients ( P ) and Spearman rank-order correlation coefficients ( S ) for various combinations of quantities investigated in Section 4. The  P values are given to the lower-left of the diagonal, for which stronger red and blue colours indicate stronger positive and negative correlations, respectively, while  S values are given to the upper-right of the diagonal, and stronger green and purple colours indicate stronger positive and negative correlations, respectively.The coefficients for correlations that are statistically significant (i.e. with -values ≤ 0.05) are highlighted with a solid box, while coefficients for correlations that have marginal significance (0.05 <  ≤ 0.10) are highlighted with a dashed box.

Figure A1 .
Figure A1.Case studies demonstrating the need to apply aperture corrections for compact sources extracted with astrodendro.Image: six Gaussian sources with varying amplitudes modelled that are: compact (5 and 6), extended (1, 2, 3, 4), and partially blended (1, 2).The white contour shows the lowest level considered to be considered as part of a source (at a value of 3), while the orange contours show the boundaries of the dendrogram leaves.The bottom three panels show 1-dimensional profiles through the image at the position of the dashed white lines labelled a-c, and the shaded region shows the areas integrated for the integrated flux measurements (see text), respectively.The dashed and dotted black lines show the profile of the individual sources on the left and right of each row, respectively, while the solid blue line shows the combined profile.The numbers to the side of each model profile show the fraction of integrated flux from the model contained within the measurements on the top; the appropriate aperture correction factor is the reciprocal of this number.

Figure A2 .
Figure A2.The aperture correction factor,  ap determined as a function of peak signal-to-noise ratio (SNR) for sources with different areas with respect to the main beam size, denoted by the different colour curves.For extended and high-SNR sources, the correction factor tends to unity.These particular curves were derived to be appropriate for the synthesised NOEMA beam of SDC18.888.

Figure C1 .
Figure C1.Various statistical quantities derived from the mwydyn fitting results to each of the seven IRDCs observed with NOEMA.Each column corresponds to a particular target, and the rows are as follows: a) Integrated intensity of all components; b) Integrated intensity-weighted mean centroid velocity; c) Integrated intensity-weighted centroid dispersion with respect to the systemic velocity; d) Integrated intensity-weighted mean linewidth (i.e.FWHM / √︁ 8 ln(2); e) Total velocity dispersion; f) Number of fitted components per spectrum.
Anderson et al. (2021) new NOEMA observations, while the bottom 6 rows correspond to the adapted ALMA observations ofAnderson et al. (2021).

Table 3 .
Quantities derived for the determination of the virial parameters as described in Section 3.2.

Table 4 .
Core formation efficiencies for the MMCs (  MMC ) and for the sum of all cores (CFE), most-massive core-mass fractions ( MMC ), and the number of cores for each target in our sample.