Why does ammonia not freeze out in the center of pre-stellar cores?

TL;DR

This study investigates why ammonia does not freeze out in the center of pre-stellar core L1544, finding that current models cannot reproduce the observed increase and suggesting chemical desorption efficiency constraints.

Contribution

The paper provides a detailed parameter-space exploration of ammonia abundance in L1544, highlighting limitations of existing models and proposing constraints on chemical desorption efficiency.

Findings

Models fail to reproduce observed ammonia increase in core center.

Chemical desorption efficiency likely below 1%.

Time-dependent effects are crucial for accurate modeling.

Abstract

We carried out a parameter-space exploration of the ammonia abundance in the pre-stellar core L1544, where it has been observed to increase toward the center of the core with no signs of freeze-out onto grain surfaces. We considered static and dynamical physical models coupled with elaborate chemical and radiative transfer calculations, and explored the effects of varying model parameters on the (ortho+para) ammonia abundance profile. None of our models are able to reproduce the inward-increasing tendency in the observed profile; ammonia depletion always occurs in the center of the core. In particular, our study shows that including the chemical desorption process, where exothermic association reactions on the grain surface can result in the immediate desorption of the product molecule, leads to ammonia abundances that are over an order of magnitude above the observed level in the innermost 15000 au of the core - at least when one employs a constant efficiency for the chemical desorption process irrespective of the ice composition. Our results seemingly constrain the chemical desorption efficiency of ammonia on water ice to below 1%. It is increasingly evident that time-dependent effects must be considered so that the results of chemical models can be reconciled with observations.

Tables

Table 1. Table 1: Initial abundances (with respect to the total hydrogen number density n H subscript 𝑛 H n_{\rm H} ) used in the chemical modelling. The initial H 2 subscript H 2 \rm H_{2} ortho/para ratio is 1.0 × 10 − 3 1.0 superscript 10 3 1.0\times 10^{-3} .

Species	Abundance
$H_2$	$5.00\times {10}^{-1}$
$\mathrm{He}$	$9.00\times {10}^{-2}$
$\mathrm{HD}$	$1.60\times {10}^{-5}$
$C^{+}$	$1.20\times {10}^{-4}$
$N$	$7.60\times {10}^{-5}$
$O$	$2.56\times {10}^{-4}$
$S^{+}$	$8.00\times {10}^{-8}$
${\mathrm{Si}}^{+}$	$8.00\times {10}^{-9}$
${\mathrm{Na}}^{+}$	$2.00\times {10}^{-9}$
${\mathrm{Mg}}^{+}$	$7.00\times {10}^{-9}$
${\mathrm{Fe}}^{+}$	$3.00\times {10}^{-9}$
$P^{+}$	$2.00\times {10}^{-10}$
${\mathrm{Cl}}^{+}$	$1.00\times {10}^{-9}$

Table 2. Table 2: Parameter variations and the associated model denominations considered in this paper. Our standard model is highlighted in boldface.

Model grid (total of 108 models)
Parameter		Values considered
Binding energy of ${\mathrm{NH}}_3$ (K)	5500	3000	1000
Sticking coefficient of atomic N	1	0.3
Photodesorption (of CO, $\left(o/p\right)H_{2}O$, ${\mathrm{CO}}_2$, $N_2$)	Yes	No
Photodesorption of $\left(o/p\right){\mathrm{NH}}_3$	Yes	No
Chemical desorption	Yes	No
External $A_V$ (mag)	5	2	1
Single models (see text for further explanations)
S1: Low elemental N abundance
S2: Multilayer ice chemistry
S3: Initial abundances from a lower-density cloud model
S4: C19 density profile
S5: Gas-phase chemistry only
S6: Initial $H_2$ o/p ratio of 0.1
S7: CR ionization rate of ${\zeta }_p={10}^{-16}{s}^{-1}$
S8: As S7, but no chemical desorption
S9: Grain radius $a_g=0.2\mu m$
S10: Grain radius $a_g=0.05\mu m$
S11: Hydrodynamical model

Equations

$$k_{\mathrm{secphot}}(i)={\zeta }_p(\mathrm{H})X({\mathrm{H}}_2)\frac{p_i}{1-\omega }[{\mathrm{s}}^{-1}],$$

Full text

††pubyear: 2018††pagerange: Why does ammonia not freeze out in the center of pre-stellar cores?–A

1 Introduction

Ammonia is observed almost ubiquitously in the interstellar medium (ISM). It serves as a useful tool for measuring the kinetic gas temperature because of its particular spectroscopic properties, and hence understanding its chemical evolution allows us to deduce important information on physical processes in the ISM. The gas-phase chemistry of ammonia is well understood, but its evolution on the surfaces of interstellar dust grains is rather poorly constrained. Furthermore, the strength of ammonia desorption from the grain surfaces, and the nature of the desorption mechanism, are still open questions.

Chemical models of star-forming regions predict that, at low temperature ($T\sim 10\,\rm K$), ammonia freezes out onto grain surfaces already at medium densities of a few times $10^{5}\,\rm cm^{-3}$, which is attributed to its high binding energy (Aikawa et al., 2012; Taquet et al., 2014; Sipilä et al., 2015b; Hily-Blant et al., 2018). These results are not in agreement with observations. Ammonia depletion has been observed toward pre-stellar cores, but it seems to occur only at very high (column) densities (Friesen et al., 2009; Ruoskanen et al., 2011; Chitsazzadeh et al., 2014). To add to the conundrum, Crapsi et al. (2007) derived an ammonia abundance profile in L1544, a well-studied pre-stellar core in Taurus, that increases toward the dust peak and shows no signs of depletion in the center of the core despite the high gas density ($n(\rm H_{2})>10^{6}\,\rm cm^{-3}$).

The discrepancy between the modelling results and observations is puzzling given that ammonia is a relatively simple molecule and its main formation and destruction pathways consist of a small number of reactions. This problem must be investigated in detail so that our understanding of the gas-grain chemical interaction can be improved. To this end, we used a comprehensive gas-grain chemical model to simulate the abundance of gas-phase ammonia in L1544. We varied several modelling parameters that are a priori expected to influence the ammonia abundance, in an effort to produce solutions where ammonia depletion either does not occur, or is mitigated to the observed levels within the uncertainties.

