NEATH III: a molecular line survey of a simulated star-forming cloud

We present synthetic line observations of a simulated molecular cloud, utilising a self-consistent treatment of the dynamics and time-dependent chemical evolution. We investigate line emission from the three most common CO isotopologues ($^{12}$CO, $^{13}$CO, C$^{18}$O) and six supposed tracers of dense gas (NH$_3$, HCN, N$_2$H$^+$, HCO$^+$, CS, HNC). Our simulation produces a range of line intensities consistent with that observed in real molecular clouds. The HCN-to-CO intensity ratio is relatively invariant with column density, making HCN (and chemically-similar species such as CS) a poor tracer of high-density material in the cloud. The ratio of N$_2$H$^+$ to HCN or CO, on the other hand, is highly selective of regions with densities above $10^{22} \, {\rm cm^{-2}}$, and the N$_2$H$^+$ line is a very good tracer of the dynamics of high volume density ($>10^4 \, {\rm cm^{-3}}$) material. Focusing on cores formed within the simulated cloud, we find good agreement with the line intensities of an observational sample of prestellar cores, including reproducing observed CS line intensities with an undepleted elemental abundance of sulphur. However, agreement between cores formed in the simulation, and models of isolated cores which have otherwise-comparable properties, is poor. The formation from and interaction with the large-scale environment has a significant impact on the line emission properties of the cores, making isolated models unsuitable for interpreting observational data.


INTRODUCTION
Star formation occurs in molecular clouds, where high gas densities and strong shielding from the external ultraviolet (UV) radiation field result in most hydrogen existing in the form of H2 (Bergin & Tafalla 2007).Temperatures in molecular clouds are too low to excite any significant line emission from H2, making the bulk of the star-forming gas essentially invisible.Studies of star formation therefore rely on indirect tracers of the gas mass, such as far-infrared (far-IR) thermal emission from dust grains (e.g.Könyves et al. 2015) or extinction measurements of background stars (Zucker et al. 2021).Of particular utility are rotational emission lines from molecular species at millimetre wavelengths.Line emission provides information on the line-of-sight velocity of the gas, crucial for assessing its dynamics and stability, and the large number of lines detectable with modern observational facilities (e.g.Pety et al. 2017;Kauffmann et al. 2017;Tafalla et al. 2021) can be used to estimate properties such as volume density and temperature, exploiting the fact that the excitation ⋆ Email: priestleyf@cardiff.ac.uk of each line has a unique response to the local physical conditions (Shirley 2015).
Analysis of molecular line data is complicated by the fact that the observed emission depends on a complex combination of the physical structure of the cloud, radiative transfer effects, and the chemical composition of the gas.Disentangling all of these factors to accurately reconstruct detailed cloud properties from observational data is not practically feasible.Many studies have therefore investigated forward modelling the problem and producing synthetic line observations of simulated clouds, which can be used to help inform a physical interpretation of real observational data.To date, most of these studies have focused on lines from important coolant species, in particular CO (e.g.Peñaloza et al. 2017;Seifried et al. 2017;Peñaloza et al. 2018;Clarke et al. 2018;Clark et al. 2019), as its chemical evolution is often followed self-consistently in modern hydrodynamical simulations.However, due to its high abundance and correspondingly high line optical depths, CO is a poor tracer of the high-density material where star formation actually occurs in molecular clouds (Clark et al. 2019;Priestley et al. 2023a).Isotopologues of CO, such as C 18 O, have much lower optical depths and hence are better tracers of dense gas, but they are also faint and hard to observe in distant clouds or extragalactic systems.
Studies investigating lines from rarer molecules such as HCN and N2H + , which can probe denser gas more effectively, generally assume that the molecular abundances are constant, or can be specified as a function of the local physical properties (e.g.Offner et al. 2008;Smith et al. 2012Smith et al. , 2013;;Chira et al. 2014;Jones et al. 2023;Jensen et al. 2023).This approach ignores the drastic variations in abundance driven by differing evolutionary histories, which become particularly important at these high densities (Priestley et al. 2023c).Obtaining reliable line emission properties from simulated molecular clouds requires the physical and chemical evolution of the clouds to be treated self-consistently.
In this paper, we achieve this by performing radiative transfer simulations of a cloud for which the molecular composition has been obtained from a full time-dependent chemical network.The abundances at each point in the cloud are thus consistent with the physical evolution of the material that ended up at that point, making the resulting line emission data genuinely comparable to that from real molecular clouds with their own complex formation histories.We use this to assess the degree of correspondance between our simulated cloud and reality, the limitations of smaller-scale models in interpreting line emission, and the accuracy of commonlyused observational diagnostics of the physical structure of clouds.

