# Non-Equilibrium Abundances Treated Holistically (NEATH): the molecular
composition of star-forming clouds

F. D. Priestley1, P. C. Clark1, S. C. O. Glover2, S. E. Ragan1, O. Fehér1, L.
R. Prole1, R. S. Klessen2,3
1School of Physics and Astronomy, Cardiff University, Queen’s Buildings, The
Parade, Cardiff CF24 3AA, UK
2Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische
Astrophysik, Albert-Ueberle-Straße 2, D-69120 Heidelberg, Germany
3Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches
Rechnen, Im Neuenheimer Feld 205, D-69120 Heidelberg, Germany
Email<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Much of what we know about molecular clouds, and by extension star formation,
comes from molecular line observations. Interpreting these correctly requires
knowledge of the underlying molecular abundances. Simulations of molecular
clouds typically only model species that are important for the gas
thermodynamics, which tend to be poor tracers of the denser material where
stars form. We construct a framework for post-processing these simulations
with a full time-dependent chemical network, allowing us to model the
behaviour of observationally-important species not present in the reduced
network used for the thermodynamics. We use this to investigate the chemical
evolution of molecular gas under realistic physical conditions. We find that
molecules can be divided into those which reach peak abundances at moderate
densities ($10^{3}\,{\rm cm}^{-3}$) and decline sharply thereafter (such as CO
and HCN), and those which peak at higher densities and then remain roughly
constant (e.g. NH3, N2H+). Evolving the chemistry with physical properties
held constant at their final values results in a significant overestimation of
gas-phase abundances for all molecules, and does not capture the drastic
variations in abundance caused by different evolutionary histories. The
dynamical evolution of molecular gas cannot be neglected when modelling its
chemistry.

###### keywords:

astrochemistry – stars: formation – ISM: molecules – ISM: clouds

††pubyear: 2023††pagerange: Non-Equilibrium Abundances Treated Holistically
(NEATH): the molecular composition of star-forming clouds–3

##  1 Introduction

Molecular clouds, as the name suggests, are primarily made up of molecules,
the most common being molecular hydrogen (H2) and carbon monoxide (CO). As the
birthplaces of stars (Bergin & Tafalla, 2007), line emission from molecular
rotational transitions with energies low enough to be excited at the $\sim
10\,{\rm K}$ temperatures of molecular clouds contains key information about
the process of star formation. Interpreting this data, or producing synthetic
observations which can be compared to it, requires knowedge of the abundances
of the relevant species. While these can be parametrized as simple functions
of the local properties (e.g. Tafalla et al., 2002; Smith et al., 2012, 2013;
Jones et al., 2023), the comparable chemical and dynamical timescales (Banerji
et al., 2009) mean that molecular abundances are dependent on the prior
history of the gas. A detailed treatment of non-equilibrium chemistry in
molecular clouds is therefore necessary in order to advance our understanding
of star formation.

Modern hydrodynamical simulations of molecular clouds (e.g. Clark et al.,
2012b; Walch et al., 2015; Smith et al., 2020) typically include time-
dependent CO chemistry self-consistently, albeit in a simplified fashion, due
to its importance as a coolant and its effect on the abundances of other major
coolants (C, C+, O). However, its high molecular abundance, and the low
critical density of its key rotational lines, make CO a poor tracer of the
high-density regions of clouds where star formation actually occurs (e.g.
Clark et al., 2019; Priestley et al., 2023b). In cloud-scale simulations, the
abundances of molecules tracing denser gas are usually inserted by hand, based
on either observational estimates or on smaller-scale simulations of
prestellar cores (Smith et al., 2012; Jones et al., 2023). These core-scale
simulations (e.g. Aikawa et al., 2005; Sipilä & Caselli, 2018) make
assumptions such as isothermality and/or spherical symmetry which are unlikely
to apply in reality, and in addition are highly sensitive to the (unknown)
properties of the cloud material exterior to the core (Kaminski et al., 2014;
Priestley et al., 2023a; Jensen et al., 2023). It is therefore desirable to
follow the chemical evolution of star-forming material from the earliest
stages of molecular cloud formation.

The abundances of molecular species other than CO (an important coolant)
typically have little effect111With the possible exception of their
contribution to the ionization state of the gas, when non-ideal
magnetohydrodynamic effects are considered (Wurster, 2021). on the thermal or
dynamical evolution of molecular gas. It is thus possible to decouple the
detailed chemical evolution from the hydrodynamics and separately post-process
this stage (e.g. Panessa et al., 2023) with very little loss of accuracy
(Ferrada-Chamorro et al., 2021), rather than taking the computationally-
expensive approach of including hundreds of molecular species directly into
the underlying hydrodynamic simulation (e.g. Lupi et al., 2021). In this
paper, we develop a framework to do this while maintaining consistency between
the chemical network incorporated in the hydrodynamics (designed for accuracy
and efficiency at low to moderate densities) and the more complex one used in
the post-processing stage (designed for completeness, but poorly constrained
at temperatures above $\sim 20\,{\rm K}$): Non-Equilibrium Abundances Treated
Holistically (NEATH).222Neath is a town in Wales roughly halfway between
Cardiff and Swansea. We use this to investigate the molecular composition of
star-forming clouds, and how a full time-dependent treatment of the physical
and chemical evolution differs from simplified approaches common in the
literature.

Table 1: Initial gas-phase elemental abundances, relative to hydrogen nuclei, used in the chemical modelling. Element | Abundance |  | Element | Abundance
---|---|---|---|---
C | $1.4\times 10^{-4}$ |  | S | $1.2\times 10^{-5}$
N | $7.6\times 10^{-5}$ |  | Si | $1.5\times 10^{-7}$
O | $3.2\times 10^{-4}$ |  | Mg | $1.4\times 10^{-7}$

Figure 1: Distribution of gas temperature (left) and shielding column density
(right) versus volume density, with the colour scale representing the number
of tracer particles. Power law relationships which reproduce the general trend
of the data (Equations 4 and 5) are shown as dashed black lines.

Figure 2: Average abundances versus gas density of H2 (left) and CO (right).
Values from arepo using the Gong et al. (2017) network are shown in red, the
NEATH values in blue. Crosses show the median values, with the 16th and 84th
percentiles as error bars.

Figure 3: Evolution of the gas density (left) and the CO (middle) and NH3
(right) abundances for the shock test, modelled with a time resolution of
$10\,{\rm kyr}$ (solid black lines) and $50\,{\rm kyr}$ (dashed red lines).
Figure 4: Median CO abundance versus density at a time resolution of
$8.8\,{\rm kyr}$ (blue crosses), $44\,{\rm kyr}$ (green squares), and
$88\,{\rm kyr}$ (red circles).

##  2 Method

###  2.1 Magnetohydrodynamical model

We simulate the dynamical evolution of molecular clouds using a modified
version of arepo (Springel, 2010), a moving-mesh magnetohydrodynamics (MHD)
code. The thermodynamics of the gas and dust are treated according to Glover &
Mac Low (2007) and Glover & Clark (2012); this represents a largely complete
account of all important processes, including radiative heating and cooling,
cosmic rays, dust-gas thermal coupling, and the energy liberated by exothermic
chemical reactions. The abundances of important gas-phase coolants (H2, C, C+,
O, CO) are calculated on-the-fly using a time-dependent chemical network based
on that of Gong et al. (2017), as described in Hunter et al. (2023). We assume
a background ultraviolet (UV) radiation field with the Mathis et al. (1983)
spectral energy distribution. Shielding from the UV field is also done on-the-
fly; column densities of H2, CO, and C for self-shielding, and the total gas
column responsible for dust extinction, are calculated using the treecol
algorithm (Clark et al., 2012a).

The on-the-fly chemical network in our arepo simulations contains only a
limited subsample of the astrophysically-relevant reactions involving key
coolant species, and neglects freeze-out of CO onto dust grains entirely.
However, benchmarking against more comprehensive networks has shown that the
effect of the missing reactions on the abundances of coolants is small (Gong
et al., 2017). Additionally, at the densities where freeze-out becomes
important ($\gtrsim 10^{4}\,{\rm cm}^{-3}$), the reduced gas-phase CO
abundance has little effect on the temperature (Goldsmith, 2001), due to a
combination of highly optically-thick emission and the increasing importance
of dust cooling. We therefore expect that the on-the-fly chemical network is
sufficient to accurately capture the thermal evolution of the gas.