The paper is organized as follows. Sect. 2 discusses our chemical code and recent updates to it. Here we also discuss the physical source models used in the paper, and our parameter-space approach to the modelling. In Sect. 3 we present our results which are discussed in Sect. 4. We give our conclusions in Sect. 5. Additionally, a benchmark of radiative codes is presented in Appendix A.

2 Model

2.1 Chemical model

We used an expanded version of the gas-grain chemical code described in Sipilä et al. (2013, 2015a, 2015b). In short, the code solves a system of rate equations connecting separate networks of gas-phase and grain-surface chemistry. The details of the fundamental chemical processes (e.g., adsorption, thermal and non-thermal desorption) including the relevant formulae are described in Sipilä et al. (2015a) and are not reproduced here for the sake of brevity. However, for the purposes of the present paper, it was necessary to expand the set of chemical processes considered in the code, as opposed to using the older model as described in Sipilä et al. (2015a). We list these additions below.

Cosmic-ray induced secondary photoreactions. The interiors of molecular clouds are well shielded from the ultraviolet (UV) photons prominent in the spectrum of the interstellar radiation field (ISRF). However, cosmic-ray(CR)-induced ionization of $\rm H_{2}$ followed by electron recombination can create an UV field of appreciable strength inside otherwise well-shielded regions (Sternberg et al., 1987; Gredel et al., 1989). This UV field may ionize and/or dissociate molecules in the gas phase and on the surfaces of grains. The rate coefficient for the ionization or dissociation of atom or molecule $i$ in this process is given by

[TABLE]

where $\zeta_{p}({\rm H})$ is the primary CR ionization rate per hydrogen atom, $X({\rm H_{2}})=n({\rm H_{2}})/n_{\rm H}$ is the fractional abundance of ${\rm H_{2}}$ ($n_{\rm H}$ is the total number density of hydrogen nuclei), $p_{i}$ is an efficiency factor for the ionization/dissociation reaction in question, and $\omega$ is the grain albedo (assumed $\omega=0.5$). We updated the list of reactions and efficiency factors included in the current release of the KIDA network (Wakelam et al. 2015; see below) using the data of Heays et al. (2017; their Table 20). The efficiency factors required here were obtained by simply dividing their photoionization/dissociation rates by $10^{-16}$ (see also Hily-Blant et al. 2018). As noted by Heays et al. (2017), the simple division of the rates leads only to an approximate agreement with the formalism of Gredel et al. (1989), but a more detailed treatment of this issue is beyond the scope of the present paper. We also note that the factor $X({\rm H_{2}})$ in Eq. (1) was originally missing in the work of Gredel et al. (1989), but is necessary for the present context as pointed out by Woodall et al. (2007) and Flower et al. (2007).

Chemical desorption. Two-body chemical reactions with one reaction product on the grain surface may lead to the immediate desorption of the reaction product, if the excess formation energy is absorbed by the grain (Williams, 1968; Watson & Salpeter, 1972a, b). This process is commonly referred to as chemical desorption, or reactive desorption. Different treatments of chemical desorption in the context of gas-grain models have been suggested in the literature. Here we adopt the approach of Garrod et al. (2007), in which exothermic surface reactions lead to desorption with a probability of $\sim$1%. Recent investigations (Dulieu et al., 2013; Minissale et al., 2016; Chuang et al., 2018) have shown that the efficiency of chemical desorption may vary significantly depending on the reaction and type of surface. These results have already been incorporated in chemical models (Vasyunin et al., 2017). However, because of the high uncertainties involved, we chose to follow the uniform $\sim$1% approach of Garrod et al. (2007).

$\rm H_{2}$ self-shielding. In molecular clouds, and in particular in starless and pre-stellar cores, hydrogen exists mostly in molecular form, and the $\rm H_{2}$ present in the outer cloud may efficiently shield the $\rm H_{2}$ in the inner cloud against UV radiation. This effect is important because of the role that H atoms play in surface chemistry. Also, $\rm H_{2}$ self-shielding affects the $\rm H_{2}$ ortho/para ratio (Sipilä & Caselli, 2018), which will in turn affect the chemistry of ammonia because the $\rm N^{+}+H_{2}\longrightarrow NH^{+}+H$ reaction that initiates ammonia formation is strongly endothermic in the presence of para-$\rm H_{2}$ (Dislaire et al., 2012). We adopted the $\rm H_{2}$ self-shielding factor from Draine & Bertoldi (1996). This factor depends on the total $\rm H_{2}$ column density (see Sect. 2.2 for details on the physical core model used here) and is included in the rate equations that involve the ionization or dissociation of (ortho or para) $\rm H_{2}$. In reality self-shielding applies to other abundant molecules such as CO (Visser et al., 2009) and $\rm N_{2}$ (Heays et al., 2014) as well, but the self-shielding factor is a function of column density which for these species is a highly time-dependent quantity111Unlike that of $\rm H_{2}$, which we assume to remain constant as we consider a static physical model; see Sect. 2.2.. Self-shielding is thus not included for species other than $\rm H_{2}$ because our model is not fully time-dependent (see Sect. 2.2). The limited inclusion of self-shielding does not affect the results presented in this paper as we concentrate on the inner, heavily shielded, areas of L1544.

Temperature-dependent sticking coefficients. One of the parameters controlling the adsorption of molecules onto dust grains is the sticking coefficient. Many chemical models assume that the sticking coefficients of the various species equal unity regardless of the temperature of the medium, and indeed this assumption was made in our previous models as well. For the present work, we updated the chemical model to include temperature-dependent sticking coefficients for selected species, namely H, $\rm H_{2}$, $\rm N_{2}$, CO, $\rm O_{2}$, $\rm CH_{4}$, and $\rm CO_{2}$, including deuterated variants whenever applicable. For atomic H, we use the parametrized expression from Cuppen et al. (2010). Their formula applies to a graphite surface, but the sticking coefficient on water ice is almost identical to that on graphite for temperatures below 100 K (Cuppen et al., 2010). We assume the same formula for atomic D. For the rest of the species listed above we adopted the sticking coefficients presented by He et al. (2016), who also derived a sticking coefficient for $\rm D_{2}$. For $\rm HD$ and deuterated methane we assume that the sticking coefficient equals that of $\rm H_{2}$ or $\rm CH_{4}$, respectively. Species for which there is no theoretical or experimental data are assumed to stick with an efficiency of unity.