Hydrodynamical simulation
We simulate the formation and subsequent evolution of a molecular cloud via the collision of two spherical, initiallyatomic gas clouds.The simulations are performed using the magnetohydrodynamic (MHD) moving-mesh code arepo (Springel 2010;Pakmor et al. 2011), modified to properly capture the thermodynamics of the gas and dust (Glover & Mac Low 2007;Glover & Clark 2012).This includes a simplified chemical network for H2 and CO (based on Gong et al. 2017, with some additions and modifications described in Hunter et al. 2023), and a self-consistent treatment of shielding from the background UV radiation field (Clark et al. 2012a).
The initial conditions of the simulation are two spherical 10 4 M⊙ clouds with radii R = 19 pc, giving an initial volume density of hydrogen nuclei nH = 10 cm −3 .The gas and dust temperatures are initialised to 300 K and 15 K respectively.The clouds are displaced in the ±x direction by R, so they are just in contact at the outset, and given a bulk motion along the x-axis of ∓7 km s −1 .Each cloud also has a virialised turbulent velocity field with 3D velocity dispersion σ = 0.95 km s −1 , and a 3 µG magnetic field is present parallel to the collision axis.In collapsing high-density regions of the simulation, sink particles are introduced (following Tress et al. 2020; see also Prole et al. 2022), with a threshold density of 2×10 −16 g cm −3 and a formation radius of 9×10 −4 pc.The simulation is run for 5.53 Myr, at which point sink particles have accreted 102 M⊙ of material (0.5% of the original cloud masses); beyond this point, feedback effects from newlyformed stars (which are not modelled) are likely to become important (Whitworth 1979;Walch et al. 2012).
As the mesh cells can be created and destroyed, and so do not correspond to coherent parcels of gas, we follow the physical evolution of representative gas parcels using Monte Carlo tracer particles (Genel et al. 2013), which provide the input for the chemical modelling described below.Global properties, shared between the MHD and chemical models, are the UV radiation field strength of 1.7 times the Habing (1968) field, the cosmic ray ionisation rate 1 of 10 −16 s −1 per H atom, the dust-to-gas mass ratio of 0.01, and the elemental abundances of carbon (1.4 × 10 −4 ) and oxygen (3.2 × 10 −4 ) from Sembach et al. (2000).The 'metal' abundance in the MHD simulation, representing heavier elements such as silicon, is set to 10 −7 , corresponding to the high levels of depletion found in the dense interstellar medium (Jenkins 2009).

Chemical modelling
The tracer particles in the MHD simulation record the density, gas temperature, and effective shielding column densities at intervals of 44 kyr.These particle histories are used to model the chemical evolution under the NEATH framework 2 (Priestley et al. 2023c), which is designed so that the returned H2 and CO abundances are consistent with those of the underlying MHD chemical network.We use a modified version of the time-dependent gas-grain code uclchem (Holdship et al. 2017) to evolve the UMIST12 (McElroy et al. 2013) reaction network, with the same UV and cosmic ray values as the MHD simulation.The assumed elemental abundances are listed in Table 1.We take carbon, nitrogen and oxygen abundances from Sembach et al. (2000), and deplete silicon and magnesium by factors of 100 from their warm neutral medium (WNM) values, again representing depletion into refractory dust grains.Unlike most astrochemical models of star-forming clouds and cores, we do not deplete sulphur from its WNM value, which is close to the undepleted Solar value; we return to this point in Section 4.2.
We follow the chemical evolution of 10 5 tracer particles chosen from the central 16.2 pc of the computational domain at the simulation's end, which contains virtually all of the molecular material.Particles are selected randomly to evenly sample densities between 10 − 10 6 cm −3 , the highest density for which the chemical evolution is converged (Priestley et al. 2023c).Lower-density material has a negligible molecular content, and higher-density material makes up a negligible fraction of the cloud mass.We confirm that this number and selection of particles is sufficient for our purposes in Appendix A.

