Modelling complex line emission in the interstellar medium ISM is a degenerate highdimensional problem. Here we present McFine a tool for automated multicomponent fitting of emission lines with complex hyperfine structure in a fully automated way. We use Markov chain Monte Carlo MCMC to efficiently explore the complex parameter space allowing for characterizing model denegeracies. This tool allows for both local thermodynamic equilibrium LTE and radiativetransfer RT models. McFine can fit individual spectra and data cubes and for cubes encourage spatial coherence between neighbouring pixels. It is also built to fit the minimum number of distinct components to avoid overfitting. We have carried out tests on synthetic spectra where in around 90 per cent of cases it fits the correct number of components otherwise slightly fewer components. Typically inlineformulatexmath idTM0001 notationLaTeXTrm extexmathinlineformula is overestimated and inlineformulatexmath idTM0002 notationLaTeXtautexmathinlineformula underestimated but accurate within the estimated uncertainties. The velocity and line widths are recovered with extremely high accuracy however. We verify McFine by applying to a large Atacama Large Millimetersubmillimeter Array ALMA Ninlineformulatexmath idTM0003 notationLaTeX2texmathinlineformulaHinlineformulatexmath idTM0004 notationLaTeXtexmathinlineformula mosaic of an highmass star forming region G316.7500.00. We find a similar quality of fit to our synthetic tests aside from in the active regions forming Ostars where the assumptions of Gaussian line profiles or LTE may break down. To show the general applicability of this code we fit COJ  21 observations of NGC 3627 a nearby starforming galaxy again obtaining excellent fit quality. McFine provides a fully automated way to analyse rich data sets from interferometric observations is open source and pipinstallable.
['Astrophysics - Astrophysics of Galaxies']
McFine PYTHONbased Monte Carlo multicomponent hyperfine structure fitting
{'McFine : python -based Monte-Carlo multi-component hyperfine structure fitting': 'Thomas G. Williams, 1 ★ and Elizabeth J. Watkins 2  \n1 Sub-department of Astrophysics, Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK  \n2 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': 'Modelling complex line emission in the interstellar medium (ISM) is a degenerate, high-dimensional problem. Here, we present McFine , a tool for automated multi-component fitting of emission lines with complex hyperfine structure, in a fully automated way. We use Markov chain Monte Carlo (MCMC) to efficiently explore the complex parameter space, allowing for characterising model denegeracies. This tool allows for both local thermodynamic equilibrium (LTE) and radiative-transfer (RT) models. McFine can fit individual spectra and data cubes, and for cubes encourage spatial coherence between neighbouring pixels. It is also built to fit the minimum number of distinct components, to avoid overfitting. We have carried out tests on synthetic spectra, where in around 90 per cent of cases it fits the correct number of components, otherwise slightly fewer components. Typically, 𝑇 ex is overestimated and 𝜏 underestimated, but accurate within the estimated uncertainties. The velocity and line widths are recovered with extremely high accuracy, however. We verify McFine by applying to a large Atacama Large Millimeter/submillimeter Array (ALMA) N 2 H + mosaic of an high-mass star forming region, G316.75-00.00. We find a similar quality of fit to our synthetic tests, aside from in the active regions forming O-stars, where the assumptions of Gaussian line profiles or LTE may break down. To show the general applicability of this code, we fit CO( J = 2-1) observations of NGC 3627, a nearby star-forming galaxy, again obtaining excellent fit quality. McFine provides a fully automated way to analyse rich datasets from interferometric observations, is open source, and pip -installable.  \nKey words: ISM: abundances - ISM: molecules - galaxies: ISM - ISM: general - methods: data analysis', '1 INTRODUCTION': "Spectral line emission is a key tool for diagnosing the physical conditions of the interstellar medium (ISM, e.g. Walter et al. 2008; Leroy et al. 2009; Rigby et al. 2016; Pety et al. 2017; Liu et al. 2019; Leroy et al. 2021b, amongst many others), but extracting useful and accurate physical parameters from observed line intensities is a complex task. If only a single line is measured (for example, CO in most extragalactic studies), we can access the velocity, gas dispersion, and with a conversion factor 1 the molecular gas mass. With different transitions of the same molecule, we can instead access different phases of the gas (e.g. in density or temperature, see Ohashi et al. 2016; Peñaloza et al. 2018; Puschnig 2020). Usually, tracing different transitions requires multiple telescope configurations covering different frequency bandwidths, but in some instances, molecules can contain hyperfine structure. This hyperfine structure exhibits as multiple lines at slightly offset frequencies, arising from the splitting of the rotational energy levels by electric quadrupole and magnetic dipole interactions, which is caused by nuclear moments of atoms with non-  \n- ★ E-mail: thomas.williams@physics.ox.ac.uk (TGW)\n- 1 This CO-conversion factor, 𝛼 CO , is either assumed or estimated from some other parameter (primarily the metallicity, see Bolatto et al. 2013, and references therein).  \nzero spin. Therefore, we can trace multiple transitions, and thus gas properties such as excitation temperature and optical depth with a single frequency configuration, as modern radio telescopes allow us to resolve these hyperfine lines. These properties offer more insight into the gas conditions than simply measuring the line intensity and velocity since they allow for a detailed study of the local physical conditions of the gas, such as the gas density.  \nThe most straightforward method of calculating the excitation and optical depth is assuming that the kinetic and excitation temperature are equal, known as Local Thermodynamic Equilibrium (LTE), in which case all lines have the same temperature and opacity (e.g. Ohashi et al. 2016). This assumption holds true in highdensity regimes (i.e. significantly above the critical density of the molecules in question), but breaks down in regions of lower density (see Mangum & Shirley 2015). A next step towards a more physically accurate, but complex model, involves knowledge of the collisional properties of the molecules. In this case, a kinetic temperature is modelled, and the excitation temperature and opacity of each line may be different. Local Velocity Gradient (LVG) methods take advantage of this, as well as the popular /r.pc/a.pc/d.pc/e.pc/x.pc tool (van der Tak et al. 2007, which we use for the radiative transfer models in this work). More complex models remove the requirement of local excitation and rather solve the radiative transfer as functions of depth into the cloud (e.g. van Zadelhoff et al. 2002).  \nThe spectra along any given line-of-sight can also show complexity in multiple ways. Many lines will split into multiple hyperfine components, which can be blended together or separated by a small velocity offset. The ISM is continuous and so there may be multiple distinct regions of emission along any given line-of-sight. It is these multiple distinct regions that we will specifically refer to under the term of 'spectral complexity' throughout this work, though we do acknowledge that many individual spectra can also show complex profiles. Spectral complexity within datasets is ubiquitous (e.g. Lindner et al. 2015; Henshaw et al. 2016, 2020; Koch et al. 2021; Rigby et al. 2024), and blindly applying a single component fit to a spectrum with multiple components will at best lead to results that are, on average, correct, and at worst are completely unphysical or misleading (Koch et al. 2021). Defining a number of components per-sightline by hand for the increasing scale of datasets from state-of-the-art radio telescopes is subjective and prohibitively time-consuming, and so this step must be automated along with the fitting of each spectral component.  \nThere are a number of tools available to help automate this task, and for a review of some of these we refer the reader to Riener et al. (2019). Since that work, there also now exist codes more specific for hyperfine structure, such as MUFASA (Chen et al. 2020) and mwydyn (Rigby et al. 2024). However, all of these codes often rely on either capping the maximum number of components that can be fit, or require some human interaction to set the number of components. Given the extremely large volume of datasets now regularly obtained with radio instruments, this is becoming increasingly unwieldy. Additionally, the existing tools are highly susceptible to the degeneracies present in the fitting procedure. A tool that can fit hundreds of thousands of spectra with no manual intervention, and fully explores the complex parameter space, is vital for making progress in characterising the ISM.  \nIn this work, we present McFine , a python -based tool to automatically fit complex spectral line data within a fully Bayesian/MCMC framework. Although built to fit N 2 H + , this tool can be easily extended to other lines and transitions. McFine offers both simple LTE approximation solutions, as well as more complex radiative transfer (RT) solutions based on /r.pc/a.pc/d.pc/e.pc/x.pc calculations. McFine is highly parallelised, and is primarily designed to fit data cubes in an efficient manner, although it can also be used on individual spectra. The code is freely available at https://github.com/ thomaswilliamsastro/mcfine , and is pip -installable 2 . Given the much more complex nature of the /r.pc/a.pc/d.pc/e.pc/x.pc fitting, we will focus on LTEin this paper and defer the RT modelling to a future work. Whilst the RT functionality is present in the current version of McFine , it should be considered relatively untested.  \nThe layout of this paper is as follows. In Section 2, we describe how we build and fit our multi-component spectral model. In Section 3, we apply the fitting routine to a number of synthetic spectra to test the efficacy of the fitting approach. In Section 4, we apply McFine to real data, showing the quality of the fits achieved. Finally, we conclude in Section 5.", '2.1.1 LTE Case': 'In building a hyperfine model, we require an optical depth, 𝜏 𝑖 at the centre of each hyperfine component and an excitation temperature for each component, 𝑇 ex . In the LTE case, 𝑇 ex is the same for all hyperfine components, and the line strengths are fixed so that the sum of each individual 𝜏 𝑖 is a free parameter in the model.  \nWith an excitation temperature and optical depth for each line, we then build a spectrum. We make a number of assumptions at this point, mainly that the opacity is Gaussian in velocity, and that the line width is the same for each hyperfine line 3 . The opacity as a function of velocity is  \n𝜏 𝑣,𝑖 = 𝜏 𝑖 × exp -𝑣 - ( 𝑣 𝑖 + 𝑣 cen ) 2 𝜎 2 , (1)  \nwhere 𝑣 𝑖 is the offset velocity with respect to the velocity of the central line 𝑣 cen . The observed intensity is given by  \n𝑇 = ∑︁ 𝑖 ( 1 -exp [ 𝜏 𝑣,𝑖 ]) × ( 𝑇 𝑖, ex -𝑇 BG ) , (2)  \nwhere 𝑇 BG is the background temperature, typically 2.73 K, the value for the cosmic microwave background. At very low and high values for the opacity ( 𝜏 ≲ 0 . 1, 𝜏 ≳ 30), the 𝜏 term will drop out and the problem becomes degenerate, and (in the low limit, the line ratios do not depend on optical depth and in the high limit the exponential becomes very small), so we add a bound at these points. For each distinct line we therefore have four free parameters in the LTE case -𝑇 ex , 𝜏 , 𝑣 cen and 𝜎 . We fit 𝜏 in natural log-space, which allows the algorithm to more efficiently explore the parameter space 4 . This method is described in more detail in Section 2.2. Given that our line profiles are assumed to be Gaussian, this means we are not able to model, for example, infall, outflow, or self-absorption, which cause deviations from Gaussianity. This is likely to manifest as a poorer fit to the data, which can be explored using the statistics that McFine outputs (see Sect. 3).', '2.1.2 RT Case': 'For the /r.pc/a.pc/d.pc/e.pc/x.pc RT case, 𝑇 ex can vary between hyperfine components and so we generate a grid of models based on the free parameters within /r.pc/a.pc/d.pc/e.pc/x.pc (the kinetic temperature 𝑇 kin , H 2 density 𝑛 H 2 , the column density of the molecule 𝑁 col , and the line width 𝜎 ). From this grid, we extract 𝑇 ex and 𝜏 𝑖 for each component in the hyperfine structure, interpolating between grid values where necessary (a requirement to limit the computational cost of these calculations). To generate the grid, we use /n.pc/d.pc/r.pc/a.pc/d.pc/e.pc/x.pc 5 , a wrapper around /r.pc/a.pc/d.pc/e.pc/x.pc, which we have edited to fix an issue with parsing values from lines with identical frequencies 6 .  \nOverall, we have an extra free parameter bringing the total number to five: 𝑇 kin , 𝑁 col , 𝑛 H 2 , 𝑣 cen and 𝜎 . We enter these values into a pre-generated /r.pc/a.pc/d.pc/e.pc/x.pc grid, which computes a 𝑇 ex and 𝜏 for each  \ncomponent. These, along with 𝑣 cen and 𝜎 , are then used in Eqs. (1) and (2). Here, the column and number densities are both fit in logspace as they have a much larger dynamic range than the other free parameters.', '2.2 Model fitting': "We perform our fitting in a fully Bayesian framework using the MCMC sampler /e.pc/m.pc/c.pc/e.pc/e.pc 7 . Unlike least-squares fitting, MCMC can robustly account for covariance between parameters (particularly important for excitation temperature and opacity). It can also account for asymmetric uncertainties, which we find to be typical in our fits. To fit the data, we build up a steadily more complex model starting from zero components (i.e. a flat line). The (ln-)likelihood function is given by  \nln L = -0 . 5 ∑︁ 𝑖 GLYPH<18> 𝑇 obs -𝑇 model 𝜎 obs GLYPH<19> 2 , (3)  \nwhere 𝑇 obs and 𝑇 model are the observed and modelled intensity, respectively, and 𝜎 obs is the noise in the spectrum. We use 100 walkers (with starting positions perturbed around the initial guess) for each fit. Each walker takes 250 steps of burn-in and then 1000 steps for the main run, and both of these are variables that users can alter as desired. We chose these default parameters from inspecting where the walks converged by eye (i.e. where the walkers in the step plot appear well-mixed; see Figure 1 for an example). We will implement a more rigorous convergence criterion in a future release. We calculate a maximum likelihood from the median parameters of this fit, discarding the first half of the steps as the walkers explore the parameter space. Throughout this fitting, we make use of python 's multiprocessing capabilities either by sending walkers to separate cores in the case of a single spectrum fit, or by passing each spaxel (i.e. the spectrum corresponding to a single 𝑥 -𝑦 position in a cube) fit to a separate core in the case of a data cube.  \nThe MCMC walkers require an initial guess, and for this we find a least-squares fit using lmfit . The initial guesses are configurable we recommend setting an 𝑇 ex similar to the peak temperature in the cube, a ln ( 𝜏 ) of 0, a velocity similar to the systemic velocity of the object and a velocity dispersion of 1 km s -1 , for most cases. However, as long as these parameters are reasonable then lmfit should be able to find an acceptable guess to initialise the MCMC. By default, we use the global basinhopping method (although this can be tweaked as users require) as we find it to produce good initial estimates even into a high-dimensional parameter space. Note that the exact details of the least-squares fit are not vital, but should simply provide a good initial guess for the parameters. If the initial guesses are far from the 'true' parameters, the MCMC will typically not converge and the fit will be poor. This is especially true for more complex fits as the parameter space becomes increasingly more complicated. This is a general problem in fitting data, and not unique to our problem.  \nAfter fitting for one component, we gradually increase the complexity by adding components one at a time. When we increase the number of components, we do not use any of the parameters from the fewer component fits to guide our procedure. Each fit is done completely independently, and then compared at the end. We discriminate between the more complex and previous model using the the Bayesian information criterion (BIC). This is given by  \nBIC = 𝑘 ln ( 𝑛 ) -2 ln (L) , (4)  \nTable 1. Bounds for randomly drawn parameters in our synthetic tests. 𝑛 comp is number of components, 𝑇 ex the excitation temperature, 𝜏 the optical depth, 𝑣 𝑖 the offset from the -40 km s -1 centre of the spectrum, and 𝜎 the line width.  \nwhere 𝑘 is the number of fitted parameters, and 𝑛 the number of points that are fit to. A change in BIC, Δ BIC > 10 is typically used as indicating a strong preference for one model over another (see, e.g. Koch et al. 2021). We continue to add components until the change in BIC is less than 10, at which point we stop adding components and take the less complex fit as our end point.  \nAlthough this procedure does well at identifying the majority of the spectral features, the observational noise is sometimes modelled as a faint, wide spectral feature. We therefore take an additional 'top-down' approach where we remove components from faintest to strongest (in terms of integrated intensity). We then fit using the median value of each fitted parameter from the MCMC run as a starting point, and calculate the Δ BIC again. This assumes that the faintest component is most likely to be noise, not necessarily that it is the 'worst' fit - in some cases, these components are indeed real, and the comparison will retain such features. If the model with fewer components is not strongly preferred (i.e. Δ BIC < 10), then we stick with our original, more complex model.  \nFinally, for spectral cubes, neighbouring spaxels are not independent (both because the pixel size is selected to oversample the beam, and the gas properties should be at least somewhat correlated over short distances). We therefore have an option to compare the fit parameters to the neighbouring eight pixels. If any of these are strongly preferred by the BIC, we take the most preferred model. We loop over pixels in RA and Dec, first forwards and then backwards. This ends up replacing less than one per cent of pixels (see Sect. 4) in our tests on real data sets. Note that at this step we do not perform any refitting, but rather compare the spectrum from the median of each fitted parameter to the observed spectrum in each spaxel.", '3 TESTS WITH SYNTHETIC SPECTRA': "Totest the efficacy of our modelling approach, we generate a selection of models and apply McFine . We will focus on the LTE approach in this work for simplicity. We generate a number of synthetic spectra with a velocity resolution of 0 . 2 km s -1 , Gaussian noise with a standard deviation of 0.1 K and a velocity range between -80 and 0 km s -1 . This spectral resolution, range, and noise level are all approximately matched to the real ALMA data we will test later . However, we increase the velocity range to account for the potential wide separation in components for these tests, to more appropriately match those seen between large-scale dynamic features in galaxies (e.g. Beuther et al. 2017). We also include a 5 per cent 'calibration' uncertainty, to avoid biasing high 𝜒 2 𝜈 values for bright spectra (see Sect. 4). The number of components and parameter bounds for each component are given in Table 1. We generate and fit a total of 100 models to explore this parameter space and the limitations of the fitting. We show some of these fits in Appendix A.  \nFigure 1. Step plot for one of our synthetic spectra tests. We can see the individual walkers are well mixed, and have converged to a stable value for each parameter.  \n<!-- image -->  \nFigure 2. Distribution of reduced chi-square values for our synthetic spectra tests. The black point and error-bar show the median, and the 16 th and 84 th percentiles of the distribution (also given in the text box). 𝜒 2 𝜈 = 1 is shown as a dashed vertical line.  \n<!-- image -->  \nFigure 3. Distribution of input minus fitted components for our 100 generated spectra. We colour the histogram by the number of input components, and indicate where the number of fitted components is equal to the number of input components as a dash-black, vertical line.  \n<!-- image -->  \nv  \nFigure 4. Comparison between input and fitted parameters for the models where the correct number of parameters are fitted. The points are coloured by parameter to match those of Figure 5. In each case, the abscissa shows the input parameter value, and the ordinate the fitted parameter. In each subplot we also show the median and offset values from the truth in a text box in the top left, and the 1:1 line is shown as a dashed black line. The small inset plot shows the residual, having removed this 1:1 line, to make any deviations clearer.  \n<!-- image -->  \nOur first metric in this comparison is the reduced chi-square, 𝜒 2 𝜈 , which we calculate from the 𝜒 2 value divided by the number of free parameters ( 𝑛 -𝑘 ). The distribution of this is shown in Fig. 2, and it shows that the models perform well with a value 𝜒 2 𝜈 = 1 . 01 + 0 . 09 -0 . 07 , consistent with the expected 𝜒 2 𝜈 = 1 for a perfectly fit model. The 𝜒 2 𝜈 distribution is approximately Gaussian, indicating no particular bias towards under- or overfitting.  \nSecondly, we investigate the number of components recovered in our synthetic tests, and we summarise this in Fig 3. Around twothirds of the models recover the correct number of components, with no fits including more components than input. This indicates that McFine successfully avoids overfitting spectra. Some of the fits have fewer components than input, and when we examined these found them to be models with highly blended components. This causes the model to preferentially fit a single, brighter component than attempting to split into multiple components. Reassuringly, we find no broad components fitted to the noise. In general, the number of  \ncomponents fit by McFine should be considered a lower limit to the true spectral complexity.  \nFor the 89 models that fit the correct number of parameters, we compare the fitted parameters for each component to the input model parameters, and we show these in in Figure 4. We see a good correspondence between the input and fitted parameters, particularly with the velocity and velocity width. Typically, 𝑇 ex tends to scatter above the one-to-one line and 𝜏 scatters low, and this is driven by the degeneracy between these two parameters. We also show the distributions of these parameters normalised by their associated errors in Figure 5. We see that despite some systematic underestimating of parameters, the fitted parameters are within one standard deviation of the true value in ≃ 50 per cent of cases (we should expect this number to be 68 per cent, which is similar). The exception to this is 𝜏 , which is very slightly outside the 1𝜎 range.  \nThese synthetic tests show that McFine typically produces very good fits to input data, and so we are confident in our fitting procedure going forward and apply it to real data.  \n<!-- image -->  \nX  \nfit  \nFigure 5. Difference between fitted and input parameters over our synthetic models, normalised by the associated errors. Lines are coloured by the particular parameter, which are given in the Figure legend. The vertical, dashed black line indicates where the fitted and input parameter are equal.Figure 6. Moment 0 (integrated intensity) map of G316.75. A physical scalebar is shown in the bottom right corner.  \n<!-- image -->", '4 APPLICATION TO ALMA DATA': 'In this Section, we provide two benchmarks of the code on ALMA data. The first is G316.75, a region in the Milky Way (MW) observed in the N 2 H + line, and which we have focussed McFine on fitting. Secondly, to demonstrate the versatility of the fitting code, we apply McFine to CO( J =2-1) observations of NGC 3627 (a nearby starforming galaxy).', '4.1 G316.75': "For the first real-world benchmark for McFine , we apply the fitting to an N 2 H + ( J =1-0) mosaic of G316.75-00.00 (G316.75 hereafter), which has been fully observed with ALMA using the 12m, 7m, and total power (TP) arrays at ∼ 4 '' resolution, equivalent to a physical resolution of 0.036 pc (project codes 2013.1.01029.S and 2018.1.01851.S). These observations comprise three large 12m, 7m, and total power (TP) observations to be sensitive to emission at all an-  \nscales. The integrated intensity map for this mosaic is shown in Figure 6. G316.75-00 is a dense, high-mass, star-forming molecular filament located 2.69 ± 0.45 kpc away in the Galactic plane (Watkins et al. 2019). The filament contains ∼ 20,000 M ⊙ of molecular gas over a 14 pc length and exhibits a range of star-forming features in different evolutionary stages. More specifically, the southern half of the filament is dominated by an active H/i.pc/i.pc region powered by 2-4 O-stars. The northern half - which is just as massive - is a quiescent infra-red dark cloud (IRDC) containing large star-forming clumps and has hints of a more intricate sub-filament network fuelling their growth in dust emission.  \nSince Watkins et al. (2019) found that G316.75 has overlapping velocity components using lower resolution (0.5 pc) N 2 H + observations (Foster et al. 2011, 2013; Jackson et al. 2013), it is an excellent test case for McFine . Due to the range of evolutionary states present, we expect complex, non-Gaussian velocity motions throughout the filament, including outflows and inflows, as well as conditions that might challenge the thermalised assumption of LTE modelling. These features will allow us to evaluate McFine 's ability to extract distinct velocity structure in extreme, non-ideal environments, and gauge the underlying limitations of our models.  \nThe full details of the data reduction will be given in Watkins et al. (in preparation), but in brief, we image each mosaic separately, as they have somewhat different noise parameters and observing conditions. Broadly following the strategy of (Leroy et al. 2021a), we combine the 12m and 7m data into a single measurement set (MS) per mosaic, but instead of using tclean on the interferometric data and then feathering the TP after deconvolution, we use the experimental sdintimaging routine to include the TP imaging in the deconvolution. The full details of this algorithm are given in Rau et al. (2019), and we find this correctly recovers the total flux compared to the TP map, whereas using the more typical tclean plus feather method underestimates the flux by up to 30 per cent in the brightest regions. We note that this is typically not an issue in extragalactic observations (as deviations from tclean/feather and sdintimaging only start to arise at a beam temperature perchannel of around 2 K), but will likely be more problematic in Milky Way targets and extremely bright regions of extragalactic objects. We first perform a shallow multi-scale clean down to 4 𝜎 , followed by a deep single-scale clean to 1 𝜎 , or a relative flux change in the model flux of less than 0.1 per cent. Finally, we linearly combine the mosaics based on their noise properties, and homogenise the beam to a common, round beam.  \nTo define a mask that will contain almost all emission, we take the strict mask produced by the reduction pipeline. The details of the mask are given in section 7.3 of Leroy et al. (2021a), but briefly these grow 4𝜎 peaks of emission down to a 2𝜎 level, thus capturing most emission but maintaining a low false-positive rate. We take any spaxel that has significant emission anywhere down the line of sight as defined by the mask cube. This leads to a total of ∼ 100,000 pixels. Werun McFine onthis cube on a compute server, using 96 CPU cores. This server achieves around 1,000 fits per hour, although individual fit time is highly dependent on the model complexity (i.e. a two component model will take around a factor of two times longer to fit, as it must first fit a single component model). We then loop forth and back over RA and Dec, comparing the models of neighbouring pixels to encourage spatial coherence. We find this replaces ∼ 700 fits, less than one per cent. This would indicate that in the vast majority of cases we successfully find the best fit in our first pass. In total, this fitting takes around 120 hours. This means that each fit takes around 5 mins on average on a single CPU core, and so for large data  \nFigure 7. Position-position-velocity scatter plot of the McFine fits to the G316.75 ALMA data. Points are coloured by the measured velocity dispersion and the overall transparency is intensity-weighted. The projections along each axis are shown in greyscale, so that in the position-position plane, this is equivalent to the number of components along the line-of-sight.  \n<!-- image -->  \nvolumes we highly recommend using a high-performance compute cluster.  \nA position-position-velocity (PPV) rendering of the fits is shown in Figure 7. A rich variety of structure is seen that would be missed in traditional moment analysis; there is complex structure along the ridge of the IRDC, but also a significant amount of emission above and below the filament. This off-filament gas often has a higher velocity dispersion, indicative of infall or gas ejection (e.g. Klessen & Hennebelle 2010; Watkins et al. 2019). It is also clear that a full 3-dimensional treatment is necessary to build up associations of coherent velocity structures - any single plane of this PPV space misses crucial information.  \nWe will defer a full exploration of these fits to Watkins et al. (in preparation), but as a first demonstration of the fit quality we show the fitted 𝜒 2 𝜈 map and distribution in Fig. 8. Note that in calculating 𝜒 2 𝜈 , we include a 5 per cent systematic error added in quadrature to the RMS uncertainty, similar to the scatter in flux calibration seen in Leroy et al. (2021a). We include this to account for the fact that even if the fit quality looks good in the brightest regions, the very small RMS uncertainties lead to extremely large values of 𝜒 2 𝜈 . In general, we produce extremely good fits to the data, but there are some exceptions (see Appendix B). Particularly around the south-  \negion with high 𝜒 2 𝜈 values. This region contains young O-stars and the gas flows around these stars likely violate the LTE assumptions or Gaussian line profiles. We see in the fits that the fitted 𝜏 values are either very low or very high (reflecting the degeneracy highlighted in Eq. 2), and that often the fitted excitation temperatures will also be either very high or very low. We therefore suggest inspection of the corner plots (see Fig. B2 for an example) for fits where 𝜏 is at the prior bounds, as these are likely to be the regions where the LTE assumptions are inappropriate. These relatively poor fits, compared to the rest of the data, therefore reflect limitations in the simple profile for the line model, rather than a breakdown in the fitting routine.  \nSecondly, we investigate the number of components in Fig. 9. The majority (55 per cent) of our spectra are well-described by a single component, with a maximum of seven distinct components (0.02 per cent of the fits). We find that the spectral complexity increases towards the central filament, as well as towards the star forming cores and clumps residing there. This can cause non-Gaussian line profiles (e.g. Peretto et al. 2020; Rigby et al. 2024). We leave the full interpretation to Watkins et al. (in preparation). Just under half of the fits required more than one component, indicating that this spectral complexity is ubiquitous throughout the filament, further  \nFigure 8. Left: 𝜒 2 𝜈 map for each fit within the mask of G316.75. Right : Distribution of 𝜒 2 𝜈 values. The black point and errorbar shows the median value and 16 th /84 th percentiles (also given in the box in the top right). A vertical dashed black line shows 𝜒 2 𝜈 = 1.  \n<!-- image -->  \nFigure 9. Left: 𝑛 comp map for each fit within the mask of G316.75. Right : Histogram of 𝑛 comp values.  \n<!-- image -->  \nRA (J2000)  \nmotivating the need for automated, multi-component fitting. Despite there being no prior constraint of the number of components to fit, we find that the code typically prefers to fit a few components, rather than many weaker components.  \nTo test any potential dependence on the quality of fits with the number of fitted components, we produce violin plots of 𝜒 2 𝜈 against 𝑛 comp in Fig. 10. There is a trend towards higher 𝜒 2 𝜈 with increasing 𝑛 comp . Reassuringly, this trend is opposite to the idea that increasing the number of components may lead to overfitting (in which case, 𝜒 2 𝜈 would systematically decrease with increasing 𝑛 comp ). In the vast majority of cases then, McFine produces reasonable fits to the data. We see some evidence of underfitting (i.e. potentially missed components in the spectrum), which we diagnose by manual inspection of the fits. These are either very narrow and weaker components, or very broad, very faint components - in this case, adding an extra  \ncomponent does not substantially change the Δ BIC, and so McFine does not include this component. We prefer missing some components for the benefit of ensuring the fit components are trustworthy, although users can fine-tune this using the Δ BIC if they desire.", '4.2 NGC 3627': "Our second real-world benchmark is a nearby star-forming galaxy observed as part of the Physics at High Angular Resolution (PHANGS 8 ) survey (Leroy et al. 2021c), NGC 3627 (project code ADS/JAO.ALMA#2017.1.00886.L). The details of the data processing are given in Leroy et al. (2021a), and are generally the same  \n10  \nN  \nFigure 10. Violin plots showing 𝜒 2 𝜈 distribution with fitted number of components, 𝑛 comp . Each violin is coloured by the number of points fitted with that number of components, and in each violin we show the median value, along with the 16 th and 84 th percentiles of the distribution.  \n<!-- image -->  \nFigure 11. Moment 0 (integrated intensity) map of NGC 3627. A physical scalebar is shown in the upper right corner.  \n<!-- image -->  \nas G316.75, with the exception of the line, which is CO( J =2-1) instead of N 2 H + , the cube is imaged at 2.5 km s -1 resolution, and feather is used instead of sdintimaging . NGC 3627 is a strongly barred galaxy with complex dynamics, particularly at the bar ends (Beuther et al. 2017), making this a good test for McFine . Using the PHANGS-ALMA data from the first public data release, we fit the cube in the same way as in Section 4.1, although for the CO line there is no hyperfine splitting, simplifying the spectra considerably. Given that we only have a single component per spectrum fitted, this makes 𝑇 ex and 𝜏 completely degenerate with each other, although the velocity and dispersion remain useful and well-constrained. Again, we fit around 100,000 pixels, and encouraging spatial coherence replaces around 1 per cent of fits. As the CO line is only modelled as a single Gaussian, and we generally fit fewer components, the fitting here is faster than in Section 4.1, taking around 40 hours.  \nWe show the 𝜒 2 𝜈 and 𝑛 comp maps in Figures 12 and 13. Unlike for G316.75, we see an extremely uniform distribution of 𝜒 2 𝜈 values across the entire galactic disc, with a peak very close to one, and a very narrow spread of around 0.1 (16 th and 84 th percentiles). This shows the complexity of our model is sufficient to describe the line profiles in all ∼ 100 pc regions of these CO observations. The peak of the 𝜒 2 𝜈 distribution is in the bar ends of this galaxy, which is the most dynamically active region with the highest level of spectral complexity. We also see that the majority of the CO lines-of-sight can be fit with a single component, but the spectral complexity increases towards the galaxy centre and the bar ends. Even in these relatively low-resolution CO observations, in dynamically complex regions, a tool like McFine is required to properly model the gas dynamics.  \nTo highlight the importance of accounting for multiple velocity components along each line-of-sight, we show a comparison between the available PHANGS moment 1 (mom1; intensity-weighted average velocity) and effective width (a proxy for the velocity dispersion, see Leroy et al. 2021a, for a discussion of this quantity) to our calculated values in Figure 14. For the PHANGS maps, we use the 'strict' masked maps, which trades flux completeness for a lower false-positive rate. In the McFine case, when we have multiple components along the line-of-sight, we perform a peak temperatureweighted average of all the values, having removed any parameters that are not detected above a 5 𝜎 limit, to reduce the noise. We also only show sightlines where we fit multiple components, as this is where we expect to see the majority of the deviations between the two methods. We can see that the velocity values are relatively similar, as expected as the mom1 map is also intensity averaged. However, the effective widths are significantly biased along sightlines with higher spectral complexity, with the effective width being significantly higher than the average velocity dispersion McFine calculates. This is because the effective width blindly takes all emission along a sightline into account, leading to extremely large effective widths for regions with multiple components more widely spaced in velocity. In dynamically complex regions, it is therefore crucial to account for this spectral complexity. Henshaw et al. (in preparation) perform a similar comparison across the full PHANGS sample, finding very comparable results.  \nAs one further test, we compare McFine to another multicomponent fitting code. In Figure 15, we show a comparison of the number of components fitting by McFine to those found by scouse , a semi-automated multi-component line fitting technique (Henshaw et al. 2016, 2019), which are presented in Henshaw et al. (in preparation). In the majority of cases, McFine and scouse fit the same number of components. In around 30 per cent of cases, scouse fits more components than McFine , but the opposite is much less fre-  \nFigure 12. As Figure 8, but for NGC 3627. To highlight the much more uniform 𝜒 2 𝜈 , we use the same colourbar scale as in Figure 8.  \n<!-- image -->  \nquent (3 per cent). In general, these two codes suggest very similar levels of spectral complexity of NGC 3627.", '5 CONCLUSIONS': "In this work, we have presented McFine , an MCMC tool for fitting complex spectral profiles in astronomical data. McFine works by building up steadily more complex models, and performing comparisons using the Bayesian Information Criterion to decide an acceptable stopping point for the complexity of the model. In the case of data cubes (e.g. from ALMA), it can also encourage spatial coherence by comparing neighbouring models. In a number of synthetic tests, we have shown that McFine overall is successful at recovering our input parameters, fitting the correct number of components in around 90 per cent of cases. We tend to recover the correct parameters (within the fitting uncertainties) for spectra where we fit the correct number of components. For fits with too few components, either the lines are highly blended and confused with a single brighter line, or the line is faint enough that it did not strongly affect the combined spectrum. Our fits well-describe the data, with a median 𝜒 2 𝜈 value of 1.01 across our synthetic tests, and a small spread of around 0.1 (16 th and 84 th percentiles).  \nWe also tested McFine on two large ALMA mosaics as a realworld example. For the first, we use LTE modelling to fit ≃ 100,000 N 2 H + spectra from a Galactic molecular filament. McFine most  \ncommonly fits a single component, with a maximum of seven components (0.02 per cent of the fits). There are a significant number of spaxels (45 per cent) that require more than one component. The 𝜒 2 𝜈 distribution is very close to a value of one, indicating that on average the fits almost perfectly describe the data. However, there are some regions that have elevated 𝜒 2 𝜈 values. These are complex, active regions of star formation where the LTE or Gaussian line assumptions are likely to break down due to the complex geometry and motions of the gas. These regions have either very high or low 𝜏 values, and are likely where full radiative transfer modelling is required. We find a slight increase in the 𝜒 2 𝜈 values with the number of fitted components, indicating that this spectral complexity is needed to help properly model the gas in this region, rather than McFine overfitting the data.  \nAlthough McFine has primarily been developed to work with N 2 H + emission, it is relatively simple to extend to other lines via simple dictionaries within the code. In our second test, we demonstrate this adaptability by fitting an ALMA CO( J =2-1) mosaic of NGC 3627, a nearby star-forming galaxy. We find a smooth 𝜒 2 𝜈 distribution across all of our fitted spectra, showing that the model successfully captures the variations in the data. The spectral complexity increases towards the ends of the bar and the galaxy centre, as we would expect in dynamically complex regions. We have shown that accounting for these multiple components is critical to obtain reasonable measures of the velocity dispersion of the gas. In comparison to  \nFigure 13. As Figure 9, but for NGC 3627.  \n<!-- image -->  \nanother multi-component fitting code, scouse , we generally fit the same number of components, with scouse fitting more components in around 30 per cent of spaxels, but with an overall similar 𝜒 2 𝜈 value. McFine fits more components than scouse in around 5 per cent of cases. This comparison gives us confidence that our methodology does not lead to wildly different results to other fitting philosophies.  \nAs radio telescopes such as ALMA and NOEMA push towards large mosaics of thousands of independent points covering multiple diagnostic emission lines, tools like McFine are required to interpret these rich datasets, rather than a more bespoke, 'hands-on' approach. As such, an open-source tool like McFine will be invaluable for the community to use in their science.", 'ACKNOWLEDGEMENTS': 'We thank the reviewer for their comments which have improved the quality and presentation of this work. The authors would also like to thank Michael Anderson, Andrew Rigby, and Nicolas Peretto for useful discussions, and to Jonathan Henshaw for sharing the scouse fits for NGC 3627. EJW acknowledges the support of STFC consolidated grant number ST/X001229/1.  \nThis paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.01029.S, ADS/JAO.ALMA#2018.1.01851.S, ADS/JAO.ALMA#2017.1.00886.L. ALMA is a partnership of  \nESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.', 'DATA AVAILABILITY': 'The source code for McFine is available at https:// github.com/thomaswilliamsastro/mcfine , and the testing scripts at https://github.com/thomaswilliamsastro/ mcfine\\_testing . McFine can be installed with pip install mcfine , although the RT capabilities will require gfortran to be in the user PATH. The PHANGS-ALMA data is publicly available online at https://www.canfar.net/storage/list/ phangs/RELEASES/PHANGS-ALMA/ . The G316.75 data will be presented in Watkins et al. (in preparation).', 'REFERENCES': 'Beuther H., Meidt S., Schinnerer E., Paladino R., Leroy A., 2017, A&A, 597, A85  \nBolatto A. D., Wolfire M., Leroy A. K., 2013, ARA&A, 51, 207 Chen M. C.-Y., et al., 2020, ApJ, 891, 84 Foster J. B., et al., 2011, ApJS, 197, 25  \n<!-- image -->  \nn  \nn  \nFigure 14. Left : Comparison between mom1 velocities from PHANGS and the peak temperature-weighted average of McFine fitted velocity components along each line-of sight. The one-to-one correspondence is shown as a black dashed line. Right : Comparison between equivalent width from PHANGS and peak temperature-weighted average of McFine fitted velocity dispersions for each component along each line-of-sight. Again, the one-to-one correspondence is shown as a black dashed line. Points are coloured by the number of fitted McFine components.Figure 15. Comparison between number of components fitted by McFine and scouse . The vertical dashed black lines indicates where the same number of components are fit by both codes. Increasingly positive (negative) values indicate relatively more components fitted by McFine ( scouse ).  \n<!-- image -->  \nFoster J. B., et al., 2013, Publ. Astron. Soc. Australia, 30, e038 Henshaw J. D., et al., 2016, MNRAS, 457, 2675 Henshaw J. D., et al., 2019, MNRAS, 485, 2457 Henshaw J. D., et al., 2020, Nature Astronomy, 4, 1064 Jackson J. M., et al., 2013, Publ. Astron. Soc. Australia, 30, e057 Klessen R. S., Hennebelle P., 2010, A&A, 520, A17 Koch E. W., et al., 2021, MNRAS, 504, 1801  \nLeroy A. K., et al., 2009, AJ, 137, 4670  \nLeroy A. K., et al., 2021a, ApJS, 255, 19 Leroy A. K., et al., 2021b, ApJS, 257, 43 Leroy A. K., et al., 2021c, ApJS, 257, 43 Lindner R. R., et al., 2015, AJ, 149, 138 Liu D., et al., 2019, ApJ, 887, 235 Mangum J. G., Shirley Y. L., 2015, PASP, 127, 266 Ohashi S., Tatematsu K., Fujii K., Sanhueza P., Nguyen Luong Q., Choi M., Hirota T., Mizuno N., 2016, PASJ, 68, 3 Peñaloza C. H., Clark P. C., Glover S. C. O., Klessen R. S., 2018, MNRAS, 475, 1508 Peretto N., et al., 2020, MNRAS, 496, 3482 Pety J., et al., 2017, A&A, 599, A98 Puschnig J., 2020, Dense Gas Toolbox - Calculate Density and Temperature from Observed Molecular Lines, doi:10.5281/zenodo.3686329 Rau U., Naik N., Braun T., 2019, AJ, 158, 3 Riener M., Kainulainen J., Henshaw J. D., Orkisz J. H., Murray C. E., Beuther H., 2019, A&A, 628, A78 Rigby A. J., et al., 2016, MNRAS, 456, 2885 Rigby A. J., et al., 2024, MNRAS, 528, 1172 Walter F., Brinks E., de Blok W. J. G., Bigiel F., Kennicutt Robert C. J., Thornley M. D., Leroy A., 2008, AJ, 136, 2563 Watkins E. J., Peretto N., Marsh K., Fuller G. A., 2019, Astronomy & Astro- physics, 628, A21 van Zadelhoff G. J., et al., 2002, A&A, 395, 373 van der Tak F. F. S., Black J. H., Schöier F. L., Jansen D. J., van Dishoeck E. F., 2007, A&A, 468, 627', 'APPENDIX A: EXAMPLE SYNTHETIC FITS': 'Here, we give the fits for the lowest, median, and highest 𝜒 2 𝜈 fits from our synthetic tests in Sect. 3. In the case for the model with the worst 𝜒 2 𝜈 , we see that McFine fits fewer components than the input spectrum, leading to a worse fit to the overall spectrum.  \nFigure A1. Fit spectrum for the synthetic model with the lowest 𝜒 2 𝜈 in our tests. The red line shows the fit to the data, shaded red region the fit uncertainty. In the upper-right corner, we give the 𝜒 2 𝜈 value, and the input and fitted number of components.  \n<!-- image -->  \nFigure A2. As Fig. A1, but for the fit with the median 𝜒 2 𝜈 .  \n<!-- image -->', 'APPENDIX B: EXAMPLE ALMA FITS': 'Here, we give the fits for the lowest, median, and highest 𝜒 2 𝜈 fits from our real-world tests to G316.75 in Sect. 4. For the fit with the median 𝜒 2 𝜈 , we also show the corner plot.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \nFigure A3. As Fig. A1, but for the fit with the worst 𝜒 2 𝜈 . The faint black lines show the individual components, and the faint red lines the fit to them. In this case, we fit one fewer component than input, but this is due to there being two relatively faint, blended components.  \n<!-- image -->  \nFigure B1. Left: Fit spectrum for the model with the lowest 𝜒 2 𝜈 in our ALMA data. The red line shows the fit to the data, shaded red region the fit uncertainty. We also show the fit and uncertainty for each individual component. In the upper-right corner, we give the 𝜒 2 𝜈 value. Right : Integrated intensity map for the region, with the location of this spaxel shown as a bright green cross.  \n<!-- image -->  \nFigure B2. As Fig. B1, but for the fit with the median 𝜒 2 𝜈 .  \n<!-- image -->  \nFigure B3. Corner plot for the fit shown in Figure B3.  \n<!-- image -->  \nFigure B4. As Fig. B1, but for the fit with the highest 𝜒 2 𝜈 .  \n<!-- image -->'}