The additions to the chemical model presented above are essential because they either directly or indirectly affect the abundance of ammonia, which is the main target of our simulations.

Our gas-phase reaction network is essentially the kida.uva.2014 network (Wakelam et al., 2015) that was deuterated and spin-state separated according to the prescriptions laid out in Sipilä et al. (2015a, b), and updated as described above. A similar procedure was applied to our base grain-surface network, which is an updated version of the one presented by Semenov et al. (2010). For reactions involving CRs or photons we assume the same rate coefficients in the gas phase and on the grain surface, with the important distinction that we only consider dissociation reactions on the grain surface (we assume that no ionic species exist in the ice). Photodesorption is also included for a limited set of species (Sipilä & Caselli, 2018). The secondary photoreactions discussed above were added to both the gas-phase and grain-surface networks, again ignoring pathways that produce ions on the grain surface. The gas-phase and grain-surface networks contain a combined total of $\sim$82000 reactions.

2.2 Physical model for L1544

By default, we use the density and (gas and dust) temperature structure from the one-dimensional (1D) L1544 source model published by Keto & Caselli (2010; hereafter K10; see also Keto et al. 2014). This model has been previously used for interpreting observations and for carrying out chemical modelling of L1544 in several studies, such as Bizzocchi et al. (2013), Sipilä et al. (2016b), Vasyunin et al. (2017). The density and temperature structures of the model core are plotted in Fig. 1.

However, recent observations by Chacón-Tanarro et al. (2019; hereafter C19) show that the central density of L1544 may in fact be a factor of $\sim$4 lower than that calculated by K10, and that the slope of the density profile outside the flat radius may be different. In C19 a new dust temperature profile for L1544 was also presented. These new density and dust temperature profiles are plotted in Fig. 1 along with the profiles of K10, and are used in one of the models introduced below (see Sect. 2.3). Recent dust continuum emission observations toward L1544 obtained with the Atacama Large (sub)Millimeter Array agree with the profiles derived by C19, but local density enhancements up to $\sim 10^{7}\,\rm cm^{-3}$ are also present within the central 1400 au (Caselli et al., 2019).

For the purposes of the present work, a gas temperature profile is required as well. To this end, we took the density and dust temperatures from C19 and used our hydrodynamical code (Sipilä & Caselli, 2018), with infall/expansion velocity forced to zero, to calculate a (representative) gas temperature, plotted in Fig. 1 as the red dashed line. The profile corresponds to $t=10^{5}\,\rm yr$ of chemical evolution, which is a reasonable timescale since the physical model is static and the abundances of the cooling molecules hardly change after $t=10^{5}\,\rm yr$ in the inner part of the core, which is our main interest presently, owing to the high density. We only used this alternative L1544 physical structure in one of our model runs, the details of which are explained in Sect. 2.3.

As in our previous works (e.g., Sipilä et al., 2016b), we derived radius-and-time-dependent chemical abundance profiles in L1544 by dividing the physical model into concentric shells, calculating the chemical evolution separately in each shell, and combining the results obtained in the different shells at a given time step. We assume that the gas is initially atomic with the exception of $\rm H_{2}$ and HD; the adopted initial abundances are presented in Table 1. These initial abundances are used for all of the models presented in this paper except when otherwise noted (see below).

2.3 Parameter-space exploration

The chemistry of ammonia is sensitive to various model parameters. For example, disregarding chemical desorption will have a direct impact on the ammonia abundance; the binding energy of ammonia is very high, and without efficient non-thermal mechanisms it will, effectively, not desorb at low temperature (close to 10 K). To investigate the sensitivity of the ammonia abundance on the model assumptions, we constructed a grid of models in which we varied several key parameters that influence the ammonia abundance, particularly on the grain surfaces. The parameters and their variations are displayed in Table 2. Our fiducial model is highlighted in boldface.

The binding energy of ammonia on water ice is uncertain. Our fiducial value of 5500 K is taken from Collings et al. (2004). However, lower values of the order of 3000 K have been suggested (see Kamp et al., 2017, and references therein). We consider both values in our model grid. Furthermore, as a test, we also consider a very low (ad hoc) value of 1000 K to explore its effect on the ammonia abundance.

We take a fiducial value of unity for the sticking coefficient of atomic N, as is most often assumed in astrochemical models. The possibility of lower values has been suggested by Flower et al. (2006), and we take a value of 0.3 as an alternative to the canonical value of unity.

By default, we consider the photodesorption of CO, (o/p)$\rm H_{2}O$, $\rm CO_{2}$, and $\rm N_{2}$. The assumed photodesorption yields are, respectively, $2.7\times 10^{-3}$ (Öberg et al., 2009a), $1.0\times 10^{-3}$ (Öberg et al., 2009b), $1.2\times 10^{-3}$ (Öberg et al., 2009a), and $3.0\times 10^{-4}$ (Öberg et al., 2009a). The $\rm CO_{2}$ yield only approximates the complex expression derived by Öberg et al. (2009a) that depends on the ice thickness, which we do not track in this paper. We assume that both ortho and para $\rm H_{2}O$ are photodesorbed with the same yield. We also tested cases where $\rm NH_{3}$ photodesorption is included, with a yield of $1.0\times 10^{-3}$ (Martín-Doménech et al., 2018). This yield is assumed to apply to both ortho and para $\rm NH_{3}$, and we note that the experimental yield has been derived for pure $\rm NH_{3}$ ice, and not for $\rm NH_{3}$ on water or an ice mixture.

We assume that the core model is embedded in a larger molecular cloud; the parameter “external $A_{\rm V}$” expresses the attenuation (in the visual) by the parent cloud. All together, the grid consists of 108 parameter combinations.