Radiative transfer
We perform radiative transfer modelling of our simulated cloud using radmc3d (Dullemond et al. 2012).The unstructured Voronoi mesh of the MHD simulation is interpolated onto a cubic adaptive mesh, with side length 32.4 pc and a base resoluton of 20 3 .Any cell which contains more than one Voronoi sampling point is refined into eight subcells, and this process is repeated until all cells hold at most one sampling point each.Cells are then assigned physical properties (density, velocity, gas and dust temperatures) from the nearest Voronoi sampling point, and molecular abundances from the nearest post-processed tracer particle.We assess the accuracy of the interpolation in Appendix A. The output positionposition-velocity (PPV) cubes have a spatial resolution of 0.06 pc, and a velocity resolution of 0.03 km s −1 for lines without hyperfine structure (lines with hyperfine structure are discussed below).
We use collisional rate data taken from the LAMDA (Schöier et al. 2005) and EMAA databases, for the molecules and collision partners listed in Table 2.Where collisional data for both ortho-and para-H2 are available, we assume an ortho:para ratio of 3 : 1, and we assume para-NH3 makes up half the total abundance of NH3.As our reaction network does not include separate isotopic chemistry, we assume, where necessary, that the ratio of molecular isotopologues is equal to the underlying isotope ratio, namely 12 C/ 13 C= 77 and 12 O/ 18 O= 550 (Wilson & Rood 1994).Deviations from these ratios under the conditions we investigate are unlikely to be larger than a factor of a few (Szűcs et al. 2014).We assume thick ice mantle dust opacities from Ossenkopf & Henning (1994) for wavelengths λ > 1 µm, and Mathis et al. (1983) opacities at shorter wavelengths (see Clark et al. 2012b).
The NH3, HCN and N2H + transitions involve multiple hyperfine components, and we use larger velocity resolutions of 0.25, 0.1 and 0.1 km s −1 respectively to capture all the structure without increasing the computational cost of the radiative transfer.While radmc3d does account for the optical depth effects of multiple overlapping hyperfine transitions, this is not included when determining level populations, meaning that the hyperfine levels are not radiatively coupled with each other.The distribution of flux between the individual components may therefore be calculated incorrectly if this coupling is significant, although we expect the integrated intensities over all hyperfine components to be valid.

Cloud-scale line properties
Figure 1 shows column density maps in total hydrogen nuclei of the simulated cloud, viewed both parallel ('face-on') and perpendicular ('edge-on') to the collision axis.The collision between the initial atomic clouds forms a layer of enhanced density at their interface, which subsequently fragments into a network of filaments and cores.The initial turbulent velocity field causes deviations from spherical symmetry, such as the 'bend' in the cloud when viewed edge-on.Figure 2 shows integrated intensity maps of lines from the molecules listed in Table 2, for the face-on cloud orientation (line emission maps of the edge-on cloud show the same qualitative behaviour).
The CO isotopologues, with their relatively high abundances and low critical excitation densities, show widespread emission throughout the cloud, although this is less pronounced for the rarest (C 18 O).This behaviour is also seen for several species commonly thought of as tracers of highdensity material, most notably HCN, CS and HNC, as is found observationally (Pety et al. 2017;Kauffmann et al. 2017;Evans et al. 2020).The line emission maps of NH3 and HCO + more closely track the actual column density of material shown in Figure 1, while N2H + is the only species investigated which selectively traces the densest regions of the cloud.These species' tendency to trace higher density gas than CO is due to chemical differences causing their abundances to peak at higher volume densities, whereas species such as CS and HCN are chemically quite similar to CO (Priestley et al. 2023c).
Figure 3 shows the relationships between integrated line intensity and column density, compared to observations of the Perseus molecular cloud.For all lines except NH3, the data are the pointed observations taken by Tafalla et al. (2021).For NH3, we use maps of the NGC 1333 subregion from Friesen et al. (2017), which have been re-reduced (Pineda et al. in prep.) and matched to the Herschel-derived column density 3 maps of Singh & Martin (2022).The noise level of the NH3 observations is ∼ 0.5 K, so the apparentlyconstant intensity at around this value below a column density of 10 22 cm −2 should be interpreted as an upper limit.We note that the commonly-used approach of fitting singletemperature modified blackbodies to far-IR data may introduce an observational bias in the measured column densities, but in Appendix B we demonstrate that this effect is modest (a factor of a few at the highest column density values), so we use the true column densities from the simulation throughout the rest of this paper.
The predicted relationships between column density and line intensity for the three CO isotopologues are in general agreement with those observed, although the 12 CO line is somewhat weaker than observed.A possible explanation for this difference is that our assumed radiation field strength of 1.7 Habing may be too low to properly represent the Perseus molecular cloud; Tafalla et al. (2023) find that regions with more active star formation produce more 12 CO emission at a given column density, which they attribute to higher gas temperatures caused by a stronger local UV field.The predicted strengths of the 13 CO and C 18 O lines, which originate from further within the cloud, are in better agreement with the data, suggesting that this effect, if relevant, is only important in the outer, less-shielded regions.
Alternatively, the difference in 12 CO intensities may be due to systematic differences in the velocity dispersions of our simulated cloud and Perseus.Tafalla et al. (2021) find typical 12 CO full-widths at half-maximum (FWHMs) of 3 − 4 km s −1 in Perseus, with little column density variation.Our simulated cloud also has a flat FWHM-column density relationship, but at a lower value of ∼ 1 km s −1 .If the velocity dispersion were to be increased to match the Perseus FWHM values, the increase in integrated intensity could plausibly reach the factor of 2 − 3 required to bring our simulated data into agreement with the Perseus observations in Figure 3.For the rarer CO isotopologues, the differences between observed and simulated FWHMs become less severe (2 versus 1 km s −1 for 13 CO, and 1 versus 0.7 km s −1 for C 18 O) and the line intensities are in correspondingly better agreement.This again suggests that discrepancies between our simulated cloud and the Perseus data primarily originate in the outer, lower-density regions of the cloud.
Line emission from molecules with similar chemical behaviour to CO (HCN, CS, HNC; Priestley et al. 2023c) differs from the Perseus data in our simulated observations.The intensities rise sharply with column density up to ∼ 10 22 cm −2 and then saturate at a near-constant value, whereas the observed behaviour is a continuous linear rise between 10 21−23 cm −2 .For species which preferentially trace denser gas (NH3, N2H + , HCO + ), the simulation is in somewhat better agreement with the data.As with the CO lines, this suggests that while our simulated cloud differs from Perseus in some respects along moderate-density sightlines (∼ 10 21 cm −2 ), we are accurately capturing the physical and chemical properties of the densest regions.We note that while the linear relationship between line intensity and column density has also been observed in the Orion A and California clouds (Tafalla et al. 2023), similar studies of Orion B (Pety et al. 2017;Santa-Maria et al. 2023) and W49 (Barnes et al. 2020) find plateau-like behaviour for HCN and other molecules, which more closely resembles our simulated cloud.We discuss interpretations of the unexpected linear relationship for optically-thick lines in Section 4.3.