The arepo grid is refined with a target mass per cell of $0.005\,{\rm
M}_{\odot}$, with the constraints that the volumes of neighbouring cells do
not differ by more than a factor of eight, and enforced minimum (maximum) cell
volumes of $5\times 10^{-11}\,{\rm pc}^{3}$ ($12\,{\rm pc}^{3}$). We also
ensure that the local Jeans length is resolved by at least 16 cells. We
include a sink particle algorithm to treat small regions collapsing under
self-gravity, as described in Tress et al. (2020) and Prole et al. (2022);
cells with a density above $2\times 10^{-16}\,{\rm g\,cm}^{-3}$ are converted
into sink particles if they are at a local minimum of the gravitational
potential, and the gas within a sink formation radius of $9\times
10^{-4}\,{\rm pc}$ is (by mass-weighted average) both moving and accelerating
towards the prospective sink particle.

###  2.2 Initial conditions

Our simulations are of spherical, colliding clouds, each with mass
$M=10^{4}\,{\rm M}_{\odot}$ and radius $R=19\,{\rm pc}$, for an initial number
density $n_{\rm H}=10\,{\rm cm}^{-3}$. The centres of the clouds are initially
separated by $2R$, so the clouds are just touching. The initial gas and dust
temperatures are set to $300\,{\rm K}$ and $15\,{\rm K}$ respectively. The
clouds are seeded with a virialised turbulent velocity field. The
computational domain is a cubic box of side length $96\,{\rm pc}$ with
periodic boundary conditions, filled with a hot, rarefied medium ($0.01\,{\rm
cm}^{-3}$, $10^{4}\,{\rm K}$). The initial velocity of each cloud is $7\,{\rm
km\,s^{-1}}$, giving a collision velocity of $14\,{\rm km\,s^{-1}}$, and we
include a $3\,{\rm\mu G}$ magnetic field parallel to the collision axis.
Simulations are run for $5.0\,{\rm Myr}$, by which point the clouds have
formed $28$ sink particles with a total mass of $24.2\,{\rm M}_{\odot}$.
Chemical parameters which are shared with the later post-processing stage
described in Section 2.3 are the elemental carbon ($1.4\times 10^{-4}$) and
oxygen ($3.2\times 10^{-4}$) abundances from Sembach et al. (2000), the
external ultraviolet (UV) radiation field ($1.7$ times the Habing (1968)
field), and the cosmic ray ionisation rate ($10^{-16}\,{\rm s^{-1}}$ primary
ionisations of H2). The ‘metal’ abundance is set to $10^{-7}$, corresponding
to the high levels of depletion seen in the dense ISM (Jenkins, 2009).

###  2.3 Chemical post-processing

We use tracer particles, as implemented in arepo by Genel et al. (2013), to
record the evolving properties of individual parcels of gas, which are used as
input for uclchem (Holdship et al., 2017), a time-dependent gas-grain chemical
code.333The source code for the version used here is available from
https://fpriestley.github.io/neath/. Each cell is initialised with five tracer
particles, which are subsequently propagated in a stochastic manner weighted
by the mass fluxes between cells. The properties of the cell the particles
find themselves in are recorded at regular intervals; any changes on shorter
timescales than this interval will not be accounted for in the chemical post-
processing. Increasing the time resolution of the model by decreasing the
interval between updates will produce more accurate results, at the expense of
increasing the memory requirements and so (for a fixed memory allocation)
reducing the number of tracers which can be post-processed. We investigate the
sensitivity of our model to the chosen time resolution in Section 2.4, and for
the remainder of the paper use a value of $44\,{\rm kyr}$, which we find
provides an acceptable trade-off between accuracy and computational cost (the
chemical models described below have a runtime of $\sim 200$ CPU-hours, and
show excellent parallel scaling, compared to $\sim 6000$ CPU-hours to run the
underlying MHD simulation).

The UV field strength and cosmic ray ionisation rate are kept the same as in
the arepo models. We use the UMIST12 (McElroy et al., 2013) chemical network,
which includes additional elements not present in that of Gong et al. (2017).
The abundances of these elements, listed in Table 1, are also taken from
Sembach et al. (2000), but we reduce the values of the refractory elements (Mg
and Si) by a factor of 100 so that the level of depletion is consistent with
the ‘metal’ abundance in arepo. Although astrochemical models often also
deplete sulphur in order to match observations of S-bearing species (e.g.
Navarro-Almaida et al., 2020), we leave its abundance at the Sembach et al.
(2000) value; there is no known solid-phase reservoir that can account for any
significant level of sulphur depletion (Jenkins, 2009; Boogert et al., 2015),
and recent work suggests that the gas-phase abundance may in fact be close to
its undepleted value (Hily-Blant et al., 2022).

The quantities required to calculate the chemical evolution of a particle are
the number density of hydrogen nuclei ($n_{\rm H}$), the gas
temperature444uclchem assumes that the gas and dust temperatures are equal.
This is not necessarily the case, particularly at low densities, but the dust
temperature is only used to calculate the rates of diffusion reactions on
grain surfaces. These become significant only at relatively high densities
($>10^{4}\,{\rm cm}^{-3}$), where dust and gas temperatures are close to
equal, so the distinction is unlikely to have much impact on our present
results. We plan to extend the model to follow gas and dust temperatures
separately in future work. ($T$), and the effective shielding column densities
(distinct from the line-of-sight columns; Clark & Glover 2014) of hydrogen
nuclei ($N_{\rm H}^{\rm shield}$), molecular hydrogen ($N_{\rm H_{2}}^{\rm
shield}$), and CO ($N_{\rm CO}^{\rm shield}$). For $N_{\rm H}^{\rm shield}$,
we use the same approach as in arepo, and take a weighted average over the
$n_{\rm PIX}$ treecol rays,

$\exp\left(-2.5\,A_{\rm V}\right)=\frac{1}{n_{\rm
PIX}}\displaystyle\sum_{i}^{n_{\rm PIX}}\exp\left(-2.5\,A_{\rm V,i}\right),$
(1)

where $A_{\rm V}=5.3\times 10^{-22}\left(N_{\rm H}^{\rm shield}/\,{\rm
cm}^{-2}\right)\,{\rm mag}$ (Bohlin et al., 1978). The factor of $2.5$ is
representative of the key photodissociation reactions in the network (Glover &
Mac Low, 2007).

The self-shielding of H2 and CO is more complicated, as the lines involved may
saturate, and broadening and Doppler shifting of the frequencies can be
important (Glover & Mac Low, 2007). While arepo calculates the self-shielding
factor for each ray individually, we find that the harmonic mean of the column
density of the relevant molecule,

$(N_{X}^{\rm shield})^{-1}=\frac{1}{n_{\rm PIX}}\displaystyle\sum_{i}^{n_{\rm
PIX}}(N_{X,i}^{\rm shield})^{-1},$ (2)

provides an acceptable approximation. Self-shielding factors are then
calculated following Federman et al. (1979) for H2 and van Dishoeck & Black
(1988) for CO. We also replace the default formation rate of H2 on grain
surfaces with

$R_{\rm H_{2}}=3\times 10^{-18}\,\sqrt{T}\,\exp\left(-\frac{T}{1000\,{\rm
K}}\right)\,{\rm cm^{3}\,s^{-1}}$ (3)

from de Jong (1977), to suppress its formation at high temperatures.

We randomly select $10^{5}$ tracer particles to post-process, which at the
simulation endpoint are located within $16.2\,{\rm pc}$ of the centre of the
computational domain, distributed so that we evenly sample555We achieve this
by randomly selecting $5000$ particles from each of $20$ logarithmically-
spaced density bins. The distribution of particles over these density bins is
not uniform, leading to the ‘striping’ visible in Figure 1, where particles
are more densely sampled toward one edge of a bin than the other. densities in
the range $10-10^{7}\,{\rm cm}^{-3}$. Figure 1 shows the distributions of $T$
and $N_{\rm H}^{\rm shield}$ with $n_{\rm H}$ for these particles; the chosen
density range corresponds to temperatures $\sim 5-200\,{\rm K}$ and visual
extinctions $\sim 0.1-30$ mag. The general trend of the temperature
distribution can be reproduced by the equation of state

$T=120\left(n_{\rm H}/10\,{\rm cm}^{-3}\right)^{-0.8}+10\left(n_{\rm
H}/10^{5}\,{\rm cm}^{-3}\right)^{-0.15}\,{\rm K},$ (4)

