Quasi-equilibrium chemical evolution in starless cores
Jonathan Rawlings Eric Keto Paola Caselli


firstpage–lastpage
Revisiting Out-of-Distribution Detection in LiDAR-based 3D Object Detection

    Michael Kösel^1, Marcel Schreiber^2, Michael Ulrich^2, Claudius Gläser^2 and Klaus Dietmayer^1
	^1Institute of Measurement, Control, and Microtechnology, Ulm University, Germany, {first.last}@uni-ulm.de
    ^2Robert Bosch GmbH, Corporate Research, 71272 Renningen, Germany, {first.last}@de.bosch.com and michael.ulrich2@bosch.com

    April 29, 2024
==============================================================================================================================================================================================================================================================================================================================================




The chemistry of H_2O, CO and other small molecular species in 
an isolated pre-stellar core, L1544, has been assessed in 
the context of a comprehensive gas-grain chemical model, coupled to 
an empirically constrained physical/dynamical model.

Our main findings are (i) 
that the chemical network remains in near equilibrium as the 
core evolves towards star formation and the molecular abundances change 
in response to the evolving physical conditions. 
The gas-phase abundances at any time can be calculated accurately
with equilibrium chemistry, and the concept of chemical clocks is
meaningless in molecular clouds with similar conditions and dynamical 
time scales, and  
(ii) A comparison of the results of complex and simple chemical networks
indicates that the abundances of the dominant oxygen and carbon species, 
H_2O, CO, C, and C^+ are reasonably approximated by simple networks. 
In chemical equilibrium, the time-dependent differential terms vanish and 
a simple network reduces to a few algebraic equations. This allows rapid 
calculation of the abundances most responsible for spectral line radiative 
cooling in molecular clouds with long dynamical time scales. 

The dust ice mantles are highly structured and the ice 
layers retain a memory of the gas-phase abundances at the time of their
deposition.
A complex (gas-phase and gas-grain) chemical structure therefore exists, 
with cosmic-ray induced processes dominating in the inner regions. 

The inferred H_2O abundance profiles for L1544 require that the outer parts 
of the core and also any medium exterior to the core are 
essentially transparent to the interstellar radiation field.



astrochemistry  –  molecular processes  –  
ISM: clouds  –  ISM: molecules  –  stars:formation
 –  stars:protostars








§ INTRODUCTION



In this paper, we seek to improve our modelling and understanding of the chemistry 
of quiescent molecular gas with our astrochemistry code, STARCHEM
<cit.>. 
The STARCHEM code is designed to integrate chemical evolution along Lagrangian streamlines of fluid flow. In this paper, we apply STARCHEM in conjunction with 
a previously developed dynamical model for the gravitationally contracting molecular 
cloud L1544 <cit.>. 
Our goals in this study are twofold.

Firstly, we test the hypothesis that the chemical evolution can be treated as being
in quasi-equilibrium, so that the gas-phase abundance distributions are simply 
determined by the instantaneous physical conditions. 
As part of this analysis, we pay particular attention to the chemical stratification 
within the dust grain ice mantles, and so differentiate between the gas-phase 
chemistry and that of the ices, which record the chemical evolution of the core.

Secondly, we assess the accuracy and applicability of a radically simplified 
chemical model developed in a series of earlier papers; KC10 and <cit.> 
(hereafter KC08 and KRC14, respectively). 
This model is one of several proposals in the literature for reduced chemical 
networks that are also useful in rapid computation of molecular gas temperatures 
<cit.>.
They directly compute the abundances of the important coolants such as CO, C^+, 
C and O. Abundances of other molecular species can be estimated from those 
calculated to develop a general cooling function <cit.>. 
Because of their computational speed, these reduced 
chemical networks are useful in numerical, hydrodynamic simulations of the 
ISM which are already computationally intensive without the additional burden 
of a time-dependent complex chemistry with many reactions. An earlier 
version of our reduced chemical network has been incorporated into one 
such code <cit.>.

We choose L1544 as an example because we have an observationally verified model 
for the dynamics of this cloud that provides specific combinations of gas and 
dust temperatures,  gas density, and UV radiation for our simulations of the astrochemistry. In addition, the different radial zones in the cloud, outer, 
central, and transition, provide examples where the chemistry operates in very 
different ways. Whilst L1544 is a particular cloud, its chemistry is 
representative of quiescent molecular gas in the ISM meaning gas that is not 
strongly disturbed by shocks, bipolar outflows, intense radiation from 
nearby stars, or processes in accretion disks.

In the following section we discuss the structure and dynamics of the L1544 
pre-stellar core. Our dynamical/chemical model is described in section <ref>. 
In section <ref> we present the results discussing (a) the chemical 
structure of the  cloud (b) the structure of the ice mantles (c) the chemical 
time-dependence, and (d) the validity of reduced chemical networks. 
Our conclusions are summarized in section section <ref>. A proposed 
revised form of our reduced chemical network is given in 
Appendix <ref> and the main results of STARCHEM applicable to 
starless and pre-stellar cores are described in Appendix <ref>.



§ THE PHYSICAL STRUCTURE AND DYNAMICS OF L1544



The pre-stellar core L1544 is located in the Taurus Molecular Cloud, at a 
distance of ∼170pc <cit.>. Infrared observations show no evidence 
of an embedded protostar; however, split, self-absorbed spectral line profiles 
of several molecules that are tracers of dense gas indicate that the central 
region is in the early stages of gravitational collapse toward star formation. 
On a larger scale, molecular line observations show patchy emission around the 
core suggesting that the core is isolated within the Taurus cloud with no 
evidence for surrounding filamentary structure or other cores overlapping 
along the line of sight. The overall morphology is roughly elliptical with 
an aspect ratio of 2:1; however contours of constant density become more 
spherical with increasing density suggesting that the dynamics can be 
understood in the context of a spherical model <cit.>. 
With these properties, L1544 has become a prototypical example of a small 
molecular cloud that is just beginning gravitational collapse towards the 
formation of a solar-mass scale star. 

The chemistry within L1544 is both relatively simple and interestingly 
complex. The absence of a protostar implies the absence of the chemistry 
associated with shocks, a bipolar outflow, and an accretion disk. The chemistry 
is consistent with that of quiescent molecular gas. With no embedded protostar 
or nearby bright stars, the core is illuminated only by the external radiation 
field from Galactic starlight. This sets up a radial structure within the 
core of varied chemistry. 
Photochemistry dominates the outer region, but the inner core is cold 
<cit.> and dark, so that the chemistry is mostly 
restricted to grain surfaces and the 
gas-grain interactions of freeze-out and desorption. L1544 is in a family of 
starless cores in which CO has been highly depleted <cit.>.
Indeed, recent high resolution observations have shown almost completele 
depletion in the central regions of the core 
<cit.>.

Along with its advantages for theoretical study, L1544 is observationally 
favorable. At 170 pc distant, the core is extended over arc minute scales 
and well matched to the angular resolution of single dish radio telescopes 
at millimeter and submillimeter frequencies. A number of spectral lines has 
been observed with varying optical depths from molecules with different 
critical densities for collisional de-excitation which provide information 
on both the inner and outer regions of the core. While the 
dynamics and the chemistry can be studied independently, it is not possible 
to understand the observations without consideration of both. The most 
obvious effect is that the molecular abundances and gas velocity are both 
functions of radius resulting in differences in the self-absorbed line 
profiles of different molecules. Once understood, this effect provides precise 
estimates of the emission and absorption regions for each molecule. 

In a series of earlier papers we developed a model for the chemistry and 
dynamics of L1544 that has been successful in both matching and predicting 
the spectral line shapes of several molecules. Originally based on 
observations of CO and N_2H^+, the model was able without modification to 
predict the blue-shifted self-absorbed emission of H_2O and match the 
complex pattern of the self-absorbed hyperfine lines of NH_3 
<cit.>.

The physical model of L1544 is described in KC10, as resulting from 
the quasi-static contraction of a thermally supercritical 10M_⊙ cloud, 
of radius 0.75pc that is initially configured as a Bonnor-Ebert sphere, in 
unstable dynamical equilibrium, with a central density of 
2× 10^4cm^-3. 
After contraction to a central density of 2× 10^7cm^-3, the 
structure of the model describes the observed structure of L1544. Subsequent 
observations suggest minor differences <cit.>.
<cit.> proposed that the `flat central region' is more extended and 
has a lower density (∼ 2× 10^6cm^-3) than that of KC10 - based 
on 1.1mm Large Millimetre Telescope and 3.3mm Green Bank Observatory data.
Based on 1.3mm ALMA data, <cit.> found that it is somewhat smaller, 
elongated, but consistent with a central density of ∼ 1× 10^7cm^-3.
Observationally, the density in the center of the core cannot be determined
directly, only the column density through the centre. The suggested 
differences are consistent with observational uncertainty.






Figure <ref> shows the temperature (gas and dust), density, 
velocity and extinction profiles; T_ gas(r), T_ dust(r), 
n(r), v(r) and A_ v(r), at several times between 0.30Myr 
(t_ min) and 0.61Myr (t_ max) since the start of the 
contraction from the initial configuration. The profiles for 
t=t_ max (which are labelled and indicate the direction of time) 
correspond to those published as the best fits in 
<cit.> (hereafter KRC14) and subsequently used in other studies.

The A_ v(r,t) profiles give the, astrochemically 
relevant, mean direction-averaged extinction at any particular point 
and time within the spherical model of the core.
[Since starless cores have no internal sources of 
radiation, the radiation field at any point in the core is the radiation 
field exterior to the core diminished by the mean Av averaged over all 
directions to the boundary of the core. The mathematical definition of 
the average can be found in Appendix 1 of <cit.>.]
This is not the same as the, more observationally relevant, line of
sight extinction - which is what is depicted in, for example, figure 3 of KRC14.

From Figure <ref> we can see that the density profiles progress 
towards n∝ r^-2
as the central region of constant density defined by the radius where 
the crossing time equals the free-fall time decreases in size.
For L1544, the profiles for t=t_ max (and, in particular, the density
and velocity profiles, and the value of the central density 
∼ 10^7cm^-3) are empirically constrained by the observed 
H_2O 576GHz (1_10-1_11) line profile. This density corresponds to a 
specific time (t=t_ max∼0.61Myr) in the dynamical evolution, 
as can be seen from the time-dependence of the velocity profiles. At this time 
the dynamical evolution of the cloud is beginning to accelerate, with 
infall velocities just becoming supersonic in the inner region 
(r0.01pc, A_ v 23). 

It is also evident that the gas and dust temperatures are reasonably 
well-coupled within ∼0.15pc (A_ v 0.96), and are 
very low throughout most of the cloud. 
This has significant implications for the chemistry which we discuss in
section <ref>.






Although the core evolves very slowly, the physical model nonetheless 
implies an important time-dependence, so that the best fit to the 
observations and the implied conditions representing the present status of 
a core are defined by a specific dynamical age. 






In general, there are essentially two possibilities: 
(i) the timescale for chemical evolution (τ_ chem) is less than 
that for the dynamical evolution (τ_ dyn), 
in which case quasi-equilibrium chemical evolution will hold - with the 
abundances responding to the changes in the instantaneous physical conditions, or 
(ii) τ_ chem>τ_ dyn, in which case the chemistry will 
depend on both the dynamical evolution and the initial conditions.
In either case, a static chemistry is not a realistic representation.





The fastest timescale in which 
a Bonnor-Ebert configuration in hydrostatic equilibrium can be assembled 
via sub-sonic quasi-static contraction is of the order of the free-fall 
collapse timescale:

    t_ ff = √(3π/32Gρ) = 3.1× 
    10^5(10^4 cm^-3/n_ H_2)^1/2  years

where n_ H_2 is the (assumed molecular) gas density and we have assumed 
that  He/H =0.1.

Furthermore, as the core must be in near-hydrostatic equilibrium
(KCR15) its age must be greater than the sound-crossing time 
(t_ cross). 
Even if we only consider the inner 0.1pc then, for a core composed of 
molecular gas, with mean temperature ∼10K; t_ cross>1Myr. 

Shorter timescales than these could only be obtained as a result of 
supra-gravitational, supersonic contraction, 
which would occur only if the core were out of equilibrium everywhere. 
As shown in KCR15 an increase of the density by ∼ 10% over its 
equilibrium values along the profile would result in global cloud collapse 
and a velocity structure that is incompatible with the subsonic spectral 
line profiles of all observed molecules other than H_2O.

Moreover, as can be inferred from Fig.1, if the dynamical age were only 
slightly larger than this, then much higher central densities, together with 
larger inflow velocities throughout the envelope would result.

It has taken t_ form∼0.6Myr to 
construct the existing configuration from an initial state in dynamical 
equilibrium that was marginally unstable 
and (from the implied velocity profile) that it is now evolving on a near 
gravitational free-fall time-scale of <0.3Myr.
Thus t_ form is a fixed parameter. Any model of the 
chemistry that identifies a nominal `chemical age' that is independent of 
the physical evolution would be quite inconsistent with this, or any 
other realistic, dynamical model. The inapplicability of hydrodynamically
static chemical evolution was also noted by <cit.> in their 
combined chemical-hydrodynamical study of the velocity structure in L1544. 