Tracers of dense gas
Figure 4 shows maps of the intensity ratios of the 12 CO, HCN and N2H + lines seen face-on, and Figure 5 shows how their values vary with column density.Edge-on maps, as with the total intensity maps, are qualitatively similar to their faceon equivalents.The HCN/ 12 CO ratio is frequently used as an indicator of the fraction of dense gas (with the exact interpretation of 'dense' varying) on extragalactic scales, but in our simulated cloud, its value varies by only a factor of two over two orders of magnitude in column density, suggesting it is a poor indicator of genuinely-dense material.Ratios of other lines, such as CS and HCO + , also show relatively little variation over the entire range of column densities where they are detectable (above 3 × 10 21 cm −2 ).As argued by several observational studies (Pety et al. 2017;Kauffmann et al. 2017;Tafalla et al. 2021), ratios involving N2H + are much more sensitive to the presense of material above a column density of 10 22 cm −2 , thought to be related to the onset of star formation (Lada et al. 2010), although effects not considered in our simulation such as protostellar heating may complicate this interpretation (Fehér et al. 2024).
Of particular note, the N2H + /HCN ratio in our cloud rises by a factor of ten over a relatively small range in column density, making it a sensitive probe of high-density material.This line ratio has become observationally accessible on extragalactic scales in recent years, with a typical value of ∼ 0.2 (Jiménez-Donaire et al. 2023; Stuber et al. 2023).On the scale of our simulated cloud, this value corresponds closely to the 10 22 cm −2 star formation threshold found by Lada et al. (2010), suggesting that a significant fraction of the kpc-scale beam area in these extragalactic studies is made up of actively star-forming gas.
The inadequacy of HCN as a dense gas tracer is also apparent in velocity-resolved data.Figure 6 shows the face-on line profiles for 12 CO, the central HCN hyperfine component, and the isolated N2H + component, averaged over the central 16.2 pc region, as compared to the actual line-of-sight velocity distribution of the mass in this region.Both the 12 CO and HCN lines are quite accurate tracers of the whole-cloud dynamics, and in fact have almost indistinguishable profile shapes.The N2H + profile is significantly narrower, being a much better (although still imperfect) tracer of the dynamics of the high-density (> 10 4 cm −3 ) material.
The difference in the velocity dispersions measured via different tracers has been noted previously in observations (Goodman et al. 1998;Pineda et al. 2010;Ragan et al. 2015;Peretto et al. 2023), and appears to represent a real reduction in the magnitude of turbulent motions in the densest material, possibly due to dissipation in accretion shocks (Klessen et al. 2005;Priestley et al. 2023b).This means that the observational definition of velocity dispersion, used to determine quantities such as virial stability (e.g.Rigby et al. 2024) and magnetic field strength (Wang et al. 2024), is tracerdependent in a highly non-trivial manner.A structure's virial ratio determined via the HCN line will be different to the one calculated using the N2H + line, and the 'correct' value to use will depend entirely on the underlying science question.
The association between a line's velocity structure and the underlying gas dynamics should be made with considerable caution.