although this fails to capture the plateau at $\sim 20\,{\rm K}$ between
$10^{3}-10^{4}\,{\rm cm}^{-3}$. Similarly, the typical column density is
correlated with density as

$N_{\rm H}^{\rm shield}=6\times 10^{20}\left(n_{\rm H}/10\,{\rm
cm}^{-3}\right)^{0.3}\,{\rm cm}^{-2},$ (5)

but with $\sim 0.5$ dex of scatter between particles of the same density. This
is close to the $N_{\rm H}^{\rm shield}\propto n_{\rm H}^{0.5}$ scaling that
would result from shielding by a Jeans length of material at the same volume
density (as adopted in Priestley et al. 2023b); the slightly lower exponent
likely represents the fact that the average density within a Jeans length of a
point tends to be lower than the density at that point.

Figure 2 shows how the post-processed abundances of H2 and CO compare to the
internal arepo abundances. The H2 abundances predicted by NEATH are consistent
with the arepo values over the entire $10-10^{7}\,{\rm cm}^{-3}$ density
range; while the values for individual particles are not always identical, the
differences are comparable to or smaller than the scatter between different
particles of the same density, and the overall trend of H2 abundance with
density is reproduced. The CO abundances are similarly consistent between
NEATH and arepo between $10^{2}-10^{4}\,{\rm cm}^{-3}$. At lower densities,
NEATH appears to underpredict the amount of CO, likely because our self-
shielding column densities from Equation 2 are underestimating the true
‘effective’ value. At higher densities, NEATH again predicts lower CO
abundances than arepo, with increasing severity as the density increases, but
in this regime the difference is due to the depletion of CO onto grain
surfaces. Freeze-out is not modelled in the Gong et al. (2017) network, so
this behaviour is expected, and the post-processed NEATH abundances are
overall in good agreement with those in the underlying MHD simulation.

###  2.4 Effects of time resolution

We investigate the sensitivity of the chemistry to unresolved changes in the
gas properties with a simplified model of a transient density enhancement
(hereafter ‘shock’, although we note that this is different from genuine
shocks such as in protostellar outflows), treated as a temporary and
instantaneous change in the gas density. The temperature and column density
are assumed to follow the relationships given by Equations 4 and 5
respectively. Column densities for self-shielding are obtained by multiplying
Equation 5 by the local abundance of the relevant molecule. The initial
density is $10^{3}\,{\rm cm}^{-3}$, which is enhanced to $10^{5}\,{\rm
cm}^{-3}$ after $0.515\,{\rm Myr}$ of evolution (sufficient to establish
chemical near-equilbrium) for a duration of $40\,{\rm kyr}$, after which it
reverts to the original value and evolution continues up to an end point of
$1\,{\rm Myr}$. We perform the chemical modelling with timesteps of $10\,{\rm
kyr}$ and $50\,{\rm kyr}$; in the latter case, the shock is missed completely,
and the chemical model assumes all physical properties remain constant.

Figure 3 shows the evolution of the gas density and the abundances of CO and
NH3, chosen for their different chemical behaviour (see Section 3.2). With a
time resolution of $50\,{\rm kyr}$, the transient density enhancement is not
captured, and the abundances remain nearly constant (there is a slight
increase for both molecules as they have not quite reached their equilibrium
values). CO is barely affected by the shock; the increased density causes a
slight decrease (around a factor of two) due to depletion, but the duration is
not long enough for this to have much effect, and once the initial density is
restored the CO abundance returns to its original value within a few
$100\,{\rm kyr}$.

The shock has a much greater impact on the NH3 abundance, which increases
abruptly by several orders of magnitude during the density enhancement.
However, it declines similarly abruptly once the density enhancement ends.
While the abundance does not return to its pre-shock value over the $\,{\rm
Myr}$ duration of the model, the difference is only a factor of a few compared
to the case where the shock was not captured at all. These differences are
negligible compared to those caused by different tracer particle histories
and/or the current physical conditions. Temporary changes in the gas
properties on timescales too short to be captured are unlikely to have any
significant impact on our results.

While unresolved transient changes do not appear to be an issue, poor
resolution of rapid evolution is more problematic. Figure 4 shows the median
CO abundance versus density for a smaller sample of $10^{4}$ tracer particles,
with the time resolution ranging from $8.8$ to $88\,{\rm kyr}$. The chemical
model is well-converged up to densities of $\sim 10^{4}\,{\rm cm}^{-3}$, but
at higher densities the gas-phase CO abundance decreases as the time
resolution increases. Much of the material in this density regime is
gravitationally collapsing, and so evolves on approximately the freefall
timescale $t_{\rm ff}\sim 0.8\,(n_{\rm H}/10^{4}\,{\rm cm}^{-3})^{-0.5}\,{\rm
Myr}$. If the interval between updates of physical properties does not fully
resolve the freefall time, the chemistry spends too long at too low a density,
and so undergoes less freeze-out than it should, leading to the higher CO
abundances in the lower resolution models.

This effect becomes increasingly significant as the density increases and the
freefall timescales decreases; while our fiducial $44\,{\rm kyr}$ timestep is
converged up to a density of $10^{5}\,{\rm cm}^{-3}$, this is clearly not the
case at $10^{7}\,{\rm cm}^{-3}$. Ferrada-Chamorro et al. (2021) suggest a
chemical timestep below $0.01\,t_{\rm ff}$ is required for percent-level
accuracy in the molecular abundances, which corresponds to $0.3\,{\rm kyr}$ at
a number density of $10^{7}\,{\rm cm}^{-3}$. However, increasing the time
resolution by a factor of $5$ from our fiducial value only results in a factor
of a few decrease in the CO abundance at densities up to $10^{6}\,{\rm
cm}^{-3}$, suggesting that we are reasonably close to convergence in this
regime, and the mass fraction of gas with higher densities is extremely small
(about $10^{-4}$). Reaction rates in the chemical network are rarely, if ever,
known to better precision than a factor of a few (e.g. Rocha et al., 2023), so
a percent-level numerical accuracy in abundance would be dwarfed by systematic
uncertainties. Based on the results in Figure 4, we suggest that a time
resolution of $\sim 0.3\,t_{\rm ff}$ is sufficient to obtain abundances
accurate to within a factor of two. The fiducial time resolution of $44\,{\rm
kyr}$ is likely to be adequate for investigating the overall chemical
composition of the clouds, being effectively converged between number
densities of $10^{2}-10^{6}\,{\rm cm}^{-3}$.

Panessa et al. (2023) deal with the issue of poorly-resolved changes in
physical properties by ‘subcycling’; using multiple shorter timesteps for the
chemical evolution and interpolating between the initial and final physical
properties if the difference between these is above a certain threshold. With
this approach, it is not clear how exactly the interpolation should be done.
Panessa et al. (2023) use a linear interpolation, but the gas density of a
gravitationally collapsing region grows exponentially (Larson, 1969), so a
linear interpolation will tend to overestimate the amount of time spent at
high density and thus the extent to which molecules are depleted. This is
unlikely to be an issue in their case, as their simulations do not extent far
beyond a density of $10^{4}\,{\rm cm}^{-3}$, but the method is not easily
applicable to our smaller-scale, higher-density simulations.

Figure 5: Average abundances versus gas density of CO (left), HCN (centre) and
N2H+ (right), for dynamic (blue) and static (red) models. Crosses show the
median values, with the 16th and 84th percentiles as error bars.

Figure 6: Average value of the overdensity parameter (left), gas temperature
(centre) and shielding column density (right), as a function of gas density
and CO abundance. Figure 7: Column densities in the $x-y$ plane of total
hydrogen nuclei and H2 (top panel), and of CO, NH3, HCN, N2H+, HCO+, CS, HNC,
CH3OH, CN and C2H (bottom panel). Note that all post-processed tracer
particles are located within a sphere of radius $16.2\,{\rm pc}$, so the
abundances and column densities of all molecules are set to zero outside of
this region. Figure 8: Box plot of CO column density versus $N_{\rm H}$. Boxes
show the median and the 25th/75th percentiles; whiskers show the 10th and 90th
percentiles. The dashed blue line shows the column density expected if all
carbon is in the form of CO.