We should also consider the physical and chemical 
evolution of the cloud prior to the formation of the initial 
Bonnor-Ebert sphere, since parts of this cloud will already have relatively 
high extinctions and low temperatures, implying that they must have previously 
evolved through a PDR interface.

Moreover, in the circumstances where the chemical evolution is slow 
(τ_ chem>τ_ dyn), then a larger effective chemical age could
be obtained if it is assumed that the core maintains hydrostatic equilibrium 
for some time before quasi-static
contraction starts, so that the chemical evolution occurs at relatively 
low densities (with correspondingly longer chemical timescales).



§ THE CHEMICAL MODEL OF STARLESS CORES



To understand the chemical structure and evolution of isolated starless
cores we utilise
the STARCHEM dynamical/chemical model as described in Appendix B.


When modelling the astrochemistry of pre- and proto-stellar cores, 
the usual approach is to study the time-dependent
chemical evolution, but freezing the physical conditions at t=t_ max, 
with the inclusion of some chemical pre-evolution at constant density.
As is clear from the discussion in the previous section, this is quite arbitrary 
and incompatible with our understanding of how the source has evolved.






This has been noted in previous studies (KC10, KRC14) in the context 
of the time-dependences of the abundance profiles for CO and H_2O in L1544, 
this important fact has not been considered in most subsequent modelling efforts.

Moreover, even if we make the assumption that chemical quasi-equilibrium holds
(i.e. in response to the changing physical conditions), then the abundances
will be determined by the equilibrium balance between the gas-phase chemistry 
and the chemistry in the surface layers of the ice mantles on dust grains.
In conditions where significant freeze-out is occurring then the sub-surface 
layers of the ices - which are effectively chemically inert - will be determined 
by the entire evolutionary history of the cloud/core.  

We go some way towards addressing these inconsistencies and 
assumptions by following the combined dynamical and chemical evolution between 
t_ min and t_max, although we also make the simplification
that the core `instantly assembles' in the configuration specified at
t=t_ min.
Indeed, even if the initial configuration is taken to be a sphere 
of uniform density, that sphere will itself have position-dependent 
temperature and extinction profiles that will result in significant, 
and time-dependent, chemical structure. A fully self-consistent approach
would involve a knowledge of the density, temperature and extinction 
profiles from the very earliest stages.






STARCHEM considers a full multi-point description of the time-dependent
chemistry and is configured to co-integrate the chemical and 
physical evolution in Lagrangian streamlines of 100 grid points. These  
are (initially) equally spaced between 1000au (r_ min) and 0.3pc 
(r_ max) with an assumed spherical (1D) symmetry
(corresponding to 0.33 A_ v 40). 
The value of r_ min is chosen so that r>0 at the end of the 
simulation, whilst r_ max corresponds to the physical size of 
the core.

The physical conditions are defined according to the hydrodynmic model of 
KC10 and consist of the density n(r,t), inflow velocity v(r,t), 
gas temperature T_ gas(r,t), dust temperature T_ dust(r,t) 
and the mean visual extinction to the surface of the cloud A_ v(r,t). 
The look-up tables give data as functions of positions and time for 1000 radial 
points (defined in the range 6.2au < r < 0.75pc, corresponding to
0.07 A_ v 200) at 62 times (in the range 
3.05<t<6.15 years). 

For times prior to t_ min=35 years the core evolves very slowly 
and, we take it to be dynamically static so that
we set the parameters to their boundary values, 
i.e. n(r)=n(r,t_ min),v(r)=0, T(r)=T(r,t_ min), and 
A_ v(r)=A_ v(r,t_ min).
In effect, there are two dynamical phases: (i) hydrostatic equilibrium from 
t=0 - t_ min, then (ii) quasi-static contraction from 
t=t_ min to t=t_ max. This is a satisfactory approach, so long as 
the chemical timescale is shorter than the dynamical timescale.
As this is debatable, we have included the possibility of an additional 
dynamical `pause' (t_ pause) at the beginning of the simulation to 
test the assumption.

Gas-phase reactions, freeze-out, various continuous desorption mechanisms 
and surface chemical reactions are all important processes in starless cores.

We include a comprehensive gas-phase chemistry and full descriptions
of the desorption mechanisms, coupled to a fully time and position dependent
analysis of the compositional structure of the ices.
Photochemistry plays an important role in the chemical evolution
of starless cores, both through gas-phase photolysis and the photodesorption
mechanisms, that depend on the spatial
variations of the radiation field.
The outer parts of the core (r 0.07pc, A_ v 3.1) 
are effectively in a photon dominated region (PDR), where the mutual shielding 
of the H_2, C and CO photodissociation and photoionization continua, together 
with dust absorption 
attenuates the radiation field as it penetrates the cloud. 
Allowances for these various attenuation effects are included in the model.
The model incorporates a simplified form of the surface chemistry, explained 
and justified in the appendix.