Core-scale line properties
To investigate the line emission properties on the scales of prestellar cores, we identify nine peaks in the face-on column density map, shown in Figure 7, and extract line profiles from 0.1 pc apertures around these peaks, approximating the effects of beam size for single-dish observations of nearby clouds (e.g. Lee et al. 1999;Tafalla et al. 2002).Figure 8  component N2H + profiles on similar spatial scales do occur in more distant infrared-dark clouds (Rigby et al. 2024).
Figure 9 shows the peak line intensity of CS versus that of N2H + for the simulated cloud core sample, compared to observational data from Lee et al. (1999) and Tafalla et al. (2002).The ratio of these lines was identified by Yin et al. (2021) as a diagnostic of the degree of magnetic support, based on simulations of isolated, initially-spherical cores.Priestley et al. (2022) found that model cores with supercritical mass-to-flux ratios (Mouschovias & Spitzer 1976), also shown in Figure 9, were incapable of matching the observed distribution of this ratio, as in the absence of magnetic support collapse proceeds too rapidly to significantly deplete CS from the gas phase, and its J = 2 − 1 line is subsequently stronger than observed.By contrast, the cores formed in our simulation agree well with the observational data, despite the initial mass-to-flux ratio of the simulation being 2.4 times the critical value (i.e.moderately supercritical).The complex formation histories of the cores in our simulation, involving material being repeatedly cycled in and out of highdensity phases before undergoing runaway gravitational collapse (Priestley et al. 2023c), results in them having very different chemical properties to equivalent models of isolated cores, which instead experience a monotonic increase in density with time.
Another problematic aspect of supercritical core models is their predicted line widths, which are inevitably much broader than the observed subsonic values due to the supersonic infall velocities of unimpeded gravitational collape (Larson 1969).Figure 9 shows the FWHMs4 of the CS and N2H + lines of the cores from our simulation and from the Priestley et al. (2022) supercritical sample, compared to observational data from Tafalla et al. (2002).In this case, forming cores self-consistently does not improve the agreement with observations: simulated cores still have much broader lines (> 0.5 km s −1 ) than real ones.Previous studies (Offner et al. 2008;Smith et al. 2012;Priestley et al. 2023b) have reproduced subsonic line widths in cores and filaments formed via supersonic converging flows, which efficiently dissipate the initial kinetic energy (Klessen et al. 2005), so the discrepancy here may be due to our initial turbulent velocity field being somewhat subsonic, even though the cloud collision itself is highly supersonic.Alternatively, simulations may still require a substantial degree of magnetic support to reproduce the observed properties of real objects, but it is clearly impossible to make this argument based on models of isolated spheres, which are far too idealised compared to the complex formation histories of cores in turbulent molecular clouds.