Figure 9: Box plots of molecular column density versus $N_{\rm H}$ for; NH3,
HCN and N2H+ (top row); HCO+, CS and HNC (middle row); CH3OH, CN and C2H
(bottom row). Boxes show the median and the 25th/75th percentiles; whiskers
show the 10th and 90th percentiles. Figure 10: Box plot of CO abundance versus
gas density. Boxes show the median and the 25th/75th percentiles; whiskers
show the 10th and 90th percentiles. The dashed blue line shows a fifth-order
polynomial fit to the median abundances.

Figure 11: Box plots of molecular abundance versus gas density for; NH3, HCN
and N2H+ (top row); HCO+, CS and HNC (middle row); CH3OH, CN and C2H (bottom
row). Boxes show the median and the 25th/75th percentiles; whiskers show the
10th and 90th percentiles. The dashed blue lines show fifth-order polynomial
fits to the median abundances. For HCN, N2H+, HCO+ and CS, the abundance-
density relationships assumed by Smith et al. (2012, 2013) are shown as dashed
red lines.

Figure 12: Time evolution of the gas density (left) and the CO (middle) and
N2H+ (right) abundances for five tracer particles with a final density of
$\sim 10^{5}\,{\rm cm}^{-3}$. Solid lines show the abundances for a dynamic
physical model, dotted lines for a static model with physical properties fixed
at their final values.

Figure 13: Average abundances of CO (left), HCN (middle) and N2H+ (right)
versus average density, within a radius of $0.05\,{\rm pc}$ of each of the
$28$ sink particles at $5.0\,{\rm Myr}$.

##  3 Results

###  3.1 Dynamic versus static chemical evolution

The impact of the dynamical evolution of the gas on the time-dependent
chemical evolution is shown in Figure 5 for CO, HCN and N2H+, comparing the
molecular abundances calculated with evolving gas properties (hereafter
‘dynamic’) to values calculated by evolving the chemistry for the same amount
of time ($5.0\,{\rm Myr}$), but holding all physical properties constant at
their final values (hereafter ‘static’). At low to moderate densities
($\lesssim 10^{4}\,{\rm cm}^{-3}$), the differences are minor, with static
chemistry tending to slightly overestimate the abundances of all species.
Under these conditions, the chemistry can react on short enough timescales to
effectively adopt the ‘appropriate’ abundances for the relevant physical
conditions (Holdship & Viti, 2022), although there is more scatter in
abundance at fixed density in the dynamic case due to the different tracer
particle histories.

At higher densities ($\gtrsim 10^{4}\,{\rm cm}^{-3}$), the static and dynamic
abundances diverge significantly. Holding the physical properties constant for
several $\,{\rm Myr}$ allows the gas-phase abundances to reach equilibrium
between freeze-out onto grain surfaces and desorption of molecules from ice
mantles back into the gas-phase. As such, the abundances of all molecules
decrease monotonically with density, because the rate of depletion increases
as $n_{\rm H}^{2}$ compared to $n_{\rm H}$ for desorption processes. There is
essentially no scatter in the abundances at fixed density in this regime for
the static model; differences in the other gas properties have a negligible
impact on the chemistry.

The dynamic abundances show similar qualitative behaviour to the static ones,
but with a shallower decline which begins to plateau towards a constant
average abundance at densities above $\sim 10^{5}\,{\rm cm}^{-3}$ (although to
some extent this may be due to insufficient time resolution; Figure 4).
Particles spend only a very small fraction of their $5.0\,{\rm Myr}$ evolution
at these high densities, so there is insufficient time for an equilibrium
between freeze-out and desorption to be established. This also means the
abundances of a given tracer particle are strongly dependent on its history,
rather than being mostly established by its current properties as at lower
densities. At fixed density, the HCN abundance can vary by up to three orders
of magnitude, compared to almost no scatter predicted by the static model.
While the static abundances for HCN are within the scatter of the dynamic
values, they are located towards the extreme lower end of the distribution, so
are still poorly representative of the typical chemical conditions. For CO and
N2H+, the static abundances are significantly discrepant with the range of
dynamic values, even at densities where the chemical evolution is well-
converged.

The scatter in abundance at fixed density in the dynamic results is not
random; it is a deterministic result of the different pathways taken by
different tracer particles to reach the same final state. To investigate the
impact of evolutionary history on the distribution of abundances at fixed
density, we introduce the ‘overdensity’ parameter,

$\Gamma=\frac{1}{n_{\rm H}^{2}(0)\,t_{\rm end}}\int^{t_{\rm end}}_{0}n_{\rm
H}^{2}(t)\,{\rm d}t,$ (6)

quantifying how long each particle has spent at a higher density than the
original $10\,{\rm cm}^{-3}$ (the $n_{\rm H}^{2}$ weighting corresponds to the
density dependence of two-body reactions). Figure 6 shows how the scatter in
CO abundance at fixed density is affected by $\Gamma$, $T$ and $N_{\rm H}^{\rm
shield}$. Below a density of $10^{4}\,{\rm cm}^{-3}$, more overdense particles
tend to have higher abundances, as the higher densities are more favourable
for forming CO. The trend reverses above $10^{4}\,{\rm cm}^{-3}$, with more
overdense particles having spent longer under conditions favourable to freeze-
out, and so having lower gas-phase CO abundance.

In contrast to $\Gamma$, the CO abundance has almost no dependence on the gas
temperature at fixed density. There is a mild dependence on shielding column
density at lower densities, but by the point where the abundance peaks ($\sim
10^{3}\,{\rm cm}^{-3}$), the gas is typically well-shielded enough to make any
variations in the attenuated UV flux unimportant. The prior history of the gas
is much more important for determining its chemical composition than its
present-day physical properties (other than density). The several orders of
magnitude of scatter in abundance at fixed density, which is not present in
the static model, demonstrates that the assumption of unchanging physical
properties is completely inappropriate, particularly at the highest densities.
The physical and chemical evolution of the gas have to be considered
simultaneously.

###  3.2 The molecular structure of clouds

Figure 7 shows the column densities of total hydrogen nuclei and of various
commonly-observed molecular species, viewed in the $x-y$ plane (i.e.
perpendicular to the collision axis). Except in the lowest density regions,
the H2 column density generally traces the total gas density well; hydrogen is
mostly in molecular form. CO, by contrast, is noticeably deficient in the low-
density parts of the cloud (the so-called ‘CO dark’ gas; Wolfire et al. 2010;
Glover & Mac Low 2011). Of the other molecules we consider, a distinction can
be drawn between those which essentially follow the distribution of CO (HCN,
CS, HNC, CN, C2H), and those which are biased towards higher density regions
(NH3, N2H+, HCO+, CH3OH), in that there is a higher column density contrast
when compared to that of the total gas column density. This is particularly
apparent for N2H+ (Figure 7), one of only a few species which is found to
genuinely trace dense gas preferentially (Pety et al., 2017; Kauffmann et al.,
2017; Tafalla et al., 2021). The latter category of molecules either include
nitrogen (NH3, N2H+), are charged (N2H+, HCO+), or are formed by a relatively
complex series of grain-surface reactions and subsequently returned to the gas
phase (CH3OH).

More quantitatively, Figure 8 shows the column density of CO versus the total
column density (note that this is the column along the line of sight, which is
distinct from the effective shielding column $N_{\rm H}^{\rm shield}$
discussed in Section 2.3; Clark & Glover 2014). The CO column remains fairly
low below $N_{\rm H}\sim 10^{21}\,{\rm cm}^{-2}$, but rises quite rapidly
beyond this point until seemingly reaching a constant value of $\sim
10^{18}\,{\rm cm}^{-2}$ for the highest density sightlines. The gas-phase CO
column never reaches the value expected if all elemental carbon was in this
form ($1.4\times 10^{-4}\,N_{\rm H}$), suggesting that photodissociation,
freeze-out, or some combination of both are always in effect to some extent.
The column densities of most other species, shown in Figure 9, are
qualitatively similar to CO, but the value of $N_{\rm H}$ at which the rapid
increase in molecular column begins varies significantly; N2H+, for example,
has a negligible column density for all sightlines with $N_{\rm
H}<10^{22}\,{\rm cm}^{-2}$, consistent with the observed threshold column
density at which this molecule becomes detectable in nearby666Barnes et al.
(2020) find different behaviour in the distant high-mass star forming region
W49B, with N2H+ becoming detectable at comparable column densities to HCN
($\sim 10^{21}\,{\rm cm}^{-2}$). While this could indicate chemical
differences compared to local Gould Belt clouds, a full investigation of the
underlying reasons is beyond the scope of this work. clouds (Pety et al.,
2017; Kauffmann et al., 2017; Tafalla et al., 2021).