We follow the chemistry of some 81 gas-phase and 17 solid-state 
chemical species (of limited molecular complexity, containing no more than 
six atoms) composed of the elements H, He, C, N, O, S and Na (which we use as a 
`representative' low ionization potential metal). 

The elemental abundances relative to hydrogen are given in 
Table <ref> and are typical values for low mass star-forming regions
<cit.>. 

The interaction between gas-phase atoms and molecules and the surface of
dust grains is of crucial importance to our understanding of the chemical
structure of molecular clouds and star-forming regions.
Notwithstanding the discussion above concerning the 
physics/dynamics, this is by far the most important source of
ambiguity and possible misinterpretation of the observational data.

In general, we need to consider the processes of freeze-out, or depletion
of a species from the gas-phase, surface chemistry, and the back-process of 
desorption, or sublimation, of the species from the (surface) of dust grains, 
back into to the gas-phase.




With the exception of H_2O, the desorption properties are not well defined for
the dominant ice species, CO, NH_3, CO_2 and CH_4, as well as N_2, 
O_2 and the monomers, C, N and O, from which the hydrogenated species form.

We concentrate on CO and H_2O, whose abundance profiles 
are reasonably well-constrained. Although this could be extended to the related 
species CO_2, which is a major ice component, an in-depth analysis is
not justified due to the lack of empirical constraints for this species 
in L1544.

Desorption of ices by the cosmic ray heating of grains dominates the 
desorption in the darkest parts of the core. However, the rates are uncertain.
Recent studies <cit.> have shown that when realistic dust 
grain cooling profiles are considered the rates are typically an order
of magnitude faster than those determined in previous calculations <cit.>.
The rates may be significantly larger if the effects of sporadic desorption,
ice species co-desorption and/or whole mantle evaporation from small grains 
are included.  
In this study we use the revised rates for a standard interstellar dust grain 
size distribution <cit.> given in Table 3 of <cit.>.
KC10 found that it was necessary to use a CO desorption rate of 
k_ crdes^ CO∼ 3× 10^-13s^-1 to explain the 
high CO abundance that is observationally inferred for the central regions.
We use a value of k_ crdes^ CO=5.9× 10^-13s^-1, 
appropriate for CO molecules with a binding energy of 1100K, and a mean 
adsorbate binding energy of 1200K, <cit.>. This 
value is very close to that determined by <cit.> in their detailed 
random walk Monte Carlo model of the desorption mechanism
<cit.>.



 §.§ Model parameters



For our standard model we match the conditions in L1544, as the prototypical
example of an isolated pre-stellar core. 
The values that we use for the various parameters are given in 
Table <ref>.
Being marginally thermally supported with no internal energy source, the 
physical structure of L1544 is actually relatively simple 
so that these are fairly well-constrained and treated as fixed 
in our physical-chemical model. 

The values of some of the parameters (identified in Table <ref>)
are defined by the chemical abundance profiles and we have considered the 
sensitivity of our results to these (see section <ref>).





Since the various surface processes depend on the grain surface area, 
the rms grain radius (a_ rms) for the grain size distribution
merits discussion. 

The gas-grain interaction and the microphysics of the evolution of
the ice mantles (as described in Appendix <ref>) should be 
followed for grains of different size. Here, as in other studies, we instead 
represent the distribution by mean values of the grain radius, 
cross-section etc.

Conventionally, astrochemical models use a `classical' value of
a=0.1μm. However, the empirically constrained MRN interstellar grain
size distribution is characterised by a power-law:
n(a)∝ a^-3/5 with lower and upper limits of
a_ min∼ 0.005μm and a_ max∼ 1.0μm.
The surface area in this distribution is heavily weighted towards the small 
grain population and the resulting a_ rms∼ 0.011μm. 
These small grains also dominate the UV opacity that regulates the 
photoprocesses.

A larger value of a_ rms would be more appropriate if the effects of 
dust grain coagulation and/or ice mantle accretion are significant, as may be
appropriate in dense, cold gas <cit.>.
However, if this were the case, then there would be several other
implications to consider. These include
(a) the effects on the optical properties of the dust, specifically its
UV opacity - which affects the rates of the A_ v-dependent photoprocesses, 
(b) the effect on the grain surface area per unit volume,
(c) the electrostatic enhancement factor for the freeze-out of ions, and
(d) the relative rate at which dust grains grow, due to the accumulation 
of ice mantles (Rawlings et al., in preparation). 

Nonetheless, for the purpose of this study
and for the sake of consistency and comparison with previous models, 
we adopt the values of a=a_ rms=0.01μm as given in 
Table <ref>, and defer a full discussion of these issues to a 
future study.

The value of the dust surface area per hydrogen nucleon (σ_ H) 
depends on a number of factors, including the assumed value of the dust to gas 
ratio (by mass, d_ m), a_ min, and the mean density of the dust 
grains (ρ_ g). 
With these uncertainties, σ_ H for the MRN distribution lies
in the range 3-8× 10^-21cm^2.
Assuming that d_ m∼ 0.1 (by mass) and 
ρ_ g∼ 3.0g cm^-3 then 
σ_ H∼ 6.6× 10^-21cm^2.
σ_ H is an important parameter in our model and 
we adopt a standard value of 6.0× 10^-21cm^2 
(which is the same as used in KRC14).

The photophysics can be moderated (and the temperature structure modified) by 
either (a) a reduction in the interstellar radiation field strength 
(G_0), or (b) the inclusion of an external extinction (A_ v,ext) 
which would be appropriate if L1544 were surrounded by a 
low density background gas.
However, these are not degenerate; (a) would result in 
a reduction in the physical extent of the PDR, whilst (b) 
could result in the PDR being shifted so that it is completely outside 
the core. 

The initial conditions for the chemistry are taken to be atomic, with 
the H:H_2 ratio set to 10^-3 and the fractional abundance of 
CO, X(CO)_0=10^-6.
As explained above, a more realistic set of initial 
conditions could be obtained by evolving the chemistry to near-equilibrium.
Therefore, to test the sensitivity of the chemistry to the dynamical evolution 
prior to the intial dynamical configuration that we have adopted, we have 
included the possibility of a hydrodynamic `pause' 
(t_ pause) at the start of the simulation. 

During this period the cloud is treated as being static and the physical 
conditions are held constant.



§ RESULTS





 §.§ The chemical structure



Results from the model are shown in Figure <ref> for the values 
of the parameters given in Table <ref>, which are essentially
fixed for the specific chemical structure of L1544.
The figure shows the radial profiles of the fractional abundances of key 
species as functions of position at t=t_ max, i.e. corresponding to the 
current observational epoch. The fractional abundances 
are defined relative to the number density of hydrogen nucleons 
(n_ nuc = 2n_ H_2 + n_ H). This differs to the 
definition used in some other studies (e.g. KRC14), which is the relative 
to the H_2 number density. 

From various model calculations, and Figures <ref> and 
<ref>, we find that it is possible to identify several different
chemical regimes.
From the outside, moving inwards:


  * 
For the gas densities applicable to the outer parts of the 
core, the transition between C^+/C/CO, characteristic of a photon 
dominated region (PDR) typically occurs at a mean A_ v∼ 1-3 
magnitudes. In our model, this corresponds to a radial position of 
between ∼ 0.07-0.145pc, where T_ dust∼ 10.2-12.9K.
Even at the outer boundary of 0.75pc (A_ v∼ 0.07), 
T_ dust=17.8K.
At these very low temperatures thermal desorption from dust grains is 
negligible for nearly all chemical species of interest, implying the 
existence of a `cold PDR' with a rather unusual gas-grain chemistry.
We identify two zones within the PDR boundary layer:



  * the ice free zone; where T_ dust 15K 
(r 0.34pc, A_ v 0.28). 
Here freeze-out and surface hydrogenation are inhibited by photodesorption, 
and the chemistry is dominated by gas-phase processes, and

  * the (marginal) ice formation zone; where T_ dust 15K,
freeze-out and surface hydrogenation occurs (yielding CH_4 
and H_2O etc.), but is limited by photodesorption. 
Gas-phase chemistry still dominates in this region.


  * Interior to the PDR boundary layer there is the `dark core' with three
zones;


  * an outer zone where the gas temperature falls rapidly, although the 
density does not rise as fast, so that freeze-out is once again inhibited,

  * an intermediate zone where the gas density rises significantly and 
freeze-out and photodesorption are both important; this is the location 
of the molecular (H_2O) abundance peaks, and  

  * an inner zone, where the density and extinction are high and 
freeze-out dominates <cit.>.














The abundance profiles of the various species shown in Figure <ref> 
are in broad agreement with those that are inferred from the L1544 observations.
For H_2O the principle empirical constraint on the gas-phase abundance 
profile comes from the radiative transfer analysis of the H_2O 567GHz 
(1_10-1_11) line as observed by Herschel <cit.>.
The different chemical zones discussed above are clearly evident, with 
substantial ice mantle formation only occurring within the inner 
∼0.08pc (A_ v 2.6).
Specifically, the H_2O abundance profile is consistent with what was 
presented in KRC14. 
The radial location of the abundance peak is shifted purely as a result 
of the use of an updated A_ v(r) profile.
The NH_3 abundance profile remains high until 
innermost depletion region (r 0.02pc, A_ v 12), 
consistent with observations <cit.>.








  §.§.§ Ice composition



In Figure <ref> we show the compositional structure of the ices on the 
dust grains (in STARCHEM represented by particles of a single size). 
The figure shows the fractional composition of
the ice for each monolayer, at two positions; near the outer edge of 
the ice formation zone (where photochemistry is important) and close to 
the (dark) centre of the core.

We can see that there are variations in the ice composition as a 
function of position in the core and that there are very significant 
stratifications within the ice mantles. 
In conditions where significant freeze-out occurs, the sub-surface layers of 
the ices - which are effectively chemically inert - are determined by the entire
evolutionary history of the cloud/core. 
Whilst, as we show in section <ref>, the gas-phase abundances evolve 
quasi-statically, the ice has a long memory so that each monolayer retains the 
abundances that were in equilibrium with the gas phase abundances at their time 
of deposition. The compositional stratification of the ice mantles thus reflects 
the chemical evolution of the gas. The molecules H_2O and CO are seen 
to dominate the composition of the outer layers, as we would expect, because 
H_2O and CO are the dominant oxygen and carbon species inside the PDR. 
Moreover, the O and OH that freeze onto the grains are quickly hydrogenated to
H_2O. In contrast, in the inner layers of the ice mantles,
ammonia (NH_3) and, most especially, methane (CH_4) are major components.
The source of the CH_4 is incomplete conversion of C to CO in the earliest 
stages of the chemical evolution, resulting in the freeze-out and hydrogenation 
of C (see section <ref> below). 
A similar behaviour can explain the NH_3 and N_2 abundance 
stratification; the inner mantle samples the early time chemistry, when the gas 
has a large atomic component. Here, accreting N is converted to NH_3. 
The outer mantles sample the late time chemistry, when efficient conversion 
of N to N_2 has taken place, and there is little N in the gas-phase.

These results are quite different to those obtained by <cit.> with
the GRAINOBLE model. This is due to a number of factors; <cit.> 
adopt a different grain size distribution, omit photodesorption and desorption
driven by the enthalpy of formation and - most significantly - they only included
the accretion of H, O and CO, so that the formation of CH_4 and NH_3 are
inhibited.

Unlike some studies <cit.>, we do not consider
the differences in binding energies for bare and icy grains, other than for the
desorption driven by the enthalpy of molecule formation. Whilst this has some
implications for the cosmic-ray heating-induced desorption of the most volatile 
species, this is generally not significant as the dust temperatures are so low
that thermal desorption is negligible for nearly all species of interest.

Our results can be used to compare with recent JWST observations of ice
mantle compositions toward dense dark clouds <cit.> and to guide 
future JWST observations of pre-stellar cores.






Small variations of the dust parameters (within the range of uncertainty)
can result in an
interesting phenomenon in the behaviour of the solid-state abundances of both 
CO and H_2O in the `cold core' region. These sometimes 
show strong enhancements in the range 0.18 r 0.24 pc
(0.44 A_ v 0.7).
Although the specific form of this feature is an artefact of the physical parameters
in our model (e.g. the simplifying assumption that dust grains can be represented by
a single size population), there is a genuine - if somewhat counter-intuitive - reason
for its existence: the sharp rise in the abundances at the outer positions correspond 
to the formation of ices of more than one monolayer thickness. However, whilst the 
desorption processes are only weakly dependent on local conditions, the freeze-out 
rates are proportional to the density (n) times the square root of the gas 
temperature (T_ gas). 
There is a region in the inflow where T_ gas^0.5 falls faster than n rises 
and consequently the net freeze-out rate falls and the ices reduce to the 
sub-monolayer regime (within r 0.18pc, A_ v 0.7). 
In a real grain size population this enhancement would be smoothed out, but the 
overall effect is for dust grains to be relatively ice-free (i.e. in the 
sub-monolayer regime) right down to r∼0.09pc (A_ v∼ 2.2).



  §.§.§ The chemistry of H_2O






To understand the behaviour of H_2O in more detail, we need to identify 
the dominant formation and destruction channels.
These have been discussed at some length 
elsewhere <cit.> but to summarise; in the cold, dense environment 
of pre-stellar cores the predominant gas-phase formation route is driven by 
(cosmic-ray generated) H_3^+:



    O + H_3^+ → OH^+


    OH^+ + H_2 → H_2O^+ + H


    H_2O^+ + H_2 → H_3O^+ + H


    H_3O^+ + e^- → H_2O + H

Neutral-neutral reactions, via the hydrogenation of the OH radical,
are only significant at higher temperatures, and







in cold environments the H_2O formation is 
dominated by surface reactions on interstellar ices. 
Several pathways have been identified <cit.>, 
the most direct being the surface hydrogenation of accreted 
oxygen atoms and hydroxide (OH) radicals.
In addition, and particularly in relatively dense regions, significant 
contributions may be made via alternative 
(interconnected) channels involving O_2, either accreted from the gas-phase, 
where it is primarily formed
by the OH+O reaction, or else created via oxygen atoms reacting on the 
surface. The hydrogenation of surface O_2 leads to the formation of 
intermediate species en route to H_2O formation, such as 
HO_2 and H_2O_2. Alternatively, surface oxidation leads to the 
formation of O_3, which can be hydrogenated, again leading - ultimately 
- to H_2O formation <cit.>.

Figure <ref> shows different abundances in each
of the chemical zones identified in section <ref>; 
low molecular abundances are seen in the outer PDR region where 
photodissociation is active and also in the inner core where freeze-out
dominates. The molecular abundances peak at intermediate extinctions 
and densities.

In the outer regions (r 0.2pc, A_ v 0.6) most of the 
gas-phase H_2O derives from ice that has been desorbed by photodesorption 
and the enthalpy of surface formation of H_2O from O and OH. 
Indeed, if the latter process were suppressed, then the H_2O abundance would 
be lower by a factor of up to ∼ 10×.
The main H_2O loss mechanisms are photodissociation and reaction with C^+ 
(to form HCO^+, ultimately leading to CO formation). 
This effectively inhibits ice formation until r 0.1pc 
(A_ v 1.9).





At the location of the H_2O abundance peak (∼ 0.08pc, 
A_ v∼ 2.6) photodesorption and formation enthalpy-driven desorption
still dominates the production of gas-phase H_2O, although gas-phase 
reactions (via H_3O^+) also contribute significantly. 
Photodissociation is the dominant destruction mechanism.
In the inner core cosmic-ray induced photodesorption dominates the production 
of gas-phase H_2O. As the rates for this are rather 
ill-defined there is a considerable margin of error in the model predictions.
In the innermost regions (r 0.02pc, A_ v 12) the 
abundance declines due to efficient freeze-out at high densities - but gas-phase 
destruction by H_3^+ accounts for ∼ 30% of the loss rate. 
The surface hydrogenation of accreted OH is the main formation channel.














  §.§.§ The chemistry of CO



Unlike H_2O, which is mainly formed on the surface of grains, CO
is formed in the gas-phase.
In dense cloud conditions this chemistry is relatively slow, 
compared to other gas-phase processes and freeze-out, but may result in 
the efficient conversion of C to CO.

Several channels are possible, including:

    C^+ + H_2 → CH_2^+ + hν


    CH_2^+ + O → HCO^+ + H


    HCO^+ + e^- → CO + H

or

    C + H_2 → CH_2 + hν


    CH_2 + O → CO + 2H


However, the initiating radiative association reactions for both of these channels
are very slow;
the reactions of C and C^+ with H_2 to form CH and CH^+ have large
activation barriers, whereas the corresponding radiative association reactions
(to form CH_2 and CH_2^+) do not possess inhibitive barriers, but are inefficient.

In normal dark cloud environments, reactions involving OH may therefore be
more important:

    C + OH → CO + H

or

    C^+ + OH → CO^+ + H


    CO^+ + H_2 → HCO^+ + H

followed by reaction (<ref>),
although the efficacy of reactions (<ref>) and (<ref>) have 
recently been challenged by the non-detection of CO^+ in diffuse 
interstellar clouds <cit.>, which suggests that non-thermal reactions involving
C^+, CH^+, CH_2^+ etc. and oxygen atoms are likely to be a more 
important source of the observed HCO^+:

    C^+ + H_2 → CH^+ + H


    CH^+ + H_2 ...→ CH_2^+,CH_3^+... + O → HCO^+ + ...


But, as noted above, the first of these reactions is very slow at low 
temperatures. 

Interior to the PDR, at the densities applicable to pre-stellar cores, the 
conversion of C to CO goes to near completion (i.e. X(CO)∼ X(C_ total))
within a few× 10^5 years
and the ratio of the gas-phase to solid-state CO abundances is set by 
the balance of freeze-out and desorption.
CO has a much lower surface binding energy than H_2O so that the 
gas-phase abundance is sensitive to the assumptions concerning the
desorption efficiencies. 
Thus, although the dust temperature is low enough that thermal desorption
is inhibited throughout this region, desorption by cosmic-ray heating
is important (and dominates over photodesorption). The 
effect of chemical desorption (driven by H_2 formation) may also be 
significant <cit.>.

Further out, within the PDR zone, CO formation is incomplete and the 
abundances of C and C^+ are relatively high. Indeed, even deep inside 
a PDR, the abundances of C and C^+, do not vanish and are typically 
of the order of 10^-4-10^-3 times that of CO 
<cit.>.

The situation in a `cold PDR' is more complex; in (the oxygen rich) 
gas-phase, carbon chemistry tends to lead to the formation of the stable CO 
molecule as the dominant product (leaving free volatile oxygen to form 
H_2O, O_2 and other species), with a much smaller proportion of
the available carbon present in other organic molecules.
However, on the surface of grains, it is more likely that any accreting
carbon atoms or ions are hydrogenated to form methane (CH_4).

If the conditions are propitious for freeze-out (as is 
the case in the cold PDR in the outer regions of the core), then there will 
be a competition between 

    C  or  C^+ + H_2  or  OH ...→ CO


    CO + grain → CO(s)

also leading to methanol (CH_3OH) formation (which we do not discuss
in this study; see Appendix <ref>), and

    C  or  C^+ + grain → CH_4(s)

The very low temperatures and the chemical structure of pre-stellar cores 
imply that the gas-phase and gas-grain chemistries pull in different 
directions (leading to CO and CH_4 respectively).
This implies that the efficiency of CO formation in the pre-contraction phase 
is important in determining the balance between CO/CH_3OH and CH_4.

This balance will be strongly dependent
on the assumptions concerning the grain size distribution. Thus, if the 
grain surface area is predominantly carried by negatively charged very small 
grains (as is the case in the MRN dust size distribution), then there could be 
a strong electrostatic enhancement of the freeze-out of C^+, favouring the 
methane formation channel. 
We defer a full discussion of the importance of the grain size distribution 
to a future study.

Referring to Figure <ref> we can see that 
the qualitative behaviour of the CO abundance profile is similar to that for
H_2O. That is to say there are essentially three chemical regimes;
(i) an outer PDR zone, where X(CO) tails off at low A_ v,
(ii) a middle region where X(CO) reaches near saturation, and
(iii) an inner depletion zone where the depletion factor depends on the ratio
of the freeze-out to the cosmic-ray induced desorption rates.

We confirm the findings of KC10; that an efficient desorption mechanism is
required to maintain sufficiently high CO abundances at high 
extinctions, to explain the observed line strength. 
The most efficient mechanism for CO desorption is cosmic-ray grain heating - 
but, as deduced by KC10 from a comparison with the observed brightness of CO
in L1544, the rate would need to be ∼ 10-30× larger than given 
in <cit.> [The uncertainty of a factor of three in this estimate
comes from the non-linear scaling of the CO brightness with the cosmic-ray 
desorption, the accuracy of the radiative transfer modelling and the 
observational uncertainty].
Current studies <cit.> indicate that cosmic-ray heating induced 
desorption rates may have been grossly underestimated in previous studies
<cit.>.
There are, however, very significant uncertainties in the quantification of
these rates and their A_ v-dependence but these, empirically deduced
values, are within this range of uncertainty.

In our model we find that the formation of CO is a little more complex than 
H_2O; in the outer regions several 
formation channels are effective, including dissociative recombination of
HCO^+, but also via the photodissociation/hydrogen abstraction of 
surface-formed CH_4, yielding CH, which reacts with oxygen atoms. 
A corollary of this pathway multiplicity is that the CO abundance profile is 
fairly robust to uncertainties/variations in the rate coefficients for specific
reaction pathways.
The destruction of CO is completely dominated by photodissociation (in
the PDR) and freeze-out (in the inner core). At intermediate radii
(0.08pcr  0.15pc, 0.96 A_ v 2.6) 
reaction with H_3^+ to form HCO^+ is 
significant, although most of the HCO^+ is re-cycled to CO.
In the innermost regions, there is little gas-phase chemistry and the abundance
is essentially established by the balance between (cosmic-ray heating) desorption
and freeze-out.



  §.§.§ CO and the nitrogen chemistry



Some previous studies <cit.> have predicted gas-phase CO to be 
concentrated in a `ring' at A_ v∼1.
Part of the reason for this is due to the coupling of the chemistries for 
CO and H_2O (e.g. as is evident from reactions <ref> and 
<ref> above). But, most significantly, it was postulated that in 
the regions where H_2O freezes out but (the more volatile) CO does not, CO
is destroyed (by He^+, N_2H^+ and photodissociation at low A_ v) 
yielding O.

This can then freeze out and be converted to H_2O on grains,
so that the net result is a gradual shift of the oxygen from CO 
to H_2O, and an anticorrelation between CO and N_2H^+.

Our model results do not show this anticorrelation. Observationally,
there is some evidence of a rise in the abundance of 
N_2H^+ in high density regions where CO is depleted 
<cit.>, although this observation could be an artefact of
radiative transfer effects <cit.>. 
We defer a full description of the chemistry of N_2H^+ to a future
study, but our results indicate that the abundance profiles of CO and N_2H^+ 
are only weakly connected.

We also find that destruction by He^+ is insignificant at all positions. 

The nitrogen chemistry is fairly simple:
throughout much of the core, the chemistry of N_2 is driven by CH:

    CH + N → CN


    CN + N → N_2

and, here, reactions with He^+ are the dominant destruction mechanisms 
(ultimately providing formation channels for NH_3 and N_2H^+), although 
reaction of N_2 with H_3^+ to form N_2H^+ is significant at the 
innermost positions.
In these regions of high depletion, the hydrogenation of N^+ to NH_4^+, 
followed by dissociative recombination, is also a major formation 
channel for NH_3.
As with CO the dominant desorption mechanism for other volatile species, 
such as N_2, is found to be direct cosmic ray heating of grains.






 §.§ Time-dependence



In Figure <ref> we show the radial abundance profiles for selected
species (H_2O, H_2O(s), CO, CO(s), NH_3 and CS) obtained at three 
specific times in the dynamical/chemical evolution: 
t_ min=0.30Myr, 0.45Myr and t_ max=0.61 Myr. 






We can immediately see from this figure that there is hardly any `intrinsic' 
time-dependence of the chemical abundances. The abundances do change with 
position and time, particularly in the inner regions, but this is primarily a 
result of the contraction of the core and the adjustment of the chemical 
equilibrium to the changing physical conditions. 

To test this result we have re-calculated these profiles, but 
setting t_ pause=0.5Myr - i.e. inserting a period of chemical 
evolution, during which the physical parameters are held constant, at the
beginning of the simulation. This yields almost identical results to those
shown in Figure <ref>, other than a small reduction (by a factor of
2 in the abundance of NH_3 in the range 0.15pc<r<0.25pc,
0.42 A_ v 0.96).
This result indicates that it would not be possible to determine a 
`chemical age' for the cloud.

Although the chemical evolution 
is in quasi-equilibrium, in the parts of the core where thick ice mantles are 
present it is possible that some abundances will be affected by the prior
evolution of the core. This is because the equilibrium is determined by the balance
between gas-phase and gas-surface interactions. In our model the previously 
accreted sub-surface ice layers are assumed to be inert and effectively removed 
from the active chemistry.

With the assumption of quasi-equilibrium, we can estimate the gas-phase 
abundance of H_2O in the inner parts of the core as balance between freeze-out 
and desorption. We can also calculate the relative contribution from each of
the applicable desorption mechanisms; cosmic ray induced photodesorption, cosmic ray
heating induced desorption and H_2-formation induced desorption.

Comparing the desorption efficiency for cosmic-ray induced photodesorption
with that for H_2 formation-induced desorption, we find that the two are 
equal when

    (Y_H_2/Y_pd) = (0.85/n_H)
    (ζ/1.3× 10^-17s^-1).

Thus, in normal circumstances, this implies that values of 
Y_ H_2 10^-3 would be needed for H_2 formation-induced 
desorption to be important.

Then, comparing the desorption efficiencies for cosmic ray heating-induced desorption 
with cosmic ray induced photodesorption we find that the former dominates when 

    k_cr > ( Y_pdF_cr/N_s) 
     5×10^-15s^-1

where we have used parameter values as given in Table <ref>.
For species of interest, most estimates of k_ cr are less than this value, 
but we cannot exclude the possibility that recent revisions could imply that 
cosmic ray heating may the dominant desorption mechanism for most species in dark
environments.

Hence, balancing freeze-out with cosmic-ray induced photodesorption, and with the 
additional assumptions that T_ gas=10K and that the proportion of surface
coverage of grains by water molecules is of order unity (as indicated in 
Figure <ref>), then - again using parameter values from 
Table <ref> - the equilibrium gas-phase density of H_2O is

    n_H_2O = 1.8× 10^-3( Y_pd/10^-3)
    ( ζ/1.3× 10^-17s^-1) cm^-3.

Assuming that Y_ pd and ζ do not change as functions of 
position, this is independent of the local conditions, so that the fractional
abundance of water, X(H_2O), simply scales as the inverse of the density.
This behaviour can be seen in the inner regions (r 0.05pc, 
A_ v 4.8) in Fig <ref>. 



  §.§.§ Parameter variations



We have also investigated the sensitivity of these results to the values 
of some of the parameters (identified in Table <ref>), the chemical 
initial conditions and the choice of reaction database, within the 
ranges that are appropriate for pre-stellar cores. 
We defer a full discussion to a future paper but highlight some key
results here.

Apart from some obvious dependencies on the cosmic ray induced desorption
rates and the various desorption efficiencies (particularly in the 
inner core) the H_2O and CO abundance profiles are, in general, found 
to be fairly insensitive to most of these parameter variations. For most
species, a range of gas-phase, gas-grain and surface reactions operate at 
all positions within the core.

Changing the initial conditions (specifically, the initial fractional 
abundance of CO) has little effect on the gas-phase abundances - 
consistent with the condition of quasi-static equilibrium - but 
results in significant variations in the CH_4 to CO ratio in the ice 
mantles. This is an interesting result that will be investigated in more 
detail in future studies. 

The chemistry is, however, very sensitive to the value of the external 
extinction, as this affects both the radiation field strength and
the location and extent of the PDR. We find that the inclusion of even 
modest values of A_ v,ext.∼ 1-2 results in a shift in 
the location of the H_2O abundance peak and large abundances in the 
outer envelope that are quite incompatible with the empirically
constrained abundance profile. 

There is filamentary structure in the vicinity of L1544, and the presence 
of a parsec-scale low density (∼ 30cm^-3) medium is required to 
reproduce the HCO^+(1-0) observations <cit.>.
However, the regions that are responsible for the observed CO and H_2O 
line profiles cannot be obscured from the interstellar radiation field 
by material with significant extinction in the near-IR and UV.



 §.§ The applicability of reduced chemical networks



In this section we compare the accuracy of using a highly reduced chemical 
reaction network relative to a full description of the gas-grain chemistry.

In a cold, dark cloud most of the oxygen budget is taken up by O, O_2, 
CO, OH and H_2O. As CO is relatively chemically inert and O_2 is not a 
dominant component, KRC14 proposed that the H_2O 
chemistry could be described adequately by a highly reduced network of just
a few reactions. In this scheme the gas-phase chemistry is limited to 
the hydrogenation and photodissociation reactions:

    O + H_2 → OH + H


    OH + H_2 → H_2O + H




    H_2O + hν→ OH + H


    OH + hν→ O + H

Although the hydrogenation reactions are vanishingly slow at temperatures 
below a few hundred degrees, they are included to provide reverse reactions 
for detailed balance with the photodissociation reactions. With these reactions, 
the  reduced chemical network is also applicable in warmer gas.

The gas-grain interactions simply include the freeze-out
and (assumed complete) surface hydrogenation of O,OH (and H_2O):

    O, OH, H_2O → H_2O_(ice)

balanced by desorption:

    H_2O_(ice)→ OH and H_2O

The processes include direct 
cosmic ray heating, photodesorption and (in KRC14) cosmic-ray 
induced photodesorption. Thermal desorption, though included, is 
insignificant in starless cores because of the low dust temperatures.

KC08 and KC10 assume that the carbon chemistry could be decoupled 
from the oxygen/H_2O chemistry and be reduced to just three reactions:

    C + hν→ C^+ + e^-


    CO + hν→ C + O

and the (slow) radiative association of C^+ with H_2 
(reaction <ref>), which is the rate-limiting step for CO formation 
(followed by reactions <ref> and <ref>). 
The only gas-grain processes were taken to be the freeze-out and (cosmic 
ray-induced) desorption of CO.

















  §.§.§ Comparison of reduced with full chemical networks
 


As mentioned in the introduction, one advantage of the reduced chemical networks 
is the simple and rapid calculation of the abundances of the dominant coolants 
using a few algebraic equations. These apply to each of the H_2O and CO 
chemistries and relate the abundances to the values of the model-dependent 
reaction rates. This allows for a simple analytical approximation to establish 
the temperature of the gas <cit.> which is valid even across PDR boundaries.

Previously we compared our dynamical and chemical model of L1544 to observations 
of L1544 (KC10, KRC14). The predicted C^18O and H_2O abundance profiles for 
L1544 along with radiative transfer modelling produced C^16O, C^17O and 
C^18O (1-0) and H_2O 567GHz (1_10-1_11) line profile shapes and 
strengths consistent with those observed. The reduced chemical 
network is also responsible for the gas temperature which factors 
into the line strengths. In this sense the method is observationally verified. 

A direct comparison of the abundances produced by the reduced chemical network 
with those from a full network has the advantage that the accuracy of the 
approximation can be determined at least for the examples compared. 
To assess the accuracy of the reduced network chemistry, we have compared the 
abundances at t_ max from our time-dependent full network chemical model 
with the equilibrium abundances at t_ max also from STARCHEM but including 
only the reactions from the reduced network as defined above. 
This effectively mimics the analytical approach used in our earlier papers, 
KC08, KC10, and KRC14.

The results from this comparison show that whilst there is a good fit to the 
H_2O abundance profile, the fit for the CO profile is not as good; the 
reduced network results in too much CO in the outer parts (due 
to the omission of various gas-phase destruction channels) and too little CO
in the inner regions (due to the slowness of the CO production reactions).
The C and C^+ profiles are also problematic, with high abundances persisting
down to the centre of the core. This is again partly due to the slowness 
of the C to CO conversion reactions in the reduced network.
 
An obvious problem with the carbon chemistry in the reduced network as proposed by 
KC08 

is that C^+ recombination was not included:

    C^+ + e^- → C + hν
 
and, other than photoionization, there is no loss channel for carbon atoms.
We have included additional channels in the network to account for the conversion 
of C and C^+ to CO (as suggested in section <ref>):


    C + OH → CO + H


    C^+ + OH → CO^+ …→ CO

(assuming that every CO^+ molecule so-formed leads to CO formation, via 
reactions <ref> and <ref>).
Here, the abundance for OH is found from the reduced network for oxygen. 
The loss of OH in these reactions should not much affect the abundances in the 
oxygen network since the total abundance of oxygen is three times that of carbon, 
and we ignore this effect.

In Figure <ref> we show the results of a comparison between 
the full network and our reduced networks for the H_2O chemistry
(as in KRC14) and the CO chemistry (revised, as described in 
Appendix <ref>).
The figure shows the abundance profiles (at t=t_ max)
obtained for H_2O, H_2O(s), O, C, C^+, CO and CO(s.).

To compare `like with like' we suppress enthalpy-driven desorption and the 
surface chemistry of CO (leading to CO_2 formation etc.) in these 
simulations.

In this figure, we also show the abundances determined directly from the 
algebraic equations for chemical equilibrium with the reduced networks.






Using the revised reduced networks in the full model yields remarkably 
accurate results (typically within a factor of ∼ 2×, at all radii) 
for all species, other than CO and CO(s.) in the innermost region where 
there is some discrepancy. This is due to the large methane abundance 
obtained in this partcular model, which is not accounted for in the reduced
network.

The algebraic solution also gives excellent fits, except for CO and 
H_2O which are over-estimated in the inner regions, where the dust grains
have thick ice mantles. The reason for this is straightforward; the desorption 
rate of a species depends on the fraction of the grain surface that it covers.
The algebraic formulae take the desorption rate to be proportional to the total
ice mantle abundance of that species. In the sub-monolayer regime this is a
satisfactory approximation, but for thicker ices this over-estimates the 
desorption rate by a factor of order of the number of monolayers. 
This discrepancy has only a marginal effect on the ice abundances as they
are much greater than the gas-phase abundances in these conditions.

We conclude that the reduced networks (including the modifications for the CO
chemistry as described in Appendix <ref>) work extremely well, 
although their accuracy in regions of high depletion are contigent on 
H_2O and CO being the dominant ice mantle components.
We recommend that the revised CO chemistry network be used in future 
applications.
The algebraic solutions to the networks also work well, except for CO and 
H_2O in regions of high depletion.



  §.§.§ Utility of reduced network chemistries
 


In pre-stellar cores the main heating mechanism for the gas throughout most of 
the core is cosmic ray heating.


In the outer region, the UV flux on dust also heats the gas through the 
photoelectric effect. In the central region which is completely shielded from 
starlight, only cosmic ray heating is effective, and the gas becomes very cold, 
approaching 5K. 

The main coolants for the atomic gas are [CII] and [OI] fine structure lines, 
whilst the molecular gas cools primarily through the CO rotational lines and 
dust thermal emission. The 
most effective optical depth for cooling is about unity. Where the CO(1-0) is 
optically thick, the higher excitation lines become effective coolants as well 
as a few less abundant simple carbon species <cit.>. 
Figure <ref> shows a local maximum in the gas temperature at a 
radius of 0.03 pc. This is due to the variation of the optical depths of the 
molecular line coolants, primarily CO(1-0). As the dust temperature increases 
with radius due to the increasing UV flux, so does the gas temperature by 
collisional coupling with the dust. However, at radii ∼0.03 pc, the 
optical depth to the boundary of the cloud decreases enough that the cooling 
through the molecular lines becomes more effective. At radii > 0.3 pc, the 
gas temperature decreases despite the steady increase in dust temperature. 
The rapid rise in the gas temperature at the boundary is caused by hot 
electrons photoelectrically released from dust grains by FUV radiation that 
penetrates only a short distance into the cloud.



§ SUMMARY AND CONCLUSIONS



We have applied an empirically constrained multi-point dynamical model, coupled 
to a comprehensive description of the chemistry, to a study of the chemical 
processes and timescales in pre-stellar cores. We have applied this model to
L1544 to study the temporal and spatial evolution of the abundances of H_2O, 
CO and other key small molecular species.

The key results from our studies are as follows:



  * At the time corresponding to the epoch of observation of L1544, the chemical 
timescales are less than the dynamical timescale throughout the core. Hence, the
chemistry is in near equilibrium; that is to say, the chemical abundances at 
any position simply respond to the instantaneous local physical conditions 
(n,T,A_ v) and are, to a large extent, independent of the prior chemical 
and dynamical evolution. 
The H_2O abundance profile in this type of source only depends on the 
physical structure, and the values of certain chemical parameters.

This is a significantly different result to what is usually assumed to be the
case for young star-forming regions and implies that, for isolated pre-stellar 
cores at least, the chemical and dynamical 
evolution are not separable. Hence, the concept of a 
chemical (or depletion) `age', and the applicability of `chemical clocks',
are meaningless.
The only situations in which chemical clocks may be relevant are the 
highly unusual conditions where an observed system has assembled on a 
timescale that is short compared to the chemical timescale, but then evolves 
on a very much longer timescale.


  * The ice mantles on dust grains are highly structured, so that whilst 
CO and H_2O dominate the outer layers other species, such as CH_4 and 
NH_3, are important in the inner layers. 
As has been previously noted <cit.>, 
if it is assumed that the sub-surface ice layers are chemically inactive 
they act as a fossil record/reservoir of the prior dynamical/chemical evolution.
This implies that the quasi-equilibrium gas-phase abundances are not 
completely decoupled from the evolutionary history of the core.
For example, the degree of conversion of C to CO in the earliest 
evolutionary stages determines the overall abundance of CH_4 ice (due to
the competition between gas-phase chemistry; yielding CO, and surface chemistry;
yielding CH_4). 
A similar result applies to NH_3 ice, whose abundance is determined by
the degree of conversion of N to N_2 in the gas-phase.


  * The model matches the empirically determined abundance profiles for 
H_2O and CO in L1544. 
The chemistry of hydrogenated species, such as H_2O, NH_3 
and CH_4, is dominated by surface reactions. For many other species there
are gas-phase and surface chemistry formation channels that operate
throughout the range of physical conditions that exist within the core, so 
that their abundance profiles are remarkably robust to parameter variations.


  * The highly simplified reduced network chemistries for oxygen species 
described in KRC14, and for C/C^+/CO, described in KC08 & KC10 and revised 
here (Appendix <ref>), give a very accurate description of the 
chemistry and abundance profiles for gas-phase and solid-state H_2O and CO, 
OH, O, C and C^+. This is valid over the whole range of extinction 
(1<A_ v<1000) applicable to starless cores and quiescent molecular gas. 
The accuracy of the CO and H_2O profiles in regions of high depletion depends 
on these species being the dominant components of the ice mantles.
Algebraic solutions to these networks also give good descriptions of the abundance
profiles, but over-estimate the gas-phase abundances of CO and H_2O when 
significant freeze-out has occurred.















  * The abundance profiles for both H_2O and CO in L1544 necessitate the 
existence of a PDR in the outer regions of the core and the interstellar
radiation field must reach the outer layers without attenuation 
(i.e. A_ v,ext∼ 0). This is indicative of the presence of highly 
non-uniform density and extinction structures in this source.




  * The chemical structure is quite complex, with both gas-phase and
surface chemistry important throughout the core. In the outer regions `cold PDR'
conditions prevail, with photodesorption and enthalpy-of-formation driven 
desorption dominating the gas/ice balance.
Thermal desorption is only significant in the outermost parts of the core,
and even there is only limited to highly volatile species. 


  * The chemistry in the inner (dark) regions of the core 
(r 0.08pc, A_ v 2.6) is dominated by 
cosmic-ray induced processes. In particular, 
the desorption of ices by cosmic ray induced photodesorption (for H_2O) and 
cosmic ray grain heating (for CO etc.) are the main factors that determine the 
gas-phase abundances in the highly depleted regions.
However, in the light of recent re-calculations <cit.>, we find 
that in these conditions, desorption driven by the cosmic ray heating of grains 
is probably more significant than cosmic-ray induced photodesorption
for a number of volatile species, including CH_4 and N_2.





Future modelling efforts will include updates to include a more 
comprehensive description of the (position-dependent) grain size 
distribution and its several effects on the chemistry. 
This update would also include allowances for the ice 
composition-dependences of the binding energies and desorption 
efficiencies.

Our modelling of the H_2O and CO abundance profiles in 
a specific pre-stellar core (L1544) can be extended to study the 
abundance profiles of, and correlations between, other key molecular 
tracers, such as N_2H^+, N_2, HCO^+, H_2CO and CH_3OH.
The flexibility of the model is such that this could easily be adapted 
and applied to studies of other pre-stellar cores.



§ DATA AVAILABILITY


The data underlying this study are openly available from the published
papers that are cited in the article.
The data generated in support of this research are partly available
in the article, and will be shared on reasonable request to the corresponding
author.







[Aikawa et al.(2001)]Aik01
Aikawa Y., et al., 2001, ApJ 552, 639.

[Andersson et al.(2006)]And06
Andersson S., Al-Halabi A., Kroes G.-J., van Dishoeck E.F., 
2006, J. Chem. Phys., 124, 064715

[Bate & Keto(2015)]BK15
Bate M.R., Keto E.R., 2015, MNRAS, 449, 2643

[Bergin & Tafalla(2007)]BT07
Bergin E., Tafalla M., 2007, Ann. Rev. A&A, 45, 339

[Bizzocchi et al.(2014)]Bizz14
Bizzocchi L., Caselli P., Spezzano S., Leonardo E., 2014, A&A, 569, A27

[Bringa & Johnson(2004)]BJ04
Bringa E.M., Johnson R.E., 2004, ApJ, 603, 159

[Buch & Zhang(1991)]BZ91
Buch V., Zhang Q., 1991, ApJ, 379, 647

[Caselli et al.(1999)]Cas99 
Caselli P., Walmsley C.M., Tafalla M., Dore L., Myers P.C., 
1999, ApJ, 523, L165

[Caselli et al.(2012)]Cas12 
Caselli P., Keto E., Bergin E.A. et al., 2012, ApJL, 759, L37


[Caselli et al.(2017)]Cas17 
Caselli P., Bizzocchi L., Keto E. et al., 2017, A&A, 603, L1


[Caselli et al.(2019)]Cas19 
Caselli P., Pineda J.E., Zhao B. et al., 2019, ApJ, 874, 89


[Caselli et al.(2022)]Cas22 
Caselli P., Pineda J.E., Sipilä et al., 2022, ApJ, 929, 13

[Cecchi-Pestellini & Aiello(1992)]CPA92
Cecchi-Pestellini C., Aiello S., 1992, MNRAS, 258, 125


[Chacón-Tanarro et al.(2019)]Chac19 
Chacón-Tanarro A., Pineda J.E., Caselli P. et al., 2019, A&A, 623, A118

[Chiar et al.(2011)]Chiar11
Chiar J.E., Pendleton Y.J, Allamandola L.J. et al., 2011,
ApJ, 731, 9

[Crapsi et al.(2007)]Crap07 
Crapsi, A., Caselli P., Walmsley M.C., Tafalla M., 2007, A&A, 470, 221

[Cruz-Diaz et al.(2018)]Cruz18
Cruz-Diaz G.A., Martín-Doménech R., Moreno E., Muñoz Caro G.M., 
Chen Y-J., 2018, MNRAS, 474, 3080

[Cuppen et al.(2010)]CIRL10 Cuppen H.M., Ioppolo S., Romanzin C., 
Linnartz H., 2010, Phys. Chem. Chem. Phys., 12, 12077

[Galli et al.(2019)]Galli19
Galli P.A.B., Loinard L., Bouy H. et al., 2019, A&A, 630, A137

[Garrod & Pauly(2011)]GP11
Garrod R.T., Pauly T., 2011, ApJ, 735, 15

[Hasegawa & Herbst(1993)]HH93
Hasegawa T.I., Herbst E., 1993, MNRAS, 261, 83

[Gerin & Liszt(2021)]GL21
Gerin M., Liszt H., 2021, A&A, 648, A38

[Giers et al.(2023)]Giers23
Giers K., Spezzano S., Caselli P. et al., 2023, A&A, 676, A78

[Glover & Clark(2012)]GC12
Glover S.C.O., Clark P.C., 2012, MNRAS, 421, 116

[Goldsmith(2001)]Gold01
Goldsmith P.F., 2001, ApJ, 557, 736

[He et al.(2015)]He15
He J., Shi J., Hopkins H., Vidali G., Kaufman M.J., 2015, ApJ, 801, 120

[Herbst & Cuppen(2006)]HC06
Herbst E., Cuppen H.M., 2006, Proceedings of the National Academy of 
Science, 103, 12257

[Hocuk & Cazaux(2015)]HC15
Hocuk S., Cazaux S., 2015, A&A, 576, A49

[Hollenbach et al.(2009)]Holl09
Hollenbach D., Kaufman M.J., Bergin E.A., Melnick G.J., 2009, ApJ, 690, 1497

[Ioppolo et al.(2008)]Iop08
Ioppolo S., Cuppen H.M., Romanzin C., van Dishoeck E.F., Linnartz H.,
2008, ApJ, 686, 1474

[Kalvns & Kalnin(2020)]KK20
Kalvns J., Kalnin J.R., 2020, A&A, 633, A97

[Keto & Caselli(2008)]KC08
Keto E., Caselli P., 2008, ApJ, 638, 238 [KC08]

[Keto & Caselli(2010)]KC10
Keto E., Caselli P., 2010, MNRAS, 402, 1625 [KC10]

[Keto, Rawlings & Caselli(2014)]KRC14
Keto E., Rawlings J.M.C., Caselli P., 2014, MNRAS, 440, 2616 [KRC14]

[Keto, Caselli & Rawlings(2015)]KCR15
Keto E., Caselli P., Rawlings J.M.C., 2015, MNRAS, 446, 3731 [KCR15]

[Keto & Field(2005)]KF05
Keto E., Field G., 2005, ApJ, 635, 1151








[Mathis, Rumpl & Nordsieck(1977)]MRN77
Mathis J.S., Rumpl W., Nordsieck K.H., 1977, ApJ, 217, 425

[McClure et al.(2023)]McC23
McClure M.K., Rocha W.R.M., Pontoppidan K.M., et al., 
2023, Nature Astronomy, 7, 431

[McElroy et al.(2013)]udfa12 
McElroy D., Walsh C., Markwick A.J.,
Cordiner M.A., Smith K., Millar T.J., 2013, A&A, 550, A36

[Millar et al.(2024)]udfa22
Millar T.J., Walsh C., Van de Sande M., Markwick A.J., 2024, A&A, in press

[Minissale et al.(2016)]Min16
Minissale M., Dulieu F., Cazaux S., Hocuk S., 2016, A&A, 585, A24





[Nelson & Langer(1997)]NL97
Nelson R.P., Langer W.D., 1997, ApJ, 482, 796

[Öberg et al.(2009)]OL09 
Öberg K., Linnartz H., Visser R., van Dishoeck E.F., 2009, ApJ, 693, 1209

[Pantaleone et al.(2021)]Pant21
Pantaleone S., Enrique-Romero J., Ceccarelli C., Ferrero S., Balucani N., 
Rimila A., Ugliengo P., 2021, ApJ, 917, 49

[Pontoppidan (2006)]Pont06
Pontoppidan K.M., 2006, A&A, 453, L47

[Rawlings et al.(1992)]RHMW92
Rawlings J.M.C., Hartquist T.W., Menten K., Williams D.A., 1992,
MNRAS, 255, 471





[Rawlings, Wilson & Williams(2019)]RWW19
Rawlings J.M.C, Wilson T.G., Williams D.A., 2019, MNRAS, 486, 10

[Rawlings & Williams(2021)]RW21
Rawlings J.M.C., Williams D.A., 2021, MNRAS, 500, 5117 

[Rawlings (2022)]Raw22 
Rawlings J.M.C., 2022, MNRAS, 517, 3804





[Redaelli et al.(2022)]Red22
Redaelli E., Chacón-Tanarro A., Caselli P., Tafalla M., Pineda J.E., 
Spezzano S., Sipilä O., 2022, ApJ, 941, 168

[Riedel et al.(2023)]Ried23
Riedel W., Sipilä O., Redaelli E., Caselli P., Vasyunin A.I., 
Dulieu F., Watanabe N., 2023, arXiv 2310.08389

[Röllig et al.(2007)]Retal07 
Röllig M., Abel N.P., Bell T., et al., 2007, A&A, 467, 187

[Rollins & Rawlings(2012)]RR12
Rollins R.P., Rawlings J.M.C., 2012, MNRAS, 427, 2328

[Silsbee et al.(2020)]Sil20
Silsbee K., Ivlev A.V., Sipilä O., Caselli P., Zhao B., 2020,
A&A, 641, 39

[Silsbee et al.(2021)]Sil21
Silsbee K., Caselli P., Ivlev A.V., 2021, MNRAS, 507, 6205

[Sipilä et al.(2019)]Sip19 Sipilä O., Caselli P., 
Redaelli E., Juvela M., Bizzocchi L., 2019, MNRAS, 487, 1269

[Sipilä, Silsbee & Caselli(2021)]SSC21
Sipilä O., Silsbee K., Caselli P., 2021, ApJ, 922, 126

[Sipilä et al.(2022)]SCRS22
Sipilä O., Caselli P., Redaelli E., Spezzano S., 2022, A&A, 668, A131

[Taquet, Ceccarelli & Kahane(2012)]TCK12
Taquet V., Ceccarelli C., Kahane C., 2012, A&A, 538, A42

[van Dishoeck, Herbst & Neufeld(2013)]DHN13
van Dishoeck E.F., Herbst E., Neufeld D.A., 2013, Chem. Rev., 113, 9043

[van Dishoeck et al.(2021)]vD21
van Dishoeck E.F., Kristensen L.E., Mottram J.C. et al., 2021, A&A, 648, 24

[Vasyunin et al.(2017)]Vas17 
Vasyunin A.I., Caselli P., Dulieu F., Jiménez-Serra I., 2017, ApJ, 842, 33

[Wakelam et al.(2012)]Wak12
Wakelam V., Herbst E., Loison J.-C., et al., 2012, ApJS, 199, 21

[Willacy et al.(1994)]WRW94 
Willacy K., Rawlings J.M.C., Williams D.A., 1994, MNRAS, 269, 921

[Woodall et al.(2007)]W07 Woodall J., Agúndez M., 
Markwick–Kemper A.J. & Millar T.J., 2007, A&A, 466, 1197

[Zucconi, Walmsley & Galli(2001)]ZWG01
Zucconi A., Walmsley C.M., Galli D., 2001, A&A, 376, 650









§ SIMPLIFIED CHEMICAL NETWORKS








KC08 suggested that the equilibrium abundances of CO, C, and C^+, and CO ice 
can be approximated by a simultaneous solution of the equilibrium rate equations, and 
that this approximation is adequate for quiescent molecular clouds such as starless cores. 
KRC14 suggested a similar procedure for the abundances of O, OH, H_2O and H_2O 
ice as discussed in detail in Appendix B of that paper. From the discussion above 
and in <ref>, we can now define a revised reduced chemical network 
(figure <ref>) for carbon and CO. This allows us to simultaneously solve the 
rate equations above for the equilibrium abundances of C, C^+, CO, and CO ice 
independently of STARCHEM. To do so, we assume that the abundance of OH is given 
by the reduced network chemistry for oxygen species, and the abundance of free 
electrons is equal to the abundance of C^+. The latter results in a 
quadratic equation for C^+. 
The carbon-chemistry reaction scheme and recent values for the rate 
coefficients <cit.> are:

    C + hν→ C^+ + e^-       3.5×10^-10 e^-3.76A_v s^-1
     C^+ + e^- → C + hν       1.23×10^-11T^-0.29 
    e^17.6/T cm^3s^-1
     C + OH → CO + H       1.0×10^-10 cm^3s^-1
     C^+ + OH → CO^+…→ CO       1.33×10^-8T^-0.5 cm^3s^-1
     CO + hν→ C + O       2.4×10^-10 e^-3.88A_v s^-1
     CO → CO_(ice)       freeze-out
     CO_(ice)→ CO        desorption

The rates for the freeze-out and desorption reactions are described in 
<ref> and <ref>, respectively.

The simple description of ice desorption considers the rate to be proportional to the 
total solid-state abundance of a species and omits the detailed calculations of the 
several grain surface processes in STARCHEM, including the desorption 
due to the energy released by the formation of molecules on grain surfaces. 
Owing to the quadratic equation for C^+, we solve the set of simultaneous rate
equations by elementary algebra rather than matrix inversion as was used in KRC14.



§ THE STARCHEM MODEL



Although STARCHEM has been
applied in previously published studies <cit.> it was originally
formulated for the study of pre-stellar and protostellar cores, so
we give a fuller description of the model here. 
The model was designed to allow fully flexible descriptions of the 
chemistry and dynamics and the main features of the generic code are as follows:



  * The chemistry and physical parameters are co-integrated using the 
DLSODE integration package, along Lagrangian streamlines, for a
number of test points in the chosen geometry. The model has 
been constructed so as to be easily adaptable to any user-defined 
dynamics.
The physical parameters (such as the density, gas and dust temperatures, 
visual extinction and velocity field) can either be specified analytically,
or in look-up tables, with the continuous values 
determined though smoothed interpolation. 

  * The chemistry is defined by flexible descriptions of 
gas-phase reactions, 
gas-grain interactions (in the surface layers of the dust ices), 
and solid-state (sub-surface) processes. 
In particular, the chemical datasets (species and reaction rate coefficients) 
are dynamically switchable - allowing the sensitivity to the adopted 
chemistry to be tested.

  * Where appropriate (i.e. at relatively low values of extinction), 
the code incorporates a model of the photon dominated region (PDR) to 
describe the photochemistry.

  * Freeze-out, plus a full suite of desorption processes are included
(thermal, direct and CR-induced photodesorption, CR heating - spot and whole grain,
heating by H_2 formation and species-specific enthalpy of formation etc.), 
together with a simplified surface chemistry.

  * The chemical composition of the ices that form on the dust grains
are determined with time- and position-dependence on a 
layer-by-layer basis, in a so-called `three-phase' model, which differentiates
the chemically active surface layers and the chemically inert
sub-surface layers.


Details of these various features are given below, as well as the specific 
application to our model of L1544.



 §.§ The gas-phase chemistry



One of the features of STARCHEM, and partly to allay any concerns that the
results are artefacts of errors and inconsistencies in reaction rate data, is its
flexibility to be able to dynamically switch between different chemistries. 
We consider three databases:
(a) UDFA 2012 <cit.>,
(b) UDFA 2006 <cit.> and
(c) KIDA 2012 <cit.>.
As our basis set, we use the first of these. 
In this study, the reaction set consists of some 1250 chemical reactions 
between 98 species, involving the elements H, He, C, N, O, S and a representative
low ionization potential metal, Na. Molecules are limited to species composed 
of less than 6 atoms.
The chemistry includes ion-molecule, neutral-neutral, cosmic ray ionizations and 
photochemical reactions as applicable to molecular clouds. Secondary 
(cosmic ray induced) photodissociation is also included.

In this study we made a few alterations to the reaction data.
e.g. the rates given in UDFA Rate12 for the reaction 

    H + OH → H_2O

are probably anomalous (T.J. Millar, private communication):

    k = 5.28× 10^-18 (T/300)^-5.22exp (-90/T).

If this were correct, then due to the strong 
temperature-dependence, this reaction would be a dominant contributor to 
the H_2O formation rate (in the low temperature, gas-phase limit).
Instead, we adopt the values given in the Ohio State University (OSU) database:

    k = 4.0× 10^-18 (T/300)^-2.0




 §.§ The Photon Dominated Region



Photochemistry clearly plays an important role in the chemical evolution
of starless cores, both through gas-phase photolysis and the desorption
mechanisms described below.
In the outer parts of the core (r 0.07pc) the mutual shielding of the 
H_2, C and CO photodissociation and photoionization continua, coupled to the
dust absorption are very important.
The model includes a simulation for the photon-transmission probability
in a photon dominated region (PDR). This is particularly important for 
the determination of the rates for H_2 photodissociation, C photoionization 
and CO photodissociation, and hence the H/H_2 and C^+/C/CO transitions.
This is a simplified approach that is based on accurate analytical fits 
(as functions of extinction and radiation field strength) to the data 
presented in Figures 8, 10 and 13 of <cit.>.








The code also includes allowances for the mutual shielding effects.
Specifically, the shielding by the CO photodissociation and C photoionization 
lines/continuum has significant implications for the photochemistry of 
other species; most notably N_2 and CN <cit.>. 



 §.§ Freeze-out



Following <cit.>, the rate of change of the gas-phase fractional abundance 
of a species, i, due to freeze-out onto the surface of grains is given by:

    Ẋ_ i = -( k_ B/2π m_ H)^1/2σ_ HS_ iC( T/m_ i)^1/2n_ HX_ i

where: 
σ_ H is the mean surface area of dust grains per hydrogen nucleon
and is the grain population-averaged value of 
4π⟨ a^2⟩ .d_ g.
Here d_ g is the dust to gas ratio (by number, relative to n_ H) and
the factor of 4 in this expression implicitly assumes that the grains are
spherical.
n_ H is the hydrogen nucleon density (cm^-3),
S_ i is the sticking coefficient,
T is the gas kinetic temperature,
m_ i is the molecular mass of species i in amu,
and C is a factor which takes into account electrostatic effects for ions 
(assuming an average grain charge of -1); C=1 for neutral species, whilst for singly 
charged positive ions C = 1+ (16.71/aT).

Following laboratory measurements and quantum mechanical calculations we assume 
that S_ i∼1.0 for all species except that the possibility for the
differentiation of ions and neutrals is included.
It is also assumed that electrons impinge upon and stick to grains at a rate 
that is just sufficient to neutralise all accreted ionic species.



 §.§ Desorption



Comprehensive descriptions of the various continuous desorption mechanisms 
are included in STARCHEM. These are:
(i) thermal desorption,
(ii) direct photodesorption by the (attenuated) interstellar radiation field,
(iii) photodesorption driven by radiation field generated by the cosmic-ray
excitation of excited states of H_2,
(iv) desorption due to cosmic ray heating of whole grains or ice mantle hot spots,
(v) non-selective chemical desorption driven by the enthalpy of H_2 
formation on grains, and
(vi) species-specific chemical desorption driven by the enthalpy of their 
formation (e.g. for H_2O, NH_3 etc.).

Most desorption processes involve sublimation from the 
surface layers of the ices that depends on the
compositional structure of the ice.
An ice mantle that is predominantly composed of one molecular species 
(e.g. H_2O) may be covered by ice layers of a different composition (e.g. CO) 
that effectively protects the lower levels of the ice against desorption. 

We therefore allow for the detailed time-dependent
variations of the composition in the individual ice layers, 
and adopt the so-called `three-phase' model of the chemistry; 
differentiating between the gas-phase, the `active' surface ice layers and the 
inert ice core.
This implementation is described in more detail in subsection <ref>.



  §.§.§ Thermal desorption



The rate of thermal sublimation (per molecule) for a zeroth order process 
is given by:

    k = ν_0 e^-E_ bind/k_ BT_ dust s^-1

where E_ bind is the binding energy of the adsorbate, 
T_ dust is the dust temperature and ν_0 is the vibration 
frequency of the adsorbed molecule given by 

    ν_0 = √(2N_sE_bind/π^2m_i)

where N_ s is the number of binding sites per unit area, and m_ i 
is the mass of the adsorbed particle. Typically, ν_0 = 10^12-10^13s^-1.
This implies an evaporation rate  of 

    f_ iN_ sσ_ H n_ Hν_0 e^-E_bind/kT_dust,
 
where f_ i is the fractional surface coverage of species i and 
σ_ H is the grain surface area per hydrogen nucleon, as described above.

The exponential dependence on the dust temperature implies a critical temperature,
primarly defined by the binding energy of the adsorbate. This is rather ill-defined 
in ices of mixed composition although, for the dust temperature profiles given in 
Figure <ref>, thermal desorption is not expected to be important
except for the most volatile species.



  §.§.§ Photodesorption



The photodesorption of molecules from ices proceeds by
absorption of a far-UV photon (which typically has an energy greater than the
bond energy of the absorbing molecule, and very much greater than the surface 
binding energy of the molecule). This results either in the direct
desorption of the molecule, or else its photodissociation 
in situ. The dissociation may be followed by the desorption of one or 
more of the dissociation products. Alternatively, the molecule may recombine 
and be chemically desorbed as a result of the energy released in the bond 
formation. 


For H_2O the process is driven by solid-state photodissociation. 
This has been studied, both experimentally and theoretically, in some 
detail <cit.> 
and in our model we include two possible photodesorption channels:

    H_2O(s) + hν→ H_2O


    H_2O(s) + hν→ OH + H

although the channel

    H_2O(s) + hν→ O + H_2

may also be possible with lower efficiency.
The laboratory <cit.>, 
theoretical <cit.> 
and observational 
evidence <cit.> all suggest that the yield for dissociative photodesorption is approximately
twice that for non-dissociative desorption, and we split the total yield 
for photodesorption between these processes accordingly.




Similar branchings are included for other important ice mantle components, such as 
methane:

    CH_4(s.) + hν→ CH_2 + H_2
 

    CH_4(s.) + hν→ CH_4
 

The UV photons that cause the desorption can only penetrate 
the upper (maybe 3 or 4) monolayers of ices - with quite limited efficiency.
In general, the desorption flux (cm^-2) due to photodesorption
is given by:

    F = Y_ pd,iF_ UVf_ i

where Y_ pd,i is the photodesorption yield for species i, 
f_ i is the fraction of binding sites in the surface layer(s) that are 
occupied by species i (as determined from the model of the compositional 
structure of the ice layers) and F_ UV is the photon flux, given by:


    F_UV = G_0F_0e^-β A_ v + F_cr(ζ/ζ_0),


where ζ_0=1.3-17s^-1 is the standard value for the 
cosmic ray ionization rate.
In this expression the first term accounts for direct photodesorption; 
F_0∼ 10^8photons cm^-2 is the unshielded interstellar UV flux, with a
scaling factor G_0 and the exponential accounts for shielding by dust. 
F_ cr is the cosmic ray induced photon flux <cit.> and 
accounts for cosmic-ray induced photodesorption. 

There are therefore three major uncertainties in quantifying the 
photodesorption rates:
(i) the (total) yield for the process Y_ pd,i, and its 
dependence on the substrate where the molecule is located,
(ii) the products and branching ratios of the desorption process (if 
dissociative), and
(iii) the A_ v-dependence of the rate (β).

In the case of H_2O all three of these have been well studied.
Y_ pd has been measured in the laboratory to be 
∼ 10^-4-10^-3 for pure water ices on a flat gold substrate. 
It seems unlikely that the same results would apply to 
impure, mixed ices on amorphous and morphologically complex grain surfaces.













The photodesorption yields for other species, and particularly for ices of 
mixed composition, are poorly understood.
However, it is unlikely that the photodesorption of 
H_2O molecules would occur whilst leaving other (perhaps more volatile) species 
bound to the ices, particularly in ices that are predominantly composed of H_2O.
So, we make the simplifying assumption that the same photodesorption yield 
applies to all ice species, although we consider that yield to be a
free parameter. 

The values of β are also highly uncertain and depend on the 
(a) the photodesorption efficiency, and (b) the dust opacity, which are both
wavelength dependent.
In the absence of information, we crudely assume (in common with many other 
studies) that the same value of β=1.8 for H_2O applies to all other 
species.
We recognise that this is not well justified and could result
in significant errors in the gas-phase abundances determined in regions 
of moderate extinction (A_ v∼ 1-16). 



  §.§.§ Desorption driven by cosmic-ray heating



The interaction between cosmic rays and ice-coated grains is
complex and ill-defined. Cosmic ray strikes may lead to desorption 
in a number of ways including direct sputtering, heat-spike 
pressure-pulse desorption and via grain heating.

At one extreme, high energy cosmic rays may be capable of 
totally disrupting (small) grains, whilst at the other they 
may simply act to heat the grains.
In this study we adopt the commonly used grain heating 
model - in which cosmic ray impacts, that typically occurring once every 
Myr, result in grain heating and the rapid desorption of material from 
the surface layers, accompanied by efficient evaporative cooling. 

This is a sporadic thermal desorption process
that may result either from whole grain heating or, in some cases, 
spot heating at the entry/exit points of the cosmic ray on the 
grain surface.
The efficiency of the process depends on the binding energies of the 
molecules and the spectrum of the cosmic rays, together with the duty 
cycle of grain heating and cooling.

Recent studies have questioned the various assumptions that go into
these models. These indicate that the previously adopted values of the 
effective desorption rates may have been significantly underestimated
<cit.> and that the sporadic nature of the 
process cannot easily be reduced to continuous desorption rates 
<cit.>.

For the ease of comparison, and to test our stated
objectives, we have continued to use the updated, continuous desorption 
rates determined for whole grain heating. We add a note of 
caution concerning their validity.

With this assumption, the rates can be expressed <cit.> as:

    ṅ=R_ crn_ i

where R_ cr is the desorption rate per molecule and n_ i is
the volume abundance of solid-state species i.

There are several variations in the calculated values for the 
desorption rates. For example,
by considering more localised ice mantle `hot spot' heating, <cit.>
calculated a rate for H_2O desorption that is approximately 30 times
larger than the value given in <cit.>. 
We use these `hot spot' values for H_2O, and other species, in our model. 



  §.§.§ Desorption driven by H_2 formation



The formation of H_2 molecules on the surface of grains yields an energy
equal to the bond energy of H_2. A significant fraction of this energy may be 
deposited locally in the substrate and conducted to the rest of the grain.
This energy may be sufficient to desorb other nearby molecules, non-selectively 
<cit.>. 
As H_2 formation occurs on the surface of grains the
desorption process is limited to molecules in the surface layers.
In this case, the rate of desorption (cm^-3s^-1) of a species i is
given by

    ṅ_ i = Y_ H_2,i.ṅ_ H_2.f_ i

where Y_ H_2,i is the desorption yield for species i, 
f_ i is the surface coverage of species i, as defined above, and 
ṅ_ H_2 is the rate of (surface) 
H_2 formation, which we determine from the hydrogen atom
impingement rate on grains.



  §.§.§ Desorption driven by the enthalpy of formation



We also include the selective desorption of individual molecules, 
driven by the enthalpy of their formation on the surfaces of 
bare or ice-coated grains. This is important for several species, such as H_2CO
and - most especially - OH and H_2O. To quantify this, we use the results of
<cit.> who determined the chemical desorption efficiencies to
be ∼ 50%, 
amorphous water ice substrates respectively, and ∼ 50%, ∼80% and 
∼ 30% for the OH+H→H_2O reaction on oxidised graphite, amorphous 
silicates and amorphous water ice, respectively.
This dependence of the desorption efficiencies on the nature of the surface
is included in the model. 
Additional complexities, such as the likelihood of higher desorption efficiencies
from CO-rich ices <cit.> have not yet been included in the model.



 §.§ Surface chemistry



In the current version of STARCHEM we do not consider the thermodynamics 
and kinetics of surface chemical processes, which are dominated by the thermal 
hopping of hydrogen and other, light, reactive species. Instead we apply a limited
and very simple realisation of the surface chemistry, and consider the fractional 
conversion of species, subject to empirical constraints.

Thus, all neutral molecules and radicals are allowed to stick to grains and,
where appropriate, hydrogenate to saturation 
(e.g. C, C^+, CH, CH_2, CH_3 → CH_4 and N, N^+, NH, NH_2 
→ NH_3 etc.).
Some atomic ions (He^+, Na^+ and S^+) are allowed to recombine 
and return to the gas phase, whilst all molecular ions are assumed to 
dissociatively recombine, with the products returned to the gas-phase. 

H_2 formation is described using the same formalism as for other species.
For the temperature-dependence of the combined sticking and reaction probability 
for the conversion of H to H_2 we use the simple formula of <cit.>.

For the surface chemistry that results in the 
formation of H_2O and CO_2 we note the following:


(a) H_2O and CO_2 are not very reactive and so can be treated as being 
chemically inert on dust grains,



(b) <cit.> suggest that CO_2 
is formed on a CO ice surface via the surface reaction of O or OH with CO, 
and  

(c) empirically, for Taurus cores, H_2O is the dominant (observed) ice,
with the CO and CO_2 abundances relative to H_2O being ∼10-40% 
and ∼20% respectively <cit.>.


For the oxygen chemistry, the freeze-out and surface reactions of 
significance are (representing solid-state species by `(s)'):


    O/O^+ + grain → O(s)


    OH + grain → OH(s)


    O, O^+, OH + CO(s) → CO_2(s)


    O(s) + H → OH(s)/OH


    OH(s) + H → H_2O(s)/H_2O


    O_2 + grain → O_2(s)


    O_2(s) + H  ....→ H_2O(s)/H_2O

The chemical endpoint for the oxygen surface chemistry is thus dominated 
by water formation.

For the very low dust temperatures within starless cores, CO is much less 
mobile than hydrogen atoms on grains surfaces <cit.>. So, to allow 
for the formation of CO_2(s) (as opposed to hydrogenation products)
we assume that CO_2(s) is formed by the non-diffusive prompt 
(Eley-Rideal) reaction of gas-phase O, O^+ or OH with bound CO molecules (reaction <ref>).
The impact probability is thus proportional to the surface coverage of CO, 
f_ CO, and the reaction efficiency is f_ CO.F_ CO_2, where 
F_ CO_2 is the reaction probability per impact.

In typical simple models of gas-grain interactions F_ CO_2=0.1 is
typically adopted to account for the observed CO_2:H_2O abundance ratio. 

For the O, O^+ or OH that does not react with CO, and providing the timescale 
for the surface residence of adsorbed oxygen atoms is greater than the hydrogen 
atom migration timescale over the surface, we assume that the hydrogenation reactions
(<ref> & <ref>) are efficient, so that a fraction
(F_ reac∼ 1) is hydrogenated to form H_2O.

A proportion of the surface-formed OH and H_2O may be desorbed due to the 
enthalpy of formation (D_ OH and D_ H_2O, respectively). 
The efficiencies of these processes are discussed and quantified
in <cit.> and <cit.>, as discussed above.
Following that work, we can include formation enthalpy-driven 
desorption with D_ OH=0.5,0.25 and D_ H_2O=0.8,0.3 for formation on
bare, refractory, grains and amorphous H_2O ice, respectively.

For O_2 we assume that a fraction (F_ O_2, typically ∼1) is
ultimately hydrogenated to H_2O (reaction <ref>) via several possible 
chemical pathways, that we do not discuss and, again, a proportion (D_ H_2O) 
of this H_2O is desorbed.
The remaining O, O^+, OH and O_2 not undergoing these reactions is simply assumed 
to freeze out (reactions <ref>, <ref> & <ref>).

The binding energies (E _B/k) for O, OH, O_2 and H_2O are
taken to be 800K, 1300K, 1210K and 5770K, respectively,
although some authors claim that the binding energy of oxygen 
atoms may be as high as 1660 or 1850K on amorphous ice and bare refractory 
grains, respectively <cit.>.
Our approach allows for the various poorly-defined competing processes of 
freeze-out, desorption, hydrogenation and more complex reactions. 

For the carbon chemistry, other than CO, the dominant species in 
the freeze-out/hydrogenation reactions will be C/C^+. The intermediate 
species CH has a very low binding energy (∼650K) and so it is likely 
that enthalpy-of-formation driven desorption may be important:

    C,C^+ + grain → CH

We treat the efficiency of this process as a free parameter (f_ CH).

We should also include the hydrogenation of CO, leading to the formation of 
methanol (CH_3OH). This appears to be a very efficient mechanism on the 
surface of dust grains at low temperatures <cit.>.
We do not, however, follow the chemistry of CH_3OH in this study.



 §.§ Ice composition and layering



In STARCHEM we pay careful attention to the composition of the ices 
and trace this composition on a microscopic, layer-by-layer basis.
We recognise that both the freeze-out of gas-phase species and most 
desorption processes result in the addition, or removal, of the top 
layer(s) of the ice leaving the sub-surface layers unaffected. 

This implies that, at any given time, the nature of the gas-grain interactions
is determined by the structure and composition of only the surface layers.
This is the only chemically active part of the ice, within which surface
migration and reaction can occur; the sub-surface layers are
treated as being chemically inert and, in most situations, chemically 
inaccessible. We do not allow for any diffusion between the layers.
This is similar to the so-called `three phase' approach adopted by in other
studies <cit.>.
In regions where there is progressive freeze-out the sub-surface layers will
be a `fossil record' of the prior chemical/dynamical evolution of
the cloud. Unless the ices are disrupted or sublimated, the subsurface layers
are removed from the active chemistry and the net budget of accessible volatiles.
This implies that, even in those situations where chemical
quasi-equilibrium exists, the observed gas-phase abundances will be dependent on 
the prior evolution of the cloud.
In addition, the composition of the surface layers may be quite different to 
that of the bulk ice. In these circumstances, the effect of the outer layers
`shielding' the inner layers against surface desorption results in significant
imbalances between the 
gas-phase and solid-state abundances of species. This is an idea that has been 
explored in the context of the anomalous O_2 abundances seen in comet 67/P 
<cit.>.
Moreover, the compositional differentiation may also lead to significant errors 
in the estimation of solid-state abundances from observational data.

The determination of the composition of the ice mantles in environments
where ice layers can be added or removed is a non-trivial computational 
exercise but is one that is essential if we want to correctly quantify 
the gas-grain interactions.

The chemically significant parameter at 
any one instant is the number of surface binding sites per grain, which is 
given by

    4πa^2N_ s

where N_ s is the surface density of binding sites (cm^-2)
or, expressed as a `fractional abundance of binding sites' (i.e. the number of
binding sites per hydrogen nucleon):

    f_ sites = σ_ HN_ s.

This parameter is, of course, not fixed and will vary from layer to layer,
as the radius of the grain+ice, and hence the surface area per hydrogen nucleon 
(σ_ H) will vary as ice layers are deposited or removed.
The fractional coverage of the surface layers by species i is then given by

    f_ i = X_ i,s/f_ sites

where X_ i,s represents the fractional abundance of species i, apportioned
to the surface layers.

Although this is similar to some other studies <cit.>
we do not include some of the additional complexities of those models
(e.g. the inclusion of grain porosity/lattice defects etc., as appropriate for 
`fluffy aggregate' grain models) as there are a 
number of assumptions and simplifications in (all of) the ice layering models 
that are probably more significant.
These include:


(i) the assumption that dust grains are spherical,

(ii) the assumption that ices form in neat, sequential, concentric, layers. 
It is possible that the ice layers may not stack up neatly and there may be
some parts of a dust grain that have many more layers of ice than other parts,

(iii) whilst the segregation of layers implies that each layer
will have its own chemical composition, there is no allowance for any compositional
variations within each ice layer. We might not expect polar and 
apolar components to be fully mixed and there may be `patches' of ices with different
composition and binding properties as a result, affecting both the freeze-out and 
desorption properties of the surface layers, and

(iv) the assumption that the ice layer structure remains unaltered, other than as a
result of freeze-out or desorption. This could break down, for example, in the case
of energetic cosmic ray heating which acts to melt/partly sublimate the ices and 
hence significantly disrupt the integrity of the ice layers. 





We treat each ice layer as having a (defined) constant thickness equal to the 
mean bond length in the ice. This allows us to calculate 
how the grain radius including the ice mantle (and hence the 
grain surface area) varies with the number of accreted layers. As the freeze-out 
(and surface desorption) rates are proportional to the surface area, we can see how
- in a region where freeze-out dominates over desorption - the rate of freeze-out
can rapidly increase (the `snowball effect').
The accretion of ices may have other effects whose
implications we defer to a future study. These include, for example, the effective
removal of the smallest grains from the grain size distribution 
<cit.> and variations in the refractive index and albedo of the grains.

Using the method described above, we can determine the fractions 
of the binding sites in the upper layers that are occupied by the different species, 
with the possible inclusion of additional factors, such as allowances for
variations in the binding energy of key adsorbates as a function of surface 
layer composition.
Ultimately, this allows us to calculate the surface desorption rates and the 
efficiencies of the surface reactions, such as those that result in the formation 
of H_2, H_2O and CO_2.