DISCUSSION
4.1 The (in)accuracy of isolated core models Simulations of isolated prestellar cores are widely used to interpret molecular line observations (Keto et al. 2015;Sipilä & Caselli 2018;Young et al. 2019;Sipilä et al. 2022;Tritsis et al. 2023;Redaelli et al. 2024), under the assumption that cores are sufficiently decoupled from their environment to evolve more-or-less independently.This does not appear to be the case; Wurster & Rowan (2024) find that resimulating the regions which form cores within a larger-scale cloud produces significantly different results to the original simulation, due to the chaotic nature of star formation, while several observational studies (Peretto et al. 2020;Rigby et al. 2021;Anderson et al. 2021;Redaelli et al. 2022) have suggested that cores continue to accrete mass from their surroundings throughout their lifetimes.
The prior evolutionary history of the material making up a core, and continuing accretion onto it, can significantly alter its chemical composition (Priestley et al. 2023c), with a correspondingly large effect on the predicted line emission properties.The peak line intensities of our core samples in Figure 9 resemble those of magnetically-supported isolated models from Priestley et al. (2022), despite the opposite being the case in the simulation, while the subsonic linewidths of observed cores may be a result of their formation from cloud-scale turbulent motions (Offner et al. 2008;Smith et al. 2012;Priestley et al. 2023b) rather than the isolated quasiequilibrium scenario adopted by Keto et al. (2015) and subsequent other authors.Taking these highly-idealised models as representative of real cores risks making fundamental errors in the physical interpretation of molecular line data.
We note that a similar argument could be made regarding our idealised setup of two colliding, uniform-density spheres.Real molecular clouds are rarely spherical (Faustino Vieira et al. 2024), and may be better characterised by a continuous injection of turbulence (Zhou et al. 2022), rather than the decaying velocity field we employ.However, the generallygood agreement between the predicted line emission from our simulations and observational data, on both cloud and core scales, suggests that they are accurately capturing the physical-chemical structure of real molecular clouds (see also Priestley et al. 2023a), despite their idealised initial conditions.Simulations on the galactic disc scales necessary to self-consistently follow cloud formation (e.g.Panessa et al. 2023) typically lack the resolution to investigate chemistry in detail in pre-and protostellar cores.These objects are likely to be sufficiently decoupled from their galactic environment for cloud-scale simulations to be an entirely adequate representation of their formation and evolution.

Sulphur depletion in molecular gas
It is common practice in astrochemical modelling to reduce the elemental gas-phase sulphur abundance by a factor of ten or more from its Solar value, as otherwise many sulphurbearing molecules are predicted to be far more abundant, and their line emission far stronger, than is actually observed (Fuente et al. 2019;Navarro-Almaida et al. 2020;Fuente et al. 2023).The assumption is that most of the sulphur is locked up in some chemically-unavailable solid phase, but direct measurements of sulphur depletion in the interstellar medium find a near-Solar gas-phase abundance even along the densest sightlines (Jenkins 2009), suggesting that only a modest fraction of the total sulphur is in solid form.
Our simulation, using an undepleted sulphur abundance (Table 1), successfully reproduces the observed peak intensities of CS in cores (Figure 9), often one of the more problematic species in astrochemical modelling (Navarro-Almaida et al. 2020).As argued by Hily-Blant et al. ( 2022), the main sulphur reservoir in molecular clouds may be atomic S, rather than sulphur-bearing molecules (Holdship et al. 2019;Kushwahaa et al. 2023) and ices (Boogert et al. 2015;McClure et al. 2023), which make up only a few percent of the Solar abundance when combined.It was speculated in Priestley et al. (2022) that the sulphur abundance problem may be resolved by more realistic physical models of molecular gas, rather than an improved understanding of its chemistry.We will investigate this possibility, and the detailed behaviour of sulphur-bearing species in our chemical model, in future work.