We also ran eleven single models separately from the parameter grid in order to test some additional assumptions. These models, also collected in Table 2, were run as separate cases mainly because of calculational time constraints. It takes roughly three hours to complete one model run on a standard desktop computer (with parallelization), and every new parameter change in the grid would double the required computational time. The use of supercomputing resources for the present work is however not necessary as we are not searching for an exact fit to available observational data, and the main interest is in investigating the general trends caused by the parameter variations.

The single models S1 to S10 correspond to our fiducial model with the particular changes indicated in Table 2. In model S1, we assumed an N elemental abundance of $2.14\times 10^{-5}$ (Wakelam & Herbst 2008; their model EA1). Model S2 incorporates a multilayer (three-phase) ice model as detailed in Sipilä et al. (2016a)222Although this model naturally tracks the ice thickness, we still used the approximate $\rm CO_{2}$ photodesorption yield noted above.. For model S3, we first calculated a single-point chemical model with $n(\rm H_{2})=5\times 10^{3}\,\rm cm^{-3}$, $T_{\rm gas}=T_{\rm dust}=15\,\rm K$, $A_{\rm V}=1\,\rm mag$, and extracted the abundances from that model at $t=10^{5}\,\rm yr$ to use as initial abundances for the fiducial model. Model S4 uses the C19 density profile for L1544 (see Sect. 2.2) instead of the Keto et al. profile. Model S5 considers gas-phase chemistry only, where the formation of $\rm H_{2}$, HD, and $\rm D_{2}$ is parametrized as in Kong et al. (2015). In model S6 we adopt an initial $\rm H_{2}$ ortho/para ratio of 0.1. Models S7 and S8 incorporate an increase of the CR ionization rate over our standard value of $\zeta_{p}=1.3\times 10^{-17}\,\rm s^{-1}$, and the latter model also excludes chemical desorption. In models S9 and S10 we multiply or divide our standard grain radius (0.1 $\mu$m) by a factor of two.

Finally, model S11 represents a hydrodynamical model calculation, the details of which are given in Sect. 3.3.

3 Results

We extracted the abundance profiles of ortho and para $\rm NH_{3}$ at different time steps from the models detailed in Table 2, and compared the results against the observations of Crapsi et al. (2007). We present the results of the comparison in what follows, broken down into the results from our parameter grid and the single models S1 to S11.

3.1 Parameter grid

Figure 2 displays the abundance profiles of ortho+para $\rm NH_{3}$ obtained from all of the models comprising our parameter grid, convolved to a beam size of 4$\arcsec$ for comparison with the observations by Crapsi et al. (2007). The abundance profile from Crapsi et al. (2007) is also shown for comparison. One property of the calculated models is immediately evident: we always obtain solutions where the ammonia abundance depletes at the center of the core, and the monotonically inward-increasing profile from Crapsi et al. is never reproduced.

Many of the model solutions are virtually identical to each other and overlap. Also, the solutions appear in distinct groups. It is therefore sensible to separate the solutions based on the individual parameter values. In Fig. 3, we highlight the effect of each parameter. The effect of the $\rm NH_{3}$ binding energy is shown in the top left panel. It is immediately evident that the very low (ad hoc) value 1000 K leads only to solutions where the modeled abundance is orders of magnitude above the observed one in the central areas of the core. These solutions can therefore be discarded, and are not displayed in any of the figures from here on. Notably, using a value of 3000 K or 5000 K has very little influence to the results (the red and blue lines overlap). The two remaining sets of curves correspond to models where chemical desorption is turned on or off, respectively (see the lower middle panel of Fig. 3), and these sets are further subdivided as described below.

The sticking coefficient of atomic nitrogen does not influence the solutions in a clearly predictable manner, and solutions with both high and low overall $\rm NH_{3}$ abundances can be obtained with both of the tested values of the coefficient. Photodesorption, with or without $\rm NH_{3}$, has a negligible influence on our results. We note that Furuya & Persson (2018) have recently demonstrated that ammonia photodesorption can have an impact on the gas-phase N and $\rm N_{2}$ abundances at low visual extinctions. Our model does not reproduce this effect. We performed test calculations which indicate that the effect of photodesorption is much greater in three-phase (i.e., the ice is separated into a mantle and a bulk) models than it is in the two-phase model that we consider here (except in model S2). Also, our description of ice chemistry is different from that of Furuya & Persson (2018), which causes discrepancy in the modeling results (priv. comm. with K. Furuya). A detailed comparison is out of the scope of the present paper.

Chemical desorption influences the results greatly; we obtain solutions with very high peak ammonia abundances if chemical desorption is included, while the peak ammonia abundance never rises above a few $\times~{}10^{-9}$ if chemical desorption is turned off. This is a very strong result given that we are using the conservative value of $\sim$1% for the efficiency of the chemical desorption process (Garrod et al., 2007), which is much lower than the values derived for some reactions by Minissale et al. (2016), for example. Finally, the magnitude of external $A_{\rm V}$ changes the shape of the ammonia abundance profile; if external $A_{\rm V}$ is low, we obtain rather narrow abundance distributions, while increasing values of the external $A_{\rm V}$ yield increasingly extended distributions.

In conclusion to the above: our models do not reproduce the shape of the observed ammonia abundance profile, and we always obtain solutions where $\rm NH_{3}$ depletes near the core center. If we disregard the discrepancy in the central few thousand au, the best fit to the observations of Crapsi et al. (2007) is reached with a model where chemical desorption is excluded, the sticking coefficient of atomic N is unity, and the external $A_{\rm V}$ is higher than 2 mag.

3.2 Single models S1 to S10

Figure 4 shows the $\rm NH_{3}$ abundance profiles predicted by the single models S1 to S10. A significant spread is evident in the peak ammonia abundances depending on the model and, most notably, none of the single models provide a solution where ammonia depletion does not occur.

Model S1 presents an ammonia abundance profile that is very similar to that given by our fiducial model, except scaled down. In model S2, where a three-phase ice description is adopted, ammonia depletes very strongly because of trapping in the inert ice bulk beneath the active surface layer(s) (for more details see Sipilä et al., 2016a). Furthermore, because of the trapping, the chemical desorption process is not efficient and we obtain a lower ammonia peak abundance than in the fiducial model. If the gas is first let to evolve in a diffuse cloud environment (model S3), the end result is very close to the fiducial model, suggesting that the initial conditions no longer play a role if the gas is let to evolve for a sufficient amount of time in the pre-stellar phase.