The onset of this rise in molecular column effectively replicates the
distinction between CO-like species and those tracing higher densities
outlined above, with the latter category only existing in detectable
quantities (very roughly, molecular columns greater than $\sim 10^{10}\,{\rm
cm}^{-2}$; see e.g. Ragan et al. 2006) at higher values of $N_{\rm H}$. The
exception is NH3, which has a fairly smooth increase in molecular column
density over the entire range covered. In this case, the bias towards denser
gas apparent in Figure 7 is due to the fact that the NH3 column density
continues to increase with total column up to $10^{23}\,{\rm cm}^{-2}$,
whereas CO-like species plateau before this point due to depletion onto grain
surfaces, reducing the contrast between high- and low-density sightlines.

Figure 10 shows the distribution of the CO abundance at fixed gas density.
Scatter in the abundance reaches a minimum at a density of a few $10^{3}\,{\rm
cm}^{-3}$; in this regime, variations in other gas properties are not large
enough to significantly affect the chemistry, while freeze-out is still too
inefficient for the tracer history to have an impact. Even then, there is
still $\sim 0.5$ dex of scatter. While the median CO abundance declines quite
sharply above densities of $10^{4}\,{\rm cm}^{-3}$, the amount of scatter
increases substantially, with the interquartile range of our sample of tracer
particles spanning an order of magnitude or more.

The variation of molecular abundances with volume density for molecules other
than CO are shown in Figure 11; the median abundances and interquartile ranges
at various densities are given for all species in Table 2. As with the column
densities, there are two distinct classes of molecules: CO and molecules which
behave similarly (HCN, CS, HNC, CN), peaking in abundance at $n_{\rm H}\sim
10^{3}\,{\rm cm}^{-3}$ and then declining quite rapidly due to freeze-out; and
those which peak at $n_{\rm H}\sim 10^{4}\,{\rm cm}^{-3}$ (NH3, N2H+, HCO+,
CH3OH) and either decline more slowly toward higher densities, or maintain a
roughly constant abundance. The latter category also show significantly less
scatter in abundance at fixed density, typically less than $0.5$ dex, whereas
the abundances of CO-like molecules can vary by several orders of magnitude
beyond $10^{4}\,{\rm cm}^{-3}$. Of note is C2H, which shows unique behaviour
among the molecules investigated here. Its abundance peaks at $10^{3}\,{\rm
cm}^{-3}$, like CO, but then declines sharply and levels off to a nearly-
constant value above $10^{4}\,{\rm cm}^{-3}$, more comparable to nitrogen-
bearing species such as NH3 and N2H+. This is likely due to its sensitivity to
photodissociation, as discussed by Tafalla et al. (2021).

The abundances for HCN, N2H+, HCO+ and CS assumed by Smith et al. (2012, 2013)
are overplotted on the relevant panels in Figure 11. For all four species,
molecular abundances are seriously overestimated at low densities, and the
density-independent abundances assumed for HCN and HCO+ are much higher than
the NEATH values at high densities. The exponential depletion factor assumed
for CS, by contrast, significantly overestimates the actual extent of gas-
phase depletion for this molecule. Smith et al. (2012, 2013) used these
abundances to model low-$J$ line emission, so the poor agreement at high
densities may be of limited concern; these lines are typically optically
thick, and so the emission is dominated by relatively low-density ($\sim
10^{3}\,{\rm cm}^{-3}$) gas (Jones et al., 2023). However, the peak abundances
in our simulations are higher than the adopted values, and the constant N2H+
and HCO+ abundances are many orders of magnitude higher than the true values
at $10^{3}\,{\rm cm}^{-3}$, where significant line emission is expected to be
generated (Shirley, 2015). These differences are likely to dramatically alter
the resulting line emission properties.

The trends in median abundance versus density can be more accurately
reproduced by fifth-order polynomial fits, also overplotted on Figures 10 and
11, the coefficients for which are given in Appendix A. However, these
relationships are still insufficient to fully capture the behaviour of the
models. At high densities, abundances may differ by orders of magnitude from
the median values, and we emphasise again that this scatter is not random, but
is related to each particle’s individual evolutionary history (Figure 6).
Extracting consistent observational predictions from models requires a full
treatment of their time-dependent chemical evolution, something that we plan
to investigate in future work.

###  3.3 Stochasticity in the chemical composition of dense gas

Figure 12 shows the histories of five tracer particles with nearly identical
final densities of $\sim 10^{5}\,{\rm cm}^{-3}$. Despite the similar end
point, the density evolution varies dramatically between tracers. Some
particles undergo a relatively gradual rise from the initial to final density
over the entire $5\,{\rm Myr}$ of the simulation, while others remain at low
density until the last $\,{\rm Myr}$. All particles undergo at least one
period of shock compression. Some are shocked multiple times, and the density
evolution is not monotonic. Tracers can oscillate between $10^{2}-10^{3}\,{\rm
cm}^{-3}$ for long periods of time, although those which exceed the upper end
of this range seem to undergo irreversible collapse or compression. The five
tracers in Figure 12 were selected randomly, rather than being chosen as
extreme examples; this complex evolution is an inevitable consequence of the
dynamics of star-forming clouds.

The highly variable tracer histories result in similarly large variations in
their chemical evolution. The formation of significant quantities of CO seems
to be tied to the first instance of shock compression causing the gas density
to rise above $100\,{\rm cm}^{-3}$. The abundance then roughly follows the
evolution of the gas density, with rarefaction episodes tied to decreasing CO
abundance and vice versa, up until the density approaches $10^{5}\,{\rm
cm}^{-3}$ and freeze-out results in a rapid decline for all five tracers. The
onset and extent of this decline vary significantly, however, leading to a
difference of nearly two orders of magnitude in final CO abundance despite the
near-identical final densities. N2H+ behaves somewhat differently, forming
relatively late for all tracers, and roughly at the point where the density
exceeds $10^{3}\,{\rm cm}^{-3}$ for the first time. Compared to CO, the
abundance evolution is closer to a monotonic increase, but there are still
sharp, temporary declines (related to the CO abundance reaching a maximum, as
N2H+ is destroyed by reactions with CO), and for some tracers evidence of
late-time freeze-out.

Figure 12 also shows the evolution of the tracer abundances under a static
model, where the physical conditions are held constant at their final values.
This completely fails to reproduce the actual behaviour of the time-dependent
chemistry. Both CO and N2H+ effectively reach chemical equilibrium within
$\sim 1\,{\rm Myr}$, at abundances noticeably lower than their true final
values, and there are almost no chemical differences between the different
tracers at any point. The complete lack of chemical evolution over most of the
$5\,{\rm Myr}$ duration of the static model, compared to the highly complex
variations seen when physical evolution is taken into account, demonstrates
that any concept of chemical age based on static models will have little, if
any, correspondance to the true physical age of an object.

This variation in history affects the chemical composition of the gas in the
immediate surroundings of the sink particles, corresponding to the material in
prestellar or protostellar cores. Figure 13 shows the average abundances of
CO, HCN and N2H+ versus average gas density for the tracers within $0.05\,{\rm
pc}$ of each of the $28$ sink particles that have formed by the end of the
simulation.777Strictly speaking, these averages should be mass-weighted, but
as we do not have any information on the gas mass represented by each tracer,
we perform an unweighted average over all tracers. Even for cores with the
same (average) density, there are order-of-magnitude differences in the amount
of CO and HCN present in the gas-phase. Although N2H+ is more consistent
between sinks, it is clear that the chemical composition of the material being
delivered to pre- or protostellar cores can vary significantly from object to
object.