Linear intensity scaling in optically-thick lines
The linear relationship between line intensity and column density for molecules such as HCN in Perseus (Tafalla et al.  7. Green circles show observational data taken from Lee et al. (1999) and Tafalla et al. (2002).Orange triangles show the results of initially-spherical core models from Priestley et al. (2022).2021) and some other molecular clouds (Tafalla et al. 2023) is unexpected, as these lines are optically thick (typical linecentre opacities > 10 for column densities > 10 21 cm −2 ).The amount of emitting material along the line-of-sight should then have little impact on the strength of the observed emission.Explanations for this behaviour (Tafalla et al. 2021;Dame & Lada 2023) typically rely on the presence of a positive correlation between volume and column densities; if the characteristic volume density is larger along higher-column sightlines, and the level populations are subthermally-excited, then an increasing column density can produce an increase in line intensity even in the opticallythick limit.
Simulated molecular clouds typically do show such a relationship (Clark & Glover 2014;Bisbas et al. 2019), but not one capable of reproducing the linear behaviour observed by Tafalla et al. (2021).Figure 10 shows the distribution of the volume and line-of-sight column densities of arepo cells in our simulation.The correlation between column and vol-ume density is much flatter than the NH ∝ n 4/3 H proposed by Tafalla et al. (2021), and more significantly, departs severely from a one-to-one relationship; a given line-of-sight column density can comprise gas covering several orders of magnitude in volume density.It seems probable that at least the latter behaviour is shared with real molecular clouds, making the excitation temperature argument for the linear intensity scaling unviable; there should be little, if any, correlation between the excitation temperature of a transition, a complicated mass-and optical depth-weighted average over material of all densities along the line of sight, and the observed column density.We note that despite the face-on and edge-on cloud orientations having very different volume density distributions for line-of-sight columns above 10 22 cm −2 , the resulting line intensities in this regime in Figure 3 are effectively indistinguishable, suggesting that the characteristic volume density of a line of sight has a negligible effect on the resulting line emission.An alternative explanation for the increasing intensity of optically-thick lines is if the line-of-sight velocity dispersion is higher along high-density sightlines, providing more 'bandwidth' available for emission (Whitworth & Jaffa 2018).Simulations in Priestley et al. (2023a) where the velocity dispersion increases with column density produced linearlyincreasing HCN line intensities, whereas simulations without a rise in velocity dispersion resulted in an intensity plateau similar to those in Figure 3.Our simulation has an almostconstant velocity dispersion over the range of column densities where HCN is detectable (Figure 11), so no additional velocity channels become available for emission, and once the HCN line becomes optically thick, its intensity saturates.
The distinction between rising and flat velocity dispersions in Priestley et al. (2023a) was attributed to simulations of isolated versus colliding clouds.However, we do not reproduce this behaviour here.Figure 11 shows the velocity dispersion and HCN intensity relationships for a simulation where we have reduced the collision velocity from 7 to 0.5 km s −1 ; the subsonic collision velocity (also below the average velocity of the initial turbulent field) results in behaviour comparable to an isolated, turbulent cloud as studied in Priestley et al. (2023a).Rather than increasing, the velocity dispersion actually declines over the range of column densities relevant for HCN emission, resulting in line intensities which are indistinguishable from the colliding case.We attribute this to our use of a lower initial density compared to Priestley et al. (2023a) (10 versus 250 cm −3 ), so that the initial turbulent velocity field is dissipated before the cloud material has had time to collapse to the higher densities where HCN emission can be produced.Establishing the relevance of this distinction to the interpretation of observational data is left to future work.

CONCLUSIONS
We have combined an MHD simulation of the formation and dynamical evolution of a molecular cloud with timedependent chemistry and radiative transfer modelling, producing self-consistent synthetic observations in the J = 1 − 0 rotational transitions of 12 CO, 13 CO, C 18 O, HCN, N2H + , HCO + and HNC, the J = 2 − 1 transition of CS, and the (1, 1) inversion transition of para-NH3.Our main results are as follows: • Our simulated cloud reproduces the observed range of intensities for most lines well, with better agreement for species which preferentially trace high-density regions.This suggests that despite the idealised initial conditions and modelling assumptions, the simulation has produced a molecular cloud with similar structural properties to real ones.
• HCN, along with many other molecules often described as dense-gas tracers, is actually a poor tracer of dense gas, being detectable down to low column density.The velocity information derived from the HCN line profile mostly traces the bulk dynamics of the cloud.N2H + is a much more selective tracer of high-density (> 10 22 cm −2 ) regions, and its line profile shape is strongly correlated with the dynamics of material above 10 4 cm −3 .
• Cores formed from the large-scale cloud dynamics in the simulation agree well with the distribution of line intensities from observational samples, although their line widths are somewhat broader.They do not agree with the predictions of an equivalent sample of isolated core models, which neglect the formation of the cores themselves from their more diffuse environment.Using these isolated core models to interpret molecular line observations is likely to result in erroneous conclusions.
• We obtain CS line intensities in good agreement with observed values using an undepleted elemental abundance of sulphur.This suggests that the common modelling practice of reducing this quantity by several orders of magnitude to match observations of sulphur-bearing species may be a sign of the limitations of the physical models employed, rather than genuine evidence for sulphur depletion in molecular clouds.
Our work demonstrates the importance of simulating the large-scale dynamics of molecular clouds when investigating their line emission properties, even when the subject of interest is on much smaller scales.This serves as a starting point for future work to fully exploit the ability of molecular line emission to advance our understanding of star formation.