When the density profile from C19 is used (model S4), the ammonia depletion zone is larger than in our fiducial model. This is because the central high-density area, where ammonia depletes efficiently, is broader in the C19 model than in that of Keto et al., even though the density at the very center of the core is a factor of $\sim$4 lower in the former (see Fig. 1). The results from model S4 and the fiducial model are again qualitatively similar despite the rather large difference in the density profiles.

The ammonia abundance profile does not display an inward-increasing trend even in the gas-phase model S5. This effect is tied to the electron fraction. The production chain of ammonia starts with $\rm N^{+}$, which is produced in reactions between $\rm He^{+}$ and $\rm N_{2}$. At high density and low temperature, the rate of electron impacts with $\rm He^{+}$ is high, which inhibits the production of $\rm N^{+}$ (in contrast with gas-grain models where the abundance of $\rm He^{+}$ tends to increase with freeze-out). On the other hand $\rm N_{2}$ is not produced efficiently at low density, so that the formation of ammonia is (in a gas-phase model) the most efficient at a density of $\sim 10^{5}\,\rm cm^{-3}$. We note that previous gas-phase models have predicted higher ammonia abundances of the order of $10^{-8}$ even at high density (e.g., Le Gal et al., 2014; Roueff et al., 2015). These models however incorporated rate coefficients for some important reactions in the ammonia formation network that are much higher than the up-to-date values included in kida.uva.2014 and hence in the present model. The rate coefficient revisions lead to generally lower ammonia abundances, highlighting the great sensitivity of chemical models to uncertainties in input data.

In model S6 we tested a higher initial $\rm H_{2}$ ortho/para ratio. Evidently, we obtain an abundance profile that is very similar to the fiducial model, providing further evidence of the insensitivity of the ammonia abundance to the initial conditions.

Models S7 and S8 explore the effect of the CR ionization rate. In addition to testing the effect of simply increasing the ionization rate (S7), we also tested a case where the chemical desorption process is additionally turned off (S8). Evidently, the enhancement in the CR ionization rate helps to maintain a higher abundance of ammonia in the central core, although depletion still occurs at the highest densities. However, the effect of CRs is so strong that it counteracts to some degree the exclusion of chemical desorption, and we obtain even in model S8 an ammonia abundance that is much higher than the observed one.

Finally, in models S9 and S10 we varied the grain radius. Increasing or decreasing the grain radius does not modify the abundance profile in a significant way as compared to the fiducial model, and the same trends as in the majority of our models are displayed here as well.

3.3 Single model S11: hydrodynamics

We have recently demonstrated that chemistry plays a very important role in determining the dynamics of the collapse of a star-forming cloud, and that using a static model for the physical structure of the core – as in the present paper so far – either overestimates or underestimates chemical abundances in a collapsing core, depending on both radius and time (Sipilä & Caselli, 2018). That study was partly motivated by the previous observations of ammonia toward L1544, and therefore it is logical to consider if including dynamics could solve the present problem of excessively strong ammonia depletion appearing in the models, or at least to alleviate it.

To this end, we re-ran the hydrodynamical simulation described in Sipilä & Caselli (2018), i.e., we started with an unstable Bonnor-Ebert sphere with central density $n_{\rm c}({\rm H_{2}})=2\times 10^{4}\,\rm cm^{-3}$, temperature $T=10\,\rm K$, and mass $M\sim 7.15\,M_{\sun}$, but used the updated chemical networks described above. In Sipilä & Caselli (2018), the termination condition of the hydrodynamical code was set to correspond to the time step when the infall flow becomes supersonic, and in that paper this occurred at $t=7.19\times 10^{5}\,\rm yr$. In the present case the termination of the code occurred at $t=1.20\times 10^{6}\,\rm yr$. There are three reasons for this difference. First, the changes to the chemical setup introduced here affect the chemistry at low density in particular, and this is reflected on the infall velocity profile and hence on the collapse timescale. Second, we have fixed a minor coding error in the expression that compares the infall velocity and sound speed to determine the termination time; the old version of the code calculated the sound speed inadvertently using the maximum of the gas temperature instead of its local value (the effect of this error is very small). Third, and most importantly, the cooling efficiency of HCN is clearly lower in the present paper than in Sipilä & Caselli (2018; see also below). In fact, we have determined through extensive testing that the high HCN cooling efficiency presented in Sipilä & Caselli (2018) is a numerical error, the exact cause of which we have however not been able to identify. It remains unknown why only HCN was affected while the cooling powers of the other coolants are similar in both our current and earlier model calculations.

L1544 may be somewhat more massive than the Bonnor-Ebert sphere used here (for example K10 employed a Bonnor-Ebert sphere with mass $M~{}\sim~{}10\,M_{\sun}$), but we chose to use the same initial core configuration as in Sipilä & Caselli (2018) in order to easily track down the causes for the differences in the physical and chemical evolution due to sources other than the initial physical model.

The results of the hydrodynamical simulation are shown in Figure 5. Unlike in the other abundance plots presented in this paper which concentrate on the inner 25000 au, we plot here the results up to 50000 au to facilitate easier comparison to Figs. 3 and 4 of Sipilä & Caselli (2018). First, it is strikingly evident that the HCN cooling efficiency is very low in the current model, while other cooling efficiencies are similar, as compared to Sipilä & Caselli (2018). In particular, we still obtain a very strong contribution from NO in dense gas at late times. Comparison of the solution at the final time step ($t=1.20\times 10^{6}\,\rm yr$) to the present fiducial (static) chemical model shows that ammonia is very strongly depleted even in the hydrodynamical model. The depletion factor is smaller than in the static model, but this is only because the central density is less than $10^{6}\,\rm cm^{-3}$ at the time of the termination of the calculation. If the calculation was continued beyond this point, ammonia would deplete as strongly as it does in the static case. The ammonia depletion zone is seemingly smaller in the static model, but this is only because the K10 physical model is more centrally concentrated (Fig. 1) than the hydrodynamical solution at the final time step. We note that the abrupt changes in some of the cooling functions near the origin are only transient radiative transfer artifacts that do not affect the overall evolution of the core.