Table 2: Median molecular abundances and interquartile ranges at different gas densities. |  $n_{\rm H}$ / $\,{\rm cm}^{-3}$
---|---
| $10^{3}$ | $10^{4}$ | $10^{5}$ | $10^{6}$ | $10^{7}$
Molecule | Median
CO | $2.9\times 10^{-5}$ | $7.0\times 10^{-5}$ | $2.4\times 10^{-6}$ | $1.0\times 10^{-7}$ | $5.0\times 10^{-8}$
NH3 | $3.0\times 10^{-9}$ | $2.9\times 10^{-8}$ | $1.1\times 10^{-7}$ | $3.1\times 10^{-8}$ | $1.6\times 10^{-8}$
HCN | $1.5\times 10^{-9}$ | $9.0\times 10^{-9}$ | $4.9\times 10^{-9}$ | $2.2\times 10^{-10}$ | $6.7\times 10^{-11}$
N2H+ | $7.5\times 10^{-14}$ | $1.1\times 10^{-10}$ | $2.6\times 10^{-10}$ | $7.3\times 10^{-11}$ | $7.0\times 10^{-11}$
HCO+ | $1.5\times 10^{-11}$ | $3.2\times 10^{-9}$ | $4.4\times 10^{-9}$ | $8.6\times 10^{-10}$ | $5.9\times 10^{-10}$
CS | $4.4\times 10^{-9}$ | $1.1\times 10^{-8}$ | $8.2\times 10^{-10}$ | $4.2\times 10^{-11}$ | $1.2\times 10^{-11}$
HNC | $1.9\times 10^{-10}$ | $2.5\times 10^{-9}$ | $3.0\times 10^{-9}$ | $2.0\times 10^{-10}$ | $6.3\times 10^{-11}$
CH3OH | $2.7\times 10^{-11}$ | $3.0\times 10^{-9}$ | $4.6\times 10^{-10}$ | $8.6\times 10^{-11}$ | $4.7\times 10^{-11}$
CN | $9.6\times 10^{-8}$ | $3.5\times 10^{-8}$ | $8.6\times 10^{-9}$ | $6.9\times 10^{-10}$ | $2.3\times 10^{-10}$
C2H | $1.3\times 10^{-8}$ | $9.6\times 10^{-10}$ | $1.1\times 10^{-9}$ | $3.2\times 10^{-10}$ | $1.1\times 10^{-10}$
Molecule | Range [$x\times 10^{y}$ is represented as $x(y)$]
CO |  $5.1(-6)$ \- $7.7(-5)$ |  $4.0(-5)$ \- $9.4(-5)$ |  $5.4(-7)$ \- $1.1(-5)$ |  $4.4(-8)$ \- $5.2(-7)$ |  $2.0(-8)$ \- $1.7(-7)$
NH3 |  $2.1(-9)$ \- $4.6(-9)$ |  $1.5(-8)$ \- $6.1(-8)$ |  $7.7(-8)$ \- $1.3(-7)$ |  $1.3(-8)$ \- $6.0(-8)$ |  $5.2(-9)$ \- $4.4(-8)$
HCN |  $7.4(-10)$ \- $3.1(-9)$ |  $6.7(-9)$ \- $1.2(-8)$ |  $1.8(-9)$ \- $8.5(-9)$ |  $5.0(-11)$ \- $1.7(-9)$ |  $8.8(-12)$ \- $4.3(-10)$
N2H+ |  $3.4(-14)$ \- $2.4(-13)$ |  $9.9(-12)$ \- $3.9(-10)$ |  $1.0(-10)$ \- $4.4(-10)$ |  $4.9(-11)$ \- $1.6(-10)$ |  $4.9(-11)$ \- $1.3(-10)$
HCO+ |  $1.0(-11)$ \- $3.5(-11)$ |  $5.4(-10)$ \- $7.8(-9)$ |  $2.2(-9)$ \- $6.9(-9)$ |  $5.9(-10)$ \- $2.0(-9)$ |  $4.2(-10)$ \- $1.1(-9)$
CS |  $1.8(-9)$ \- $8.6(-9)$ |  $5.2(-9)$ \- $2.8(-8)$ |  $2.6(-10)$ \- $2.9(-9)$ |  $9.8(-12)$ \- $2.3(-10)$ |  $1.8(-12)$ \- $8.6(-11)$
HNC |  $1.1(-10)$ \- $4.5(-10)$ |  $1.7(-9)$ \- $4.2(-9)$ |  $1.4(-9)$ \- $5.0(-9)$ |  $4.7(-11)$ \- $1.3(-9)$ |  $8.4(-12)$ \- $3.8(-10)$
CH3OH |  $3.4(-12)$ \- $1.4(-10)$ |  $1.6(-9)$ \- $4.3(-9)$ |  $2.5(-10)$ \- $9.2(-10)$ |  $4.2(-11)$ \- $2.0(-10)$ |  $2.0(-11)$ \- $1.4(-10)$
CN |  $5.0(-8)$ \- $1.3(-7)$ |  $2.3(-8)$ \- $5.1(-8)$ |  $3.8(-9)$ \- $1.6(-8)$ |  $1.8(-10)$ \- $3.9(-9)$ |  $3.3(-11)$ \- $1.3(-9)$
C2H |  $7.6(-9)$ \- $1.8(-8)$ |  $4.4(-10)$ \- $2.3(-9)$ |  $5.7(-10)$ \- $1.7(-9)$ |  $8.9(-11)$ \- $8.5(-10)$ |  $1.6(-11)$ \- $5.3(-10)$

##  4 Discussion

###  4.1 Chemical tracers in molecular clouds

Using MHD simulations with an essentially complete treatment of gas dynamics
and the relevant thermochemistry, we are able to investigate how molecular
abundances evolve under realistic conditions, as opposed to the idealised
models commonly adopted in previous studies. Of the species we investigate
here, there appears to be a distinction between CO-like molecules, which are
effectively depleted from the gas-phase beyond densities of $10^{4}\,{\rm
cm}^{-3}$, and those where some competing process mitigates freeze-out onto
grain surfaces. The latter (NH3-like) molecules may show this behaviour due to
reduced destruction by gas-phase CO (e.g. N2H+), their requirement to be
formed via depletion of said CO (CH3OH), or some other reason. In all cases
they present nearly constant abundances at high density, and relatively little
scatter compared to CO-like species. They are thus relatively insensitive to
the prior evolution of the gas, and may represent good targets for tracing the
present-day conditions of denser structures within molecular clouds. The
abundances of CO-like molecules, on the other hand, are strongly dependent on
the evolutionary history (Figure 6), and so would make better tracers of, for
example, core ages. However, the complex and non-monotonic chemical evolution
of all species (Figure 12) makes either of these uses challenging. Diagnostics
based on simplified evolutionary models should be considered speculative at
best.

###  4.2 Dynamic physical evolution in prestellar cores

Chemical models of prestellar cores often assume a static density profile, and
subsequently evolve the chemistry at different points along this profile for
several $\,{\rm Myr}$ in order to calculate molecular abundances (e.g. Caselli
et al., 2022; Megías et al., 2022). Figure 5 shows that this will inevitably
overestimate the amount of depletion in the high-density central regions of
cores, even for species which are supposedly resistant to freeze-out such as
N2H+, because it is unrealistic for gas to remain at these densities for this
length of time (orders of magnitude longer than the free-fall timescale). The
material at high density will likely have spent only a very small fraction of
time at this density, and so will retain much of its gas-phase molecular
content. The chemical makeup also depends critically on the core’s prior
evolutionary history; cores with identical density structures can (and most
likely will) have completely different molecular compositions, which depend on
the details of their formation. While we have focused on the gas-phase
abundances, Clément et al. (2023) find similar history-dependent effects for
the composition of ices in prestellar cores. These effects may have some
bearing on otherwise-unexpected phenomena such as the lack of depletion of
nitrogen-bearing molecules (Sipilä et al., 2019) and the ubiquitous detection
of complex organic molecules (Jiménez-Serra et al., 2016; Jiménez-Serra et
al., 2021; Scibelli & Shirley, 2020; Scibelli et al., 2021; Ambrose et al.,
2021) in cold prestellar environments.

###  4.3 Comparison with previous work

Most previous studies using large chemical networks on the scales of molecular
clouds adopt a simplified treatment of the thermal physics (e.g. Ruaud et al.,
2018; Wakelam et al., 2019, 2020; Bovino et al., 2021; Jensen et al., 2021;
Ferrada-Chamorro et al., 2021; Gómez et al., 2022; Priestley et al., 2023b).
To our knowledge, the only other peer-reviewed study on these scales888Ilee et
al. (2017) also post-process hydrodynamical models which include a detailed
treatment of the thermal evolution, but on the much smaller scales of
protoplanetary discs. Preliminary results from a comparable study to ours were
presented in a conference proceedings by Szucs et al. (2015). which uses a
hydrodynamical simulation with self-consistent thermodynamical evolution as a
base for chemical post-processing is that of Panessa et al. (2023), who
primarily focus on the abundance of HCO+. They find a similar peak HCO+ column
density to ours ($\sim 10^{15}\,{\rm cm}^{-2}$), and a similar threshold value
of $N_{\rm H}$ beyond which the abundance sharply rises ($\sim 10^{22}\,{\rm
cm}^{-2}$). While this agrees with our classification of HCO+ as an NH3-like
molecule, we note that this behaviour is in conflict with its observed line
emission properties (Tafalla et al., 2021), suggesting that some part of the
reaction network for this molecule is incorrect (Priestley et al., 2023b).