APPENDIX A: RESOLUTION TESTS
The MHD, chemical and radiative transfer simulations in this paper each have mutually-incompatible concepts of resolution, which results in some losses when properties are mapped between them.The cubic adaptive mesh used by radmc3d loses detail when compared to the unstructured Voronoi mesh used by arepo, and with only 10 5 post-processed tracer particles compared to several million grid cells, the information on molecular abundances is inevitably 'smeared out' to some extent.We demonstrate below that neither of these effects is large enough to substantially alter our results.
Figure A1 shows the relationship between the true column density obtained from the arepo Voronoi mesh, and the column density from the radmc3d grid used for the radiative transfer.We determine the latter by fixing the dust temperature to 10 K in all cells, and producing an image of the cloud at 850 µm, where the dust optical depths are ≲ 10 −4 .As the dust opacity and gas-to-dust ratio are both constant, the intensity is directly proportional to the column density of the radmc3d grid.After converting from intensity to NH, the radmc3d column densities are almost identical to the true values, with the exception of the highest end of the range near 10 23 cm −2 .Even here, the difference is less than a factor of two, and these column densities represent a small fraction of the cloud area, so we consider the mapping between the arepo and radmc3d grids to be robust.
While the radmc3d grid appears to accurately represent the underlying structure of the MHD simulation, the molecular abundances used in the radiative transfer are taken from 10 5 post-proessed tracer particles, compared to ∼ 20 million cells in the radmc3d grid.We assess the accuracy of this approach using CO emission lines, as the CO abundance for every cell can be taken directly from the arepo simulation as a point of comparison to the NEATH tracer particle abundance.
Figure A2 shows the relationship between integrated intensity of the three CO isopologues for the two abundance determinations.All three lines show a linear correlation in intensity between the NEATH and arepo values, albeit with some scatter.The two rarer isotopologues tend to have somewhat lower line intensities for the NEATH abundances than when using the arepo values, but this can be attributed to the neglect of freeze-out in the arepo chemical model, which severely reduces the CO abundance at high volume densities (Priestley et al. 2023c).This effect is not visible for the moreabundant 12 CO line, which is too optically-thick to probe the densities where freeze-out occurs.The otherwise-good correlation for all isotopologues suggests that the number of tracer particles is adequate to correctly capture the molecular structure of the cloud.

APPENDIX B: DUST-DERIVED COLUMN DENSITIES
Figure B1 shows the column density of the cloud inferred from thermal dust emission, the most common means of assessing molecular cloud structures observationally.We produce images of the cloud at wavelengths of 100, 160, 250, 350 and 500 µm, corresponding to the central wavelengths of the Herschel filters generally used to study cold molecular gas (e.g.Ragan et al. 2012a;Könyves et al. 2015).We fit the resulting pixel-by-pixel spectral energy distributions (SEDs) with a single-temperature modified blackbody, using the same dust opacity as in the radiative transfer calculation, and the same dust-to-gas ratio to convert dust mass into total column density.Any differences between the two column densities are

Figure 1 .
Figure 1.Column density maps of the cloud seen face-on (left) and edge-on (right).

Figure 2 .
Figure 2. Integrated line intensity maps for the cloud seen face-on.

Figure 3 .
Figure 3. Integrated line intensity versus column density for the cloud seen face-on (blue) and edge-on (red).Crosses show the median pixel values, error bars the 16th/84th percentiles.Observational data from the Perseus molecular cloud (Friesen et al. 2017; Tafalla et al. 2021) are shown in black.

Figure 4 .Figure 5 .Figure 6 .
Figure 4. Maps of the ratio of integrated line intensities for HCN and 12 CO (left), N 2 H + and 12 CO (centre), and N 2 H + and HCN (right), for the cloud seen face-on.Pixels with integrated intensities of either line below 0.2 K km s −1 have been masked.

Figure 7 .
Figure 7. Positions of cores selected from the face-on column density map, highlighted with white circles.Numbers correspond to the labelling in Figure 8.For the purpose of clarity, the circles are three times larger than the 0.1 pc extraction regions.

Figure 8 .Figure 9 .
Figure 8. Line profiles of C 18 O (blue), HCN (red) and N 2 H + (green) for the cores identified from the simulated cloud in Figure 7.Text labels indicate the column density of the 0.1 pc extraction region.

Figure 10 .Figure 11 .
Figure 10.Distribution of arepo cells by volume density and column density as seen face-on (left) and edge-on (right).The dashed black lines show the relationship proposed by Tafalla et al. (2021) to explain linear correlations between line intensities and column density.

Figure A1 .
Figure A1.Column density of the radmc3d grid versus the original arepo Voronoi mesh.Crosses show the median pixel values, error bars the 16th/84th percentiles.The dashed red line shows the exone-to-one relationship.

Table 1 .
Elemental abundances used in the chemical modelling.

Table 2 .
Molecules investigated, collisional partners, and sources of collisional rate data.