3.4 Time dependence

The results presented above display a clear discrepancy with the observations of Crapsi et al. (2007). From our models we always obtain an abundance profile that decreases strongly towards the center of the core. However, one further crucial aspect that we have not explored thus far is time-dependence. We plot in Fig. 6 the results from our parameter grid at three different time steps.

The ammonia depletion timescale is very short in the central areas of the core because of the high density. Even at $t=10^{4}\,\rm yr$, which is an unrealistically short timescale for L1544 given that it already displays clear contraction motions, ammonia is heavily depleted in the innermost few thousand au. A few thousand au away from the center, the solutions with chemical desorption included are already at least an order of magnitude above the observations in peak ammonia abundance, while the solutions without chemical desorption fall short of the observations by about an order of magnitude. In the outer core the difference of the model results as compared to the observed profile is greater still. If the gas is let to evolve to $t=10^{6}\,\rm yr$, none of the solutions show an acceptable agreement with the observations even when chemical desorption is excluded from the model.

We can deduce from the results presented in this Section that: 1) we cannot obtain with our physical and chemical models an ammonia abundance profile that does not present heavy ammonia depletion in the central areas of L1544; 2) the observations of Crapsi et al. (2007) can be reproduced to a satisfactory degree, except inside the central few thousand au, if we exclude chemical desorption from the modelling, assume a sufficient amount of visual extinction in the parent cloud, and retain the canonical value for the CR ionization rate.

4 Discussion

4.1 Ammonia line emission simulations

The comparison between the models and the observations presented above does not take into account any optical depth or excitation effects, which may affect for example the determination of the abundance from the observations. Therefore, the abundance profiles calculated here may not exactly correspond to the one derived by Crapsi et al. A more rigorous method of comparing models and observations is the reproduction of the emission lines with radiative transfer methods. To alleviate the ambiguity related to optical depth or excitation effects, we simulated the observed para-$\rm NH_{3}$ (1,1) inversion line and the ortho-$\rm NH_{3}$ ($1_{0}-0_{0}$) ground-state line with the non-local-thermal-equilibrium (LTE) radiative transfer code Cppsimu (Juvela, 1997). For the ortho-line calculations we used collisional coefficients from Bouhafs et al. (2017) which take the hyperfine splitting explicitly into account. Similar hyperfine-split collisional rates are not yet available for para-$\rm NH_{3}$, and so for the para-line calculations we used the data of Maret et al. (2009) along with the assumption that the hyperfine components are split according to LTE.

To test our radiative transfer setup, we attempted to reproduce the observations of Caselli et al. (2017) using two different source profiles. First, we took the density, temperature, and $\rm NH_{3}$ abundance profiles given by Crapsi et al. (2007) and used them as input to Cppsimu. The infall velocity profile was taken from K10, but multiplied by a factor of 1.75 (Bizzocchi et al., 2013). Second, we took the physical structure from K10 and mapped the ammonia abundance profile from Crapsi et al. (2007) onto this physical model, i.e., both models use the same parametrization for the ammonia abundance. The beam FWHM and spectral resolution were set to 40$\arcsec$ and 64 m s*-1*, respectively, corresponding to the Herschel observations of Caselli et al. (2017). The results of this comparison are shown in Fig. 7. We cannot reproduce the observed profile with either one of the two source structures. We point out that MOLLIE (Keto, 1990; Keto & Rybicki, 2010), which was used for the line simulations presented in Caselli et al. (2017), is able to match the observation. We explore this issue further in Appendix A, where we compare the results from the two codes in a couple of test cases. In what follows, we use Cppsimu.

To illustrate the emission lines associated with the abundance profiles presented in this paper, we used two different chemical schemes: 1) our fiducial model, and 2) the fiducial model with chemical desorption turned off and external $A_{\rm V}$ set to 5 mag. The latter model was chosen on the grounds that it provides a decent fit to the observed ammonia abundance profile (outside the core center) as discussed above. From here on we refer to these two models simply as CM1 and CM2. We considered 51 time steps logarithmically evenly spaced between $10^{4}$ and $10^{6}$ yr and searched for the closest match to the observed abundance profile in model CM2 using a $\chi^{2}$ analysis. The best-fit abundance profile as determined by this analysis ($t=1.74\times 10^{5}\,\rm yr$) is shown in Fig. 8.

Fig. 9 shows the results of the radiative transfer calculations for para-$\rm NH_{3}$. The Local Standard of Rest (LSR) velocity of L1544 is 7.2 km s*-1* (Tafalla et al., 1998). The spectral resolution of the simulation was set to $0.1\,\rm km\,s^{-1}$, and we limited the LSR range in the figure to encompass the three central components of the line as in Fig. 2 in Crapsi et al. (2007). The critical density of the (1,1) transition is $\sim 2.0\times 10^{3}\,\rm cm^{-3}$ at 10 K, which means that the line can be collisionally excited in a broad region inside $\sim$25000 au (see Fig. 1). This implies that the (1,1) inversion transition does not probe the central parts of the core well. If we choose model CM1, the beam size does not play a large role because the ammonia abundance is high and hence the column towards the core center is large even when smoothed to 37$\arcsec$. If we instead choose model CM2 where the ammonia abundance profile more or less follows the one derived by Crapsi et al. (2007) – except in the center of the core – the effect of the beam size is accentuated.

Neither one of the chemical models provides a good fit to the observations of Tafalla et al. (2002) and Crapsi et al. (2007). Model CM1 reproduces the main peak intensity of the Effelsberg observations but the satellite lines are too strong, which may also be due to the approximation of LTE for the distribution of the hyperfine components. The optical thickness in the simulations is clearly higher than observed. The VLA observations are not reproduced with this model either, as the intensities of the hyperfine components are overestimated by about 2 K compared to the observations (unfortunately the spectra of Crapsi et al. (2007) are not available to facilitate easy comparison). Model CM2 spectacularly fails to reproduce either observation owing to the missing ammonia in the innermost 5000 au. This shows that even though the (1,1) transition does not probe the core center specifically, emission from the central areas is still crucial when the ammonia abundance in the outer core is low.