The key difference between our simulations and those of Panessa et al. (2023)
is that we use an idealised model of colliding spherical clouds, whereas they
study zoomed-in regions from within a galactic disc-scale simulation (Walch et
al., 2015; Seifried et al., 2017). Our approach allows more flexibility in
choosing cloud parameters, rather than being restricted to whatever happens to
form in the parent simulations; theirs has the benefit of forming clouds self-
consistently from the larger-scale structure, rather than beginning with
arbitrary (and possibly unrealistic) initial conditions. We note that
idealised cloud models of the sort considered here have been shown to
accurately reproduce the line emission properties of observed molecular clouds
(Priestley et al., 2023b), suggesting that they are successfully capturing the
important details of the chemical evolution.

Lupi et al. (2021) take the alternative approach to post-processing, and
expand the internal chemical network in their underlying MHD simulations to
include more species and reactions. While this allows for a completely self-
consistent treatment of the coupled physical, thermal, and chemical evolution,
the increased computational expense puts a limit on the size of network that
can be employed; Lupi et al. (2021) do not include any nitrogen chemistry, for
example, a key focus of this work due to the observational importance of many
nitrogen-bearing species. Given that the additional chemical complexity is
unlikely to greatly affect the gas or dust temperatures (Goldsmith, 2001; Gong
et al., 2017), post-processing remains a viable alternative for modelling
time-dependent chemistry in molecular clouds.

##  5 Conclusions

We present a framework for chemically post-processing MHD simulations of
molecular clouds. Our post-processed molecular abundances are consistent with
those of the underlying MHD simulation for the small number of species
required to follow the thermal evolution of the gas, and which are therefore
included in the on-the-fly chemical evolution. Our post-processing framework
allows us to investigate the behaviour of many additional species at low
computational cost, such as nitrogen-bearing molecules which are
observationally important as dense gas tracers.

We find that most molecules fall into one of two families: CO-like molecules,
which reach peak abundances at gas densities of $10^{3}\,{\rm cm}^{-3}$,
decline toward higher densities, and have large scatter in abundance at fixed
density due to the varying evolutionary history of the gas; and NH3-like
molecules, which reach their peak abundance at higher ($\sim 10^{4}\,{\rm
cm}^{-3}$) densities, decline either slowly or not at all beyond this peak,
and have a tighter (although still not one-to-one) relationship between
density and abundance.

Evolving the chemistry with the physical properties held constant produces
results accurate to within a factor of a few up to densities of $\sim
10^{3}\,{\rm cm}^{-3}$, but diverges from the true behaviour at higher
densities. Ignoring the dynamical evolution of the gas tends to severely
overestimate the amount of depletion of molecules onto grain surfaces.
Differences in the prior evolution of tracer particles additionally result in
significant scatter in the molecular abundances at fixed gas density, in
contrast to static-density models. Chemical models which assume a static
density structure will inevitably produce erroneous (and potentially
misleading) results; the objects being modelled are neither static nor in
equilibrium, and their simultaneous physical and chemical evolution must be
taken into account.

## Acknowledgements