Fig. 10 shows the simulated lines for ortho-$\rm NH_{3}$. The critical density of the ($1_{0}-0_{0}$) line at 10 K is $\sim 3.7\times 10^{7}\,\rm cm^{-3}$, and so the line is collisionally excited only in the dense inner part of the core. This feature is especially evident in model CM2 where the ammonia abundance is low in the central area and the resulting emission is very faint. Model CM1 however produces a line that is much too bright and optically thick compared to the observations, and so once again the two chemical models fail to reproduce the observations.

Model CM2 is one example case in a family of solutions that is close to the abundance profile derived by Crapsi et al. (2007), except in the inner core. While we did not calculate $\chi^{2}$ values for all of our models, it is evident from Figs. 2 and 4 and from the analysis presented above that we do not obtain a solution that could reproduce the observed ammonia line profiles, even taking into account possible calibration errors in the observations. The presence of ammonia in the gas phase in the high-density central regions of the core is required.

4.2 Distributions of chemical species related to ammonia formation

As already alluded to in Sect. 3.2, the gas-phase formation efficiency of ammonia depends on the abundance of $\rm N_{2}$: the reaction $\rm N_{2}+He^{+}\longrightarrow N^{+}+N+He$ produces $\rm N^{+}$ which can then be converted to $\rm NH^{+}$ through $\rm N^{+}+H_{2}\longrightarrow NH^{+}+H$. The latter reaction is heavily dependent on the $\rm H_{2}$ ortho/para ratio because it is strongly endothermic when the reaction partner is para-$\rm H_{2}$, but close to thermoneutral with ortho-$\rm H_{2}$ (Dislaire et al., 2012). Therefore, ammonia formation should be the most efficient when there is a large amount of $\rm N_{2}$ in the gas and when the $\rm H_{2}$ ortho/para ratio is high. $\rm N_{2}$ itself cannot be observed in pre-stellar cores, but $\rm N_{2}H^{+}$ can be used as a proxy of its abundance distribution. We plot in Fig. 11 the distributions of $\rm N_{2}H^{+}$ and $\rm N_{2}$, and the $\rm H_{2}$ ortho/para ratio, in our fiducial model at three different time steps. $\rm N_{2}$ freezes out onto the grain surfaces in the central core in a relatively short timescale, and this behavior is clearly reflected in the abundance of $\rm N_{2}H^{+}$. We note that the shape of the $\rm N_{2}H^{+}$ abundance profile does not strictly follow that of $\rm N_{2}$, because $\rm N_{2}H^{+}$ formation requires $\rm H_{3}^{+}$ which in turn reacts preferentially with CO while it is still available in the gas phase. The initial $\rm H_{2}$ ortho/para ratio is already low ($10^{-3}$; Table 1) which means that the production of $\rm NH^{+}$ and hence of ammonia is hindered. The ratio can go as low as $\sim 10^{-5}$ in the central core. It is possible that we are underestimating the initial $\rm H_{2}\,\rm o/p$ ratio, but the single model S3 shows that the initial value has little bearing on the ammonia abundance at high density and low temperature.

The $\rm H_{2}$ o/p ratio that our model predicts is in line with previous observations and models (for targets other than L1544; e.g., Brünken et al., 2014; Furuya et al., 2016; Bovino et al., 2017), and hence there is little reason to believe that an unexpectedly high o/p ratio would be present in the center of L1544, boosting the formation of ammonia beyond the levels obtained with our current model. We also note that $\rm N_{2}H^{+}$ and, to a lesser extent, $\rm N_{2}D^{+}$ depletion has been recently observed in L1544 (Redaelli et al., subm.). It is all the more puzzling why some species related to $\rm N_{2}$ show signs of freeze-out ($\rm N_{2}H^{+}$), while others do not (ammonia).

4.3 Outlook on future modelling efforts

Our chemical model does not take into account some effects that may influence the ammonia abundance. One of these is the treatment of multilayer ice chemistry coupled with dynamic chemical desorption efficiencies depending on the composition of the surface, such as in the recent study of Vasyunin et al. (2017). Their description of chemical desorption is based on the work of Minissale et al. (2016), who unfortunately did not estimate desorption efficiencies for the reactions involved in the ammonia formation through hydrogenation (starting from N + H). While we cannot make a reliable test of the effect of chemical desorption on the relevant hydrogenation reactions owing to the lack of quantitative experimental and theoretical data, we did test the effect of the $\rm N+N\longrightarrow N_{2}$ reaction in our fiducial model by setting the desorption efficiency to the theoretical value of 89% (for bare grains) given by Minissale et al. (2016). We find that the N + N reaction is too slow for the enhanced $\rm N_{2}$ chemical desorption to be of consequence, and indeed the influence of this change on the gas-phase ammonia abundance is negligible. A more complete test can be carried out once experimentally or theoretically derived chemical desorption efficiencies for the appropriate hydrogenation reactions become available.

Another issue is the question of dynamic binding energies on the grain surface. At early times the grains will be covered mostly with water ice, but later as CO starts to freeze out, the grains will be covered with CO and other species. This will change the binding energies of the various species on the surface, and since ammonia is a late-type molecule, it is conceivable that the inclusion of this effect would lead to a decreased depletion factor for ammonia since the binding energy of ammonia on an apolar molecule such as CO will be lower than on water. It is however difficult to formulate this issue in a chemical model in a physically meaningful way, and experimental data is also lacking. Further studies, both theoretical and experimental, are certainly called for.

We demonstrated that considering a self-consistent hydrodynamical treatment of the core collapse coupled with chemistry and radiative transfer (an update of the model presented in Sipilä & Caselli 2018) does not solve the problem of strong ammonia depletion in the central core. Nevertheless it is evident that the dynamics needs to be taken into account in order to track the time-evolution of the abundances of the various species properly.