FDP, PCC, SER and OF acknowledge the support of a consolidated grant
(ST/K00926/1) from the UK Science and Technology Facilities Council (STFC).
SCOG and RSK acknowledge funding from the European Research Council (ERC) via
the ERC Synergy Grant “ECOGAL-Understanding our Galactic ecosystem: From the
disk of the Milky Way to the formation sites of stars and planets” (project ID
855130, from the Heidelberg Cluster of Excellence (EXC 2181 - 390900948)
“STRUCTURES: A unifying approach to emergent phenomena in the physical world,
mathematics, and complex data”, funded by the German Excellence Strategy, and
from the German Ministry for Economic Affairs and Climate Action in project
“MAINN” (funding ID 50OO2206). The team in Heidelberg also thanks for
computing resources provided by The Länd through bwHPC and DFG through grant
INST 35/1134-1 FUGG and for data storage at SDS@hd through grant INST
35/1314-1 FUGG. This research was undertaken using the supercomputing
facilities at Cardiff University operated by Advanced Research Computing at
Cardiff (ARCCA) on behalf of the Cardiff Supercomputing Facility and the
Supercomputing Wales (SCW) project. We acknowledge the support of the latter,
which is part-funded by the European Regional Development Fund (ERDF) via the
Welsh Government.

## Data Availability

The data underlying this article will be shared on request.

## References

  * Aikawa et al. (2005) Aikawa Y., Herbst E., Roberts H., Caselli P., 2005, ApJ, 620, 330
  * Ambrose et al. (2021) Ambrose H. E., Shirley Y. L., Scibelli S., 2021, MNRAS, 501, 347
  * Banerji et al. (2009) Banerji M., Viti S., Williams D. A., Rawlings J. M. C., 2009, ApJ, 692, 283
  * Barnes et al. (2020) Barnes A. T., et al., 2020, MNRAS, 497, 1972
  * Bergin & Tafalla (2007) Bergin E. A., Tafalla M., 2007, ARA&A, 45, 339
  * Bohlin et al. (1978) Bohlin R. C., Savage B. D., Drake J. F., 1978, ApJ, 224, 132
  * Boogert et al. (2015) Boogert A. C. A., Gerakines P. A., Whittet D. C. B., 2015, ARA&A, 53, 541
  * Bovino et al. (2021) Bovino S., Lupi A., Giannetti A., Sabatini G., Schleicher D. R. G., Wyrowski F., Menten K. M., 2021, A&A, 654, A34
  * Caselli et al. (2022) Caselli P., et al., 2022, ApJ, 929, 13
  * Clark & Glover (2014) Clark P. C., Glover S. C. O., 2014, MNRAS, 444, 2396
  * Clark et al. (2012a) Clark P. C., Glover S. C. O., Klessen R. S., 2012a, MNRAS, 420, 745
  * Clark et al. (2012b) Clark P. C., Glover S. C. O., Klessen R. S., Bonnell I. A., 2012b, MNRAS, 424, 2599
  * Clark et al. (2019) Clark P. C., Glover S. C. O., Ragan S. E., Duarte-Cabral A., 2019, MNRAS, 486, 4622
  * Clément et al. (2023) Clément A., et al., 2023, arXiv e-prints, p. arXiv:2306.08346
  * Federman et al. (1979) Federman S. R., Glassgold A. E., Kwan J., 1979, ApJ, 227, 466
  * Ferrada-Chamorro et al. (2021) Ferrada-Chamorro S., Lupi A., Bovino S., 2021, MNRAS, 505, 3442
  * Genel et al. (2013) Genel S., Vogelsberger M., Nelson D., Sijacki D., Springel V., Hernquist L., 2013, MNRAS, 435, 1426
  * Glover & Clark (2012) Glover S. C. O., Clark P. C., 2012, MNRAS, 421, 116
  * Glover & Mac Low (2007) Glover S. C. O., Mac Low M.-M., 2007, ApJS, 169, 239
  * Glover & Mac Low (2011) Glover S. C. O., Mac Low M. M., 2011, MNRAS, 412, 337
  * Goldsmith (2001) Goldsmith P. F., 2001, ApJ, 557, 736
  * Gómez et al. (2022) Gómez G. C., Walsh C., Palau A., 2022, MNRAS, 513, 1244
  * Gong et al. (2017) Gong M., Ostriker E. C., Wolfire M. G., 2017, ApJ, 843, 38
  * Habing (1968) Habing H. J., 1968, Bull. Astron. Inst. Netherlands, 19, 421
  * Hily-Blant et al. (2022) Hily-Blant P., Pineau des Forêts G., Faure A., Lique F., 2022, A&A, 658, A168
  * Holdship & Viti (2022) Holdship J., Viti S., 2022, A&A, 658, A103
  * Holdship et al. (2017) Holdship J., Viti S., Jiménez-Serra I., Makrymallis A., Priestley F., 2017, AJ, 154, 38
  * Hunter et al. (2023) Hunter G. H., Clark P. C., Glover S. C. O., Klessen R. S., 2023, MNRAS, 519, 4152
  * Ilee et al. (2017) Ilee J. D., et al., 2017, MNRAS, 472, 189
  * Jenkins (2009) Jenkins E. B., 2009, ApJ, 700, 1299
  * Jensen et al. (2021) Jensen S. S., Jørgensen J. K., Furuya K., Haugbølle T., Aikawa Y., 2021, A&A, 649, A66
  * Jensen et al. (2023) Jensen S. S., Spezzano S., Caselli P., Grassi T., Haugbølle T., 2023, A&A, 675, A34
  * Jiménez-Serra et al. (2016) Jiménez-Serra I., et al., 2016, ApJ, 830, L6
  * Jiménez-Serra et al. (2021) Jiménez-Serra I., Vasyunin A. I., Spezzano S., Caselli P., Cosentino G., Viti S., 2021, ApJ, 917, 44
  * Jones et al. (2023) Jones G. H., Clark P. C., Glover S. C. O., Hacar A., 2023, MNRAS, 520, 1005
  * Kaminski et al. (2014) Kaminski E., Frank A., Carroll J., Myers P., 2014, ApJ, 790, 70
  * Kauffmann et al. (2017) Kauffmann J., Goldsmith P. F., Melnick G., Tolls V., Guzman A., Menten K. M., 2017, A&A, 605, L5
  * Larson (1969) Larson R. B., 1969, MNRAS, 145, 271
  * Lupi et al. (2021) Lupi A., Bovino S., Grassi T., 2021, A&A, 654, L6
  * Mathis et al. (1983) Mathis J. S., Mezger P. G., Panagia N., 1983, A&A, 128, 212
  * McElroy et al. (2013) McElroy D., Walsh C., Markwick A. J., Cordiner M. A., Smith K., Millar T. J., 2013, A&A, 550, A36
  * Megías et al. (2022) Megías A., Jiménez-Serra I., Martín-Pintado J., Vasyunin A. I., Spezzano S., Caselli P., Cosentino G., Viti S., 2022, MNRAS,
  * Navarro-Almaida et al. (2020) Navarro-Almaida D., et al., 2020, A&A, 637, A39
  * Panessa et al. (2023) Panessa M., Seifried D., Walch S., Gaches B., Barnes A. T., Bigiel F., Neumann L., 2023, MNRAS, 523, 6138
  * Pety et al. (2017) Pety J., et al., 2017, A&A, 599, A98
  * Priestley et al. (2023a) Priestley F. D., Whitworth A. P., Fogerty E., 2023a, MNRAS, 518, 4839
  * Priestley et al. (2023b) Priestley F. D., Clark P. C., Whitworth A. P., 2023b, MNRAS, 519, 6392
  * Prole et al. (2022) Prole L. R., Clark P. C., Klessen R. S., Glover S. C. O., 2022, MNRAS, 510, 4019
  * Ragan et al. (2006) Ragan S. E., Bergin E. A., Plume R., Gibson D. L., Wilner D. J., O’Brien S., Hails E., 2006, ApJS, 166, 567
  * Rocha et al. (2023) Rocha C. M. R., et al., 2023, arXiv e-prints, p. arXiv:2307.00311
  * Ruaud et al. (2018) Ruaud M., Wakelam V., Gratier P., Bonnell I. A., 2018, A&A, 611, A96
  * Scibelli & Shirley (2020) Scibelli S., Shirley Y., 2020, ApJ, 891, 73
  * Scibelli et al. (2021) Scibelli S., Shirley Y., Vasyunin A., Launhardt R., 2021, MNRAS, 504, 5754
  * Seifried et al. (2017) Seifried D., et al., 2017, MNRAS, 472, 4797
  * Sembach et al. (2000) Sembach K. R., Howk J. C., Ryans R. S. I., Keenan F. P., 2000, ApJ, 528, 310
  * Shirley (2015) Shirley Y. L., 2015, PASP, 127, 299
  * Sipilä & Caselli (2018) Sipilä O., Caselli P., 2018, A&A, 615, A15
  * Sipilä et al. (2019) Sipilä O., Caselli P., Redaelli E., Juvela M., Bizzocchi L., 2019, MNRAS, 487, 1269
  * Smith et al. (2012) Smith R. J., Shetty R., Stutz A. M., Klessen R. S., 2012, ApJ, 750, 64
  * Smith et al. (2013) Smith R. J., Shetty R., Beuther H., Klessen R. S., Bonnell I. A., 2013, ApJ, 771, 24
  * Smith et al. (2020) Smith R. J., et al., 2020, MNRAS, 492, 1594
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Szucs et al. (2015) Szucs L., Glover S., Caselli P., 2015, in EAS Publications Series. pp 391–392, doi:10.1051/eas/1575079
  * Tafalla et al. (2002) Tafalla M., Myers P. C., Caselli P., Walmsley C. M., Comito C., 2002, ApJ, 569, 815
  * Tafalla et al. (2021) Tafalla M., Usero A., Hacar A., 2021, A&A, 646, A97
  * Tress et al. (2020) Tress R. G., Smith R. J., Sormani M. C., Glover S. C. O., Klessen R. S., Mac Low M.-M., Clark P. C., 2020, MNRAS, 492, 2973
  * Wakelam et al. (2019) Wakelam V., Ruaud M., Gratier P., Bonnell I. A., 2019, MNRAS, 486, 4198
  * Wakelam et al. (2020) Wakelam V., Iqbal W., Melisse J. P., Gratier P., Ruaud M., Bonnell I., 2020, MNRAS, 497, 2309
  * Walch et al. (2015) Walch S., et al., 2015, MNRAS, 454, 238
  * Wolfire et al. (2010) Wolfire M. G., Hollenbach D., McKee C. F., 2010, ApJ, 716, 1191
  * Wurster (2021) Wurster J., 2021, MNRAS, 501, 5873
  * de Jong (1977) de Jong T., 1977, A&A, 55, 137
  * van Dishoeck & Black (1988) van Dishoeck E. F., Black J. H., 1988, ApJ, 334, 771

##  Appendix A Analytical fits to molecular abundances

Table 3 lists the coefficients of the polynomial fits to median molecular
abundance shown in Figures 10 and 11, in logarithmic space so that

$\log_{10}(n_{X}/n_{\rm
H})=\displaystyle\sum_{i=0}^{5}a_{i}\left(\log_{10}n_{\rm H}\right)^{i}$ (7)

with $n_{X}$ and $n_{\rm H}$ in units of $\,{\rm cm}^{-3}$. While these values
are reasonably accurate in reproducing the median abundance at a given
density, there is no one-to-one connection between density and abundance for
any species, with variations particularly significant for CO-like molecules at
densities above $10^{3}\,{\rm cm}^{-3}$. In addition, it remains to be seen
how well these fits will reproduce the median abundances in molecular clouds
formed in conditions very different to those simulated here, e.g. in the Milky
Way’s Central Molecular Zone. As such, the polynomial fits should be seen as
rough approximations, suitable when order-of-magnitude level errors can be
tolerated, and not as a genuine alternative to full chemical modelling.

Table 3: Coefficients of polynomial fits to median molecular abundances as a function of density. Molecule | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $a_{5}$
---|---|---|---|---|---|---
CO | $-0.044$ | $1$ | $-8.7$ | $34$ | $-60$ | $30$
NH3 | $0.0073$ | $-0.12$ | $0.55$ | $-0.55$ | $-0.92$ | $-8.2$
HCN | $-0.0041$ | $0.12$ | $-1.2$ | $4.8$ | $-6$ | $-10$
N2H+ | $-0.018$ | $0.52$ | $-5.7$ | $28$ | $-62$ | $34$
HCO+ | $-0.022$ | $0.62$ | $-6.6$ | $33$ | $-74$ | $51$
CS | $-0.022$ | $0.5$ | $-4.2$ | $15$ | $-22$ | $-2.2$
HNC | $-0.0011$ | $0.077$ | $-1.1$ | $5.9$ | $-12$ | $-3.1$
CH3OH | $-0.065$ | $1.5$ | $-13$ | $54$ | $-100$ | $55$
CN | $0.0064$ | $-0.17$ | $1.9$ | $-11$ | $29$ | $-38$
C2H | $0.036$ | $-0.89$ | $8.6$ | $-40$ | $90$ | $-86$