L1544 is a well-studied object that presents clear contraction motions and no signs of a central source, which implies that it is a pre-stellar core in its final stages of evolution towards becoming a protostellar system. We therefore expect the central area of the object to be cold, and since it is also well-shielded from external radiation which could contribute to the chemistry via photodesorption for example, we come to the conclusion that the inability of the present chemical model to reproduce the observed ammonia abundance is due to the still limited understanding of the chemical and physical processes at play. Time-dependent models present the most promising avenue of further study into this issue. We note that especially large uncertainties pertain to the binding energies of the various species on different types of ice, and that more experimental and theoretical work on this problem is urgently needed so that chemical models of interstellar chemistry can provide reliable results.

5 Conclusions

We performed a parameter-space exploration of the abundance of (ortho+para) ammonia predicted by gas-grain chemical models with different assumptions of parameter values and included/excluded processes. Our goal was to determine if the abundance profile of ammonia observed toward the pre-stellar core L1544, which shows an inward-increasing trend, can be reproduced by considering variations of standard chemical model parameters. For the sake of simplicity and to keep the computational time at a reasonable level, we used in most of our models a static physical structure for L1544 from K10. We varied six different modelling parameters (such as the ammonia binding energy and the sticking coefficient of atomic N) that a priori should influence the ammonia abundance in a significant way, resulting in a total of 108 models. We also considered eleven other models in which we tested the influence of other important parameters such as the initial $\rm H_{2}$ ortho/para ratio, and the effect of dynamics.

Observations of the ammonia abundance toward L1544 by Crapsi et al. (2007) indicate that the abundance increases monotonically toward the center of the core. We found that irrespective of the various parameters we cannot obtain such a profile with our models. The various parameter combinations yield results with varying degrees of ammonia depletion in the central area of the core depending also on the time, but the depletion always occurs – even in a purely gas-phase model, where the effect is due to processes other than adsorption. Interestingly, the models where chemical desorption (which is here modeled assuming a uniform 1% efficiency) is taken into account are seemingly ruled out by our results, as these models lead to solutions where the ammonia abundance in the outer core is orders of magnitude above the observed one. This would constrain the chemical desorption efficiency to below 1%, at least for ammonia on water ice. We also confirmed with radiative transfer simulations that the emission lines arising from our modeled abundance profiles cannot match the observed ones.

Our results point toward a dynamic nature of the chemistry and of the underlying physical processes. On the one hand, static physical models naturally cannot account for the abundance variations caused by time-evolution of the core density profile, and this will be reflected on observable emission lines (Sipilä & Caselli, 2018). On the other hand, considering the effect of the dynamically-varying chemical composition of the ice surface can have a great impact on the efficiency of chemical desorption, as Vasyunin et al. (2017) have recently shown. In addition to chemical desorption, dynamically-varying chemical abundances affect a multitude of other processes, such as line cooling and self-shielding of molecules such as $\rm H_{2}$, CO, and $\rm N_{2}$. It is increasingly evident that simplified pseudo-time-dependent models of interstellar chemistry provide only a limited explanation of the chemical complexity that is observed in the ISM, and that future modelling efforts should concentrate on time-dependent effects. Also, laboratory measurements of the binding energy of ammonia onto ices with variable compositions are required.

Acknowledgements

We thank the anonymous referee for valuable comments. M.J. acknowledges the support of the Academy of Finland Grant No. 285769. P.C. acknowledges the support of the ERC Advanced grant PALs 320620.

Appendix A Radiative transfer code benchmark

Here we compare simulated emission lines calculated with Cppsimu or MOLLIE. We focus first on the ortho-$\rm NH_{3}$ ground-state rotational transition. Figure 12 shows the results for the Crapsi et al. (2007) physical structure (i.e., medium density and gas temperature), and Fig. 13 shows the corresponding calculation for the K10 physical structure. The Crapsi et al. (2007) density structure, which is given as a parametrized formula, was set so that it corresponds to the same outer radius as the K10 model. In both cases the ortho-$\rm NH_{3}$ abundance profile was taken from Crapsi et al. (2007). The infall velocity profile was adopted from K10; we also studied the effect of scaling up the velocity profile by a factor of 1.75 (see Bizzocchi et al., 2013).

MOLLIE fits the observed line profile fairly well when the Crapsi et al. (2007) physical structure is used, whereas Cppsimu overestimates the brightness of the strongest component by almost a factor of two. Curiously, when we switch to the K10 physical structure, MOLLIE overestimates the emission while Cppsimu underestimates it. An extra feature at $v_{\rm LSR}=7.4\,\rm km\,s^{-1}$ that is not seen in the observations is predicted by both codes when the velocity field of K10 is not scaled up. Despite the evident discrepancy, the results from the two codes are within factor of two.

We also compared the line emission profiles of the (1-0) rotational transitions of $\rm HCO^{+}$ and $\rm DCO^{+}$, excluding or including hyperfine structure for $\rm DCO^{+}$ (Fig. 14). The K10 physical structure (with 1.75 velocity scaling) was used for these tests, and we took the appropriate spectral and collisional data from LAMDA (Schöier et al., 2005). The $\rm HCO^{+}$ line profiles match well between the two codes, but for $\rm DCO^{+}$, MOLLIE produces brighter lines, as already seen with ortho-$\rm NH_{3}$ when using the K10 physical structure. The difference in line brightness does not appear to depend on the inclusion or exclusion of hyperfine structure, indicating that the differences in the ammonia spectra provided by Cppsimu and MOLLIE are not due to an incompatible treatment of hyperfines.

We have not found an obvious cause for the discrepancies in the radiative transfer simulations, although it is difficult to model well highly optically thick lines with non-LTE radiative transfer codes (see also Quénard et al., 2016). It is nevertheless very unlikely that our general conclusions on the modeled vs. observed lines, presented in the main body of the paper, would be affected by a different choice of the radiative transfer code used to produce the simulated lines (i.e., adopting a code other than Cppsimu or MOLLIE). Based on our tests we cannot conclude that MOLLIE yields the “correct” result, given that the positions and relative strengths of some spectral features are better reproduced by Cppsimu. Also, the MOLLIE results are counterintuitive in that one would – naively – expect weaker emission in the high-critical-density ortho-ammonia line when using the K10 physical model which is highly centrally concentrated (owing to self-absorption), while MOLLIE predicts the opposite.