An Efficient Method for Determining the Chemical Evolution of Gravitationally Collapsing Prestellar Cores

Dense (≥10⁴ cm⁻³), cold (∼10 K) regions in molecular clouds are the reservoir of gas for low-mass star formation. While the chemistry in these regions has been studied for many years and is broadly understood, the chemistry of star-forming cores remains controversial. The distribution of core masses within the molecular clouds is similar to the observed stellar initial mass function (Motte et al. 1998), and young low-mass stars are frequently found to be associated with dense cores (Cohen & Kuhi 1979). The process by which a dense core evolves into a star can be partially deduced by observations of molecular line emission, and additionally by dust observations later on in the evolutionary process as the gas heats up. Nevertheless, the very early stage of low-mass star formation, when collapse begins, is elusive and there is no direct measurement of the mode(s) of collapse that the core undergoes. Of course, in the absence of any pressure force, a cloud of gas collapses in one free-fall time, t_ff = $\sqrt{\tfrac{3\pi }{32G\rho }}$, but in reality, pressure gradients within a cloud will act to inhibit collapse. A common solution to the hydrostatic equilibrium equation is the well known Bonnor–Ebert (BE) sphere (Bonnor 1956), which describes the density profile of a pressure-confined, self-gravitating isothermal sphere in hydrostatic equilibrium. BE spheres are unstable to gravitational collapse if the ratio between central and outer densities exceeds a critical value, which leads to a critical mass for instability. Kandori et al. (2005) were able to match the density profiles of prestellar cores with BE spheres close to this critical mass, while Kirk et al. (2005) found that for brighter cores, the observed density profiles did not match critical BE spheres. This suggests that while cores may pass through a phase similar to an unstable, thermally supported sphere, it is not sufficient to describe their entire evolution.

In fact, a significant number of dense cores are observed to be rotating (Goodman et al. 1993), and as angular momentum must be conserved during collapse, the infalling gas would be expected to form a rotating disk rather than collapsing directly into the center. Simulations of the collapse of rotating gas clouds (Norman et al. 1980; Matsumoto et al. 1997) find that although there is a complex structure of periodic shock wave formation at the center, the overall picture is of a disk forming in the plane of rotation, with the central density and flatness increasing over time. While rotation does not seem to change the qualitative nature of the collapse, it may be important in determining the subsequent evolution of the central object. Rotation also causes the "angular momentum problem," which is that, theoretically, rotating prestellar cores have much more angular momentum than a star could contain without breaking up (Prentice & ter Haar 1971). Proposed solutions include the fragmentation of the core into multiple protostellar objects (Matsumoto & Hanawa 2003) and the removal of angular momentum by magnetic fields, known as magnetic braking (Basu & Mouschovias 1994). Simulations of the fragmentation and collapse of a magnetized filament (Nakamura et al. 1995; Tomisaka 1996) show that, as with rotation, once a cloud has begun to collapse dynamically, magnetic forces are not sufficient to halt it. However, their conclusions depend on the coupling between the magnetic field and the gas, and this is determined by the fractional ionization, which itself is controlled by the chemical evolution during the collapse and not included in these works.

The density structure of prestellar cores can be inferred from observations of either the dust continuum, or line emission from molecules. Both methods require a conversion factor to relate either the dust emissivity or the column density of the observed molecular species to the total gas density; these conversion factors introduce systematic uncertainties. Additionally, at typical distances, observations of the central regions are limited by resolution, while in the low-density outer regions, the lower signal-to-noise is a further source of uncertainty. The differences in density structure between many proposed models of collapse are often not large enough for observations to conclusively favor one over the others. However, the timescale for molecular abundances in a core to reach equilibrium levels are comparable to the timescales involved in gravitational collapse, which means that the chemical composition of a cloud should be sensitive to the hydrodynamical situation. Simulations of the chemical evolution of a core undergoing collapse (Rawlings et al. 1992; Bergin & Langer 1997; Aikawa et al. 2001) have found that this is the case, with different hydrodynamical models giving significant differences in the abundances of some molecules. This raises the possibility that molecular abundances could be used as an observational test of theories of star formation. This is our motivation.

However, there are difficulties in coupling detailed chemical models to a hydrodynamical code, especially when magnetic fields are included, and simpler chemical models cannot necessarily be relied on to give accurate molecular abundances. For example, Aikawa et al. (2001) and Rawlings et al. (1992) used analytical solutions for isothermal collapse, found by Larson (1969) and Shu (1977), respectively, but these are not necessarily applicable to real situations—the Larson (1969) solution only agrees with simulations at small radii (Hunter 1977), while the Shu (1977) solution begins from a static singular isothermal sphere, which is unlikely to be a realistic initial state. Lee et al. (2004) obtained the initial conditions for the Shu (1977) solution using a series of progressively denser BE spheres evolving quasistatically, before continuing to follow the chemical evolution throughout the subsequent dynamical collapse. Aikawa et al. (2005) improved on this earlier work and performed full hydrodynamical calculations of the collapse of BE spheres, while Li et al. (2002) used a simplified, one-dimensional model of magnetic forces, to make the simulation computationally feasible. Hincelin et al. (2013, 2016) followed the chemical evolution of tracer particles in a fully three-dimensional magnetohydrodynamical (MHD) simulation, assuming that the neutral gas is perfectly coupled to the magnetic field. Tassis et al. (2012) used the thin disk approximation to investigate the effects of nonideal MHD, where the neutral and ionized parts of the gas behave differently. Keto et al. (2015) used a simplified chemical network to investigate the effects of different collapse modes on the line profiles of CO and H₂O, finding that quasi-equilibrium contraction of a BE sphere was best able to reproduce observations.

In this paper, we propose a different approach: we parameterize the results of hydrodynamical simulations of collapsing prestellar cores to describe how the density behaves as a function of radius and time for different models, and incorporate these parameterizations into a gas-grain time-dependent chemical model. Although less accurate than a simultaneous solution of both the hydrodynamics and chemistry, this approach removes the need for simplifications in either area, while also being much less computationally expensive, and so enabling the exploration of larger regions of parameter space than has so far been feasible.

This paper is laid out as follows: In Section 2, we describe our parameterizations of hydrodynamical simulations, and we discuss the chemical model into which we incorporate them. In Section 3, we present the results of our grid of models, showing how the abundances of key molecules are affected. In Section 4, we investigate whether the inclusion of ambipolar diffusion (AD) affects star formation timescales; we briefly discuss our findings and conclude in Section 5. In Appendices A and B, we list the functions used in our parameterization of the numerical simulations.

2.1. Parameterization of Numerical Simulations

Empirical models were developed to reproduce the results from four numerical simulations of collapsing prestellar cores: Aikawa et al. (2005) used a BE sphere as the initial configuration, with the density increased by a factor of either 1.1 (model BES1) or 4 (BES4) to take the core out of equilibrium and instigate collapse; Nakamura et al. (1995) studied the collapse of a magnetically supported filament when a density perturbation is applied (MS); and Fiedler & Mouschovias (1993) followed the evolution of a core as magnetic support (MS) is removed through AD.

Each of these studies produced the density profile (number density for BES1, BES4 and AD, mass density for MS) of the core as a function of time during the collapse. We extracted data from published plots of the density profiles at different times, as shown in Figure 1 with data from Aikawa et al. (2005) as an example.

Zoom In Zoom Out Reset image size

For each density profile, a function, depending on position in the core, reproducing its shape was found. In the two models including magnetic effects (MS and AD), density profiles were given for both the radial (z = 0) and z axes. Only the density profile in the radial direction was considered, as the core rapidly collapses into a thin disk so that structure along the z-axis is less significant. For the BES1, BES4, and AD models, a function of the form

$$\begin{eqnarray}&&n(r)=\displaystyle \frac{{n}_{0}}{1+{\left(\tfrac{r}{{r}_{0}}\right)}^{a}}\end{eqnarray} \tag{ 1 }$$

with n₀, r₀, and a free parameters determining the central density, the width of the central density peak, and the slope of the profile, provided a good fit. Tafalla et al. (2002) used Equation (1) to fit the observed density profiles of prestellar cores. For the MS model, the radial profiles could not be approximated by Equation (1). The equilibrium density of the filament is given by an equation of the form

$$\begin{eqnarray}&&\rho (r)={\rho }_{0}{\left(1+{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{2}\right)}^{-a},\end{eqnarray} \tag{ 2 }$$

which was adapted to reproduce the data by changing the outer exponent a, so that the slope at large r is the same as in the simulated data.

The values of the free parameters were chosen to approximate the density profiles at each time point given. Figure 1 shows the result for the Aikawa et al. (2005) data. For each parameter, the time evolution was also approximated by a similar method. For example, the central density parameter's time evolution was found to be well reproduced by

$$\begin{eqnarray}&&\mathrm{log}{n}_{0}(t)=A{({t}_{0}-t)}^{a}+B\end{eqnarray} \tag{ 3 }$$

for all simulations, where t₀ is the simulation's duration. a was chosen such that log n₀ is approximately linear with ${({t}_{0}-t)}^{a}$, and the coefficients A and B can then by found by linear regression. Figure 2 shows the variation of the central density parameter, n₀, with time for the data from Aikawa et al. (2005), along with the resulting approximation to the time dependence of this parameter. Once all of the parameters have been approximated in this way, the density can be calculated as a function of time and space, shown in Figure 3 compared with the original data. Figures 4–6 show the corresponding density profile approximations for the BES4, MS, and AD collapses, respectively. The equations used to approximate the time dependence of the density profiles are given in Appendices A and B. The maximum discrepancies between the simulations and approximated densities are 34%, 240%, 66%, and 53% for the BES1, BES4, MS, and AD cases, respectively, while the average discrepancies are 10%, 26%, 16%, and 16%. The large (>100%) errors in the BES4 approximation occur only at late times and large radii—otherwise the agreement with the data is at a similar level to the other collapses.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our approximations give the time evolution of the density at a given radius. However, during collapse, the individual parcels of gas do not remain at a constant radius, but move inwards, leading to a different density evolution than for fixed r. For the BES1 and BES4 models, we determine the new radius of a parcel at each time step by calculating the mass interior to its initial radius at t = 0, M(<r₀), and finding the radius at the given time which encloses the same mass. The MS and AD approximations are not spherically symmetric, so this approach would require knowledge of the density at each point. Instead, we use the same methods as for the density to approximate the radial velocity profiles and use these to calculate the new parcel radius at each timestep. The gas density versus time, for parcels of differing initial radii, in the BES1 collapse is shown in Figure 7.

Zoom In Zoom Out Reset image size

2.2. Chemical Modeling

The four density approximations were used as input for the chemical code UCL_CHEM. The code is described in Viti et al. (2004) and references therein. This code has since been made public (https://uclchem.github.io/) and is fully explained in Holdship et al. (2017). Here, we briefly summarize its characteristics: UCL_CHEM is a time-dependent gas-grain chemical model that calculates the abundances of atoms and molecules in the gas and dust in the interstellar medium as a function of time under chemical and physical conditions set by the user. The original version of the code uses free-fall collapse to determine the density from the diffuse state to the final density of the gas where the star is born. Initial atomic elemental abundances are provided to UCL_CHEM, which then self-consistently calculates gas-phase chemistry, as well as sticking on to dust particles with subsequent surface processing. For the reaction network, we used the UMIST 2012 network (McElroy et al. 2013), and freeze-out and grain surface reactions as described in Holdship et al. (2017).

We chemically model our four parameterized numerical simulations: the collapse of an unstable (BES1) or highly unstable (BES4) BE sphere, collapse against MS, and collapse resulting from AD. A grid of models was run for each case to investigate the effects of changing other input parameters: the cosmic-ray ionization rate ζ (in units of ζ₀ = 1.3 × 10⁻¹⁷ s⁻¹), metallicity Z, and the desorption efficiency parameters , ϕ and Y_UV, corresponding to the number of molecules desorbed per H₂ molecule formed, per cosmic-ray impact and per UV photon (produced by cosmic rays) absorbed respectively (Roberts et al. 2007). For our fiducial desorption efficiencies, H₂ formation is the dominant desorption mechanism, as found by Roberts et al. (2007). The values adopted for each model are given in Table 1. Each model was run once for each of the density approximations. We assume an external radiation field of 1 Habing and an external extinction at the core boundary of 3 mag, the value used by Aikawa et al. (2005). The extinction from the core itself is calculated by integrating the density profile to the boundary, (0.2 $\mathrm{pc}$ for the BES1 and BES4 approximations, 0.5 pc for MS and 0.75 $\mathrm{pc}$ for the AD case) before the onset of collapse, giving maximum extinctions (at the centermost parcel) of 6.4 (BES1), 15.2 (BES4), 7.8 (MS), and 3.3 (AD) mag. We used 13 gas parcels (14 for AD) at initial radii spaced to cover the entire range of the cores, but with an emphasis on the more rapidly evolving central regions.

Table 1.  Model Input Parameters

Model	ζ/ζ0	Z/Z⊙		ϕ	YUV
A	1	1	0.01	105	0.1
B1	5	1	0.01	105	0.1
B2	10	1	0.01	105	0.1
C1	1	0.3	0.01	105	0.1
C2	1	1.5	0.01	105	0.1
D1	1	1	0.1	105	0.1
D2	1	1	1.0	105	0.1
E1	1	1	0.01	104	0.1
E2	1	1	0.01	106	0.1
F1	1	1	0.01	105	0.001
F2	1	1	0.01	105	1.0

Download table as:  ASCIITypeset image

The initial central number densities of the models are n = 2.2 × 10⁴ cm⁻³ (BES1 and MS), 8 × 10⁴ cm⁻³ (BES4) and 300 cm⁻³ (AD). The MS equations are in terms of dimensionless variables, which can be converted into physical values by choosing the initial central density, ρ_c, and the isothermal sound speed, ${c}_{s}^{2}=\tfrac{{kT}}{\mu {m}_{{\rm{H}}}}$. The gas was assumed to be at a temperature of 10 K and composed entirely of molecular hydrogen for the purposes of calculating the mean molecular mass, giving ${c}_{s}=203\,{\rm{m}}\,{{\rm{s}}}^{-1}$, and we set ${\rho }_{c}=3.67\times {10}^{-20}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ to give an initial central number density equal to the BES1 model. The initial number density at each point is then given by the relevant equation for the density profile, and we allow the chemistry to evolve for 1 Myr before the onset of collapse at this density. The models were run until the central density reached 10⁸ cm⁻³. The elemental abundances relative to H for our standard model, the solar values given by Asplund et al. (2009), are listed in Table 2—the abundances for elements other than H and He in models with varying metallicity are multiplied by the value of Z.

Table 2.  Elemental Abundances

Element	Abundance	Element	Abundance
H	1.0	N	6.8 × 10−5
He	0.085	S	1.3 × 10−5
C	2.7 × 10−4	Si	3.2 × 10−5
O	4.9 × 10−4	Cl	3.2 × 10−7

Download table as:  ASCIITypeset image

3.1. Molecular Abundances across the Different Modes of Collapse

Figure 8 shows the density profiles of the four collapse modes at the point when the central number density, n₀, reaches 2 × 10⁵ ${\mathrm{cm}}^{-3}$. The BES1 and MS profiles decrease more rapidly with distance than for the BES4 and AD approximations, which have similar densities up to the end of the BES4 core at 0.2 $\mathrm{pc}$. However, the time taken to reach this point is much shorter for the BES4 case (∼10⁵ years) than for AD (∼10⁷ years), so the chemical evolution differs significantly despite the gas density being the same. The BES1 and MS modes both reach ${n}_{0}=2\times {10}^{5}\,{\mathrm{cm}}^{-3}$ after ∼10⁶ years, and as such, the chemical evolution is similar. The densities at large (≳0.2 $\mathrm{pc}$) radii for BES4/AD are much higher than for BES1/MS, due to the more rapid collapse of the outer parts of the core in the BES4 case, and the MS of the outer regions for AD.

Zoom In Zoom Out Reset image size

Figure 9 shows the abundances of four molecules versus radii for model A for the four density approximations, at a central density ${n}_{0}=2\times {10}^{5}\,{\mathrm{cm}}^{-3}$. In all cases, the CO abundance increases from a central minimum, where freeze-out has depleted most of the molecule onto grain surfaces, to a maximum value of ∼10⁻⁴. The BES1, BES4, and MS approximations all behave similarly in reaching the maximum, although for BES4 and MS the radius at which this occurs is larger due to the higher gas densities increasing the effect of freeze-out. The CO abundance in the AD collapse increases much more slowly, only reaching 10⁻⁴ at the edge of the core, whereas in the MS case the abundance begins to fall again toward the edge. This is due to the higher densities in the outer regions for the AD collapse, as MS can still prevent material here from collapsing, unlike in the other three cases. The effect of the longer collapse duration is also apparent, as the AD abundances are far lower than the BES4 ones at comparable radii, despite the densities being very similar. The NH₃ abundance also increases from the center to a maximum, before falling with radius in all models. The decline is much more gradual for the AD case than the other collapse modes, leading to an order of magnitude difference with the MS model by r = 0.3 $\mathrm{pc}$. HCO⁺ shows similar behavior, while for HCN the AD collapse mode produces a nearly constant abundance after reaching a value of ∼10⁻⁸, in contrast to the others.

Zoom In Zoom Out Reset image size

3.2. Cosmic-Ray Ionization Rate

Raising the cosmic-ray ionization rate increases the abundances of molecules in the center of the core, regardless of the collapse mode, as the rate of desorption from grain surfaces is increased. Figure 10 shows the CO abundances for models A, B1, and B2, for the BES1 and AD density approximations, at a central density ${n}_{0}=2\times {10}^{5}\,{\mathrm{cm}}^{-3}$. For the BES1 collapse, the abundance at larger radii is mostly unaffected, whereas for AD, the A and B2 models show noticeably different behavior, with the CO abundance an order of magnitude lower at the edge of the core for the B2 model due to the higher cosmic-ray dissociation rate. Figure 11 shows the HCO⁺ abundance for the same models. The central abundances are again enhanced for models with higher ionization rates, but whereas for the AD collapse mode the abundance decreases toward the edge as with CO, for BES1 the abundance is higher throughout the cloud.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3. Metallicity

Changing the metallicity of the core usually results in a corresponding change in the molecular abundances, due to the different availability of atoms to form the molecules. However, some molecules are much less affected than others. Figure 12 shows the abundances of CO and HCN for the BES1 collapse at a central density ${n}_{0}=2\times {10}^{5}\,{\mathrm{cm}}^{-3}$, for models A, C1, and C2. Whereas the CO abundance scales nearly linearly with the metallicity, the HCN abundance is virtually unchanged between models. The main formation and destruction reactions for HCN both involve H⁺, for which the abundance increases with decreasing metallicity (and vice versa), at least partially counteracting the effect of the changing elemental abundances on the HCN abundance.

Zoom In Zoom Out Reset image size

3.4. Desorption Efficiencies

As with the cosmic-ray ionization rate, increasing the desorption efficiencies increases the molecular abundances, particularly in the denser central regions where more freeze-out has taken place. Figure 13 shows the CO abundances for models A, D1, and D2, for the MS and AD density approximations, where the H₂ formation desorption efficiency has been modified. While the abundances in the central regions are affected similarly for both collapse modes, at larger radii the effect is negligible for the MS collapse, whereas the CO abundance reaches 10⁻⁴ much more rapidly for AD—for model D2, the abundance profile looks much more similar to the other collapse modes than for model A. The BES1 and BES4 collapses behave similarly to MS for varying desorption efficiency, with very little change in the CO abundance beyond 0.1 $\mathrm{pc}$. The cosmic-ray heating desorption efficiency ϕ has very little effect on the abundance of any molecule, despite a factor of 100 difference between models E1 and E2. The cosmic-ray induced photodesorption efficiency Y_UV, however, does affect molecular abundances, in particular having a significant effect on the abundance of NH₃, which is not greatly affected by variation of the other parameters investigated.

Zoom In Zoom Out Reset image size

Only one of our modes of collapse, AD, includes the effect of AD. However, all prestellar cores are expected to be magnetized, and therefore AD could be important for all collapse models. An estimate of the timescale on which AD occurs is given by

$$\begin{eqnarray}&&{t}_{\mathrm{amb}}=4\times {10}^{5}({x}_{i}/{10}^{-8})\,\mathrm{years}\end{eqnarray} \tag{ 4 }$$

(Mouschovias 1979; Hartquist & Williams 1989). If t_amb is smaller than the free-fall timescale, magnetic pressure is unlikely to impede gravitational collapse, while if it is larger the impeding effects may be significant.

Banerji et al. (2009) showed that the AD timescale becomes very large as the fractional ionization increases, and the magnetic field is strongly coupled to the collapsing core, which, in some cases, may halt the collapse and hence the formation of the star. We calculated the AD timescale at the center of the core using Equation (4) at the beginning of the collapse, and at central densities of 10⁶ and 10⁸ ${\mathrm{cm}}^{-3}$, for the BES1, BES4, and MS approximations, for differing metallicities and cosmic-ray ionization rates, using ionization fractions calculated in our chemical simulations. We compared these timescales with the time it takes the gas to reach the final density (10⁸ cm⁻³). Table 3 shows the collapse duration and t_amb for this grid of models.

Table 3.  Collapse Duration and Ambipolar Diffusion Timescales at Increasing Density for BES1, BES4, and MS Models with Varying ζ and Z

				tamb/106 years
Model	Z/Z⊙	ζ/ζ0	tcollapse/106 years	Initial	${n}_{{\rm{H}}}={10}^{6}$ cm−3	${n}_{{\rm{H}}}={10}^{8}$ cm−3
BES1 A	1.0	1.0	1.173	0.80	0.09	0.08
BES1 B1	1.0	5.0	1.173	2.51	0.14	0.09
BES1 B2	1.0	10.0	1.173	4.70	0.19	0.10
BES1 C1	0.3	1.0	1.173	0.62	0.09	0.08
BES1 C2	1.5	1.0	1.173	0.95	0.09	0.07
BES4 A	1.0	1.0	0.184	0.22	0.08	0.07
BES4 B1	1.0	5.0	0.184	0.72	0.12	0.08
BES4 B2	1.0	10.0	0.184	1.35	0.18	0.09
BES4 C1	0.3	1.0	0.184	0.22	0.08	0.07
BES4 C2	1.5	1.0	0.184	0.22	0.08	0.07
MS A	1.0	1.0	1.393	0.73	0.09	0.05
MS B1	1.0	5.0	1.393	2.39	0.14	0.08
MS B2	1.0	10.0	1.393	4.49	0.19	0.09
MS C1	0.3	1.0	1.393	0.58	0.09	0.06
MS C2	1.5	1.0	1.393	0.85	0.09	0.06

Download table as:  ASCIITypeset image

For all models considered, the value of t_amb at the final density, ${n}_{{\rm{H}}}={10}^{8}$ cm⁻³, is significantly lower than the collapse duration, while the initial values are comparable to or larger than the collapse time, particularly for the BES4 collapse, and for the models with increased cosmic-ray ionization rates. At 10⁴ cm⁻³, the B1 and B2 models have t_amb larger than the collapse time for the BES4 case. Increasing ζ leads to larger values of t_amb, as the ionization, and therefore the coupling to the neutral gas, is increased. Metallicity has very little effect at higher densities, but the initial t_amb varies with the metallicity, as more readily ionized atoms such as carbon are present. The BES1 and MS models have lower initial values of t_amb/t_collapse than the BES4 models, suggesting that the faster collapse should be impeded more strongly by the coupling of gas to magnetic fields. These results emphasize that AD is important for diffuse material, where magnetic fields are likely to impede collapse, but once denser clumps have formed, MS will be removed too rapidly to affect the subsequent evolution.

The results of our chemical modeling show that it is difficult to disentangle the effect on molecular abundances of different collapse modes from that of varying the input parameters, which are not necessarily known a priori. Drawing information about the dynamics of star formation from molecular abundances therefore requires a full investigation of parameter space, something which would be extremely time consuming using combined hydrodynamical-chemical modeling. Our grid of models, although not large enough to draw robust conclusions about individual objects, does allow us to compare results with the general properties of prestellar cores, and determine whether particular collapse models or regions of parameter space are in conflict with observation.

Assuming the values from model A (see Table 1), all density approximations predict CO abundances away from the core center in agreement with observed values of 10⁻⁵–10⁻⁴ in starless cores (Caselli et al. 1999; Frau et al. 2012, assuming ¹⁸O/¹⁶O ≈ 10⁻³ and ¹⁷O/¹⁶O ≈ 10⁻⁴). However, the AD collapse only reaches these values at r ≳ 0.3 $\mathrm{pc}$, much larger than typical core sizes (<0.1 $\mathrm{pc}$; Frau et al. 2012). Only models with the highest desorption efficiencies investigated for H₂ formation and cosmic-ray induced photodesorption (D2 and F2) predict CO abundances of ∼10⁻⁵ at a radius of 0.1 $\mathrm{pc}$ for an AD collapse.

The BES1, BES4, and MS approximations result in similar abundance profiles for most molecules. The BES4 and MS abundances are generally more depleted in the center than the BES1 ones, as the gas densities are higher due to either MS, or higher initial densities and a more rapid collapse, causing more efficient freeze-out. However, these differences are not large enough to provide a robust observational test of the mode of collapse. The AD approximation produces significantly different profile shapes to the other three, due to both the longer collapse duration and the higher densities at large radii. The slower increase with radius of the abundance of molecules such as CO, HCN, and HCO⁺, and the subsequent slow or negligible decline beyond the peak value, are qualitatively different to the situation with the initially unstable collapse modes, suggesting spatially resolved molecular observations could be used to discriminate between them.

All density approximations predict peak NH₃ abundances of ∼10⁻⁸, consistent with observed values in prestellar cores (Tafalla et al. 2002; Johnstone et al. 2010). At ${n}_{0}=2\times {10}^{5}\,{\mathrm{cm}}^{-3}$, N₂H⁺ abundances, shown in Figure 14, are intermediate between the values found by Frau et al. (2012) (∼10⁻¹¹) and Tafalla et al. (2002) and Johnstone et al. (2010) (∼10⁻¹⁰). However, as with CO, for the AD approximation the abundances within the reported radii of the cores are much lower than the observed values. Ambipolar diffusion simulations including chemistry by Tassis et al. (2012) show similar behavior, with the inclusion of magnetic effects leading to a more extended depleted region than in unmagnetized models, suggesting that this is a genuine feature of collapse under AD, rather than being due to our approximation of the density evolution.

Zoom In Zoom Out Reset image size

Lee et al. (2004) calculated the chemical evolution of a quasistatically contracting BE sphere over 10⁶ years, finding similar abundance profiles to our BES1 approximation for CO, NH₃, and HCO⁺, although for N₂H⁺, we find much more depletion in the core centers, and higher HCN abundances overall. Given that our MS and BES4 approximations also give similar results to BES1, this suggests that the collapse timescale, rather than the specific details of the density evolution, are more important for the chemical evolution, at least for some observationally important molecules such as CO.

Our BES1 and BES4 approximations are based on hydrodynamical simulations presented in Aikawa et al. (2005), who modeled the chemical evolution of these models self-consistently, providing a test of the approximations' accuracy. Comparing to their results, we find good agreement (of the same order of magnitude) between the predicted peak molecular abundances of both approaches. However, the variation with radius differs—our approximations generally predict lower abundances at the core centers, especially for the BES4 collapse, where our results predict lower CO abundances to the BES1 case at the same r, whereas Aikawa et al. (2005) find the opposite. We attribute this to the differing initial conditions—Aikawa et al. (2005) assume the gas is entirely atomic prior to collapse, whereas we allow the abundances to evolve for 10⁶ years at the initial density before the onset of collapse, leading to significant freeze-out onto grains occurring in the denser central regions.

We conclude that despite the dependence of molecular abundances on the various parameters mentioned above, molecular observations can still be useful for discriminating between different models of collapse. While models that begin from an initially gravitationally unstable state predict similar abundances and radial variations, AD produces qualitatively different abundance profiles for many observationally important molecules, which appear to be in conflict with observations of prestellar cores, although varying the input parameters may be able to reduce this discrepancy. A more exhaustive investigation of parameter space, combined with observations of multiple species from the same source, could be used to draw much stronger conclusions on the nature of core collapse, as well as providing constraints on the values of the input parameters, which are currently assumed ad hoc. This sort of investigation would be extremely time consuming, if not impossible, using a coupled hydrodynamical-chemical system, even without the additional complications of magnetic fields. The results we have presented here demonstrate that these large grids of models are now feasible using our method of parameterizing the dynamics, while still providing a reasonable level of accuracy compared to full simulations.

F.P. is supported by the Perren fund and the IMPACT fund. We would like to thank Prof. Fred Adams for his helpful contribution, and the referee, whose comments helped to produce a greatly improved version of the original paper.

Software: UCL_CHEM (Viti et al. 2004; Holdship et al. 2017).

The BES1, BES4, and AD collapse density profiles were approximated with the function

$$\begin{eqnarray}&&n(r)=\displaystyle \frac{{n}_{0}(t)}{1+{\left(\tfrac{r}{{r}_{0}(t)}\right)}^{a(t)}}\end{eqnarray} \tag{ 5 }$$

whereas the MS collapse required a different functional form,

$$\begin{eqnarray}&&\rho (r)={\rho }_{0}(t){\left(1+{\left(\displaystyle \frac{r}{{r}_{0}(t)}\right)}^{2}\right)}^{-a(t)}\end{eqnarray} \tag{ 6 }$$

where n₀(t), ρ₀(t), r₀(t), and a(t) are functions of time since the onset of collapse given in the following subsections.

A.1. BES1

The time-dependent parameters are given by

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{n}_{0}(t)=61.8{(1.175\times {10}^{6}-t)}^{-0.01}-49.4\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{r}_{0}(t)=-28.5{(1.175\times {10}^{6}-t)}^{-0.01}+28.93\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&a(t)=2.4\end{eqnarray} \tag{ 9 }$$

where n₀ is in ${\mathrm{cm}}^{-3}$, r₀ is in au and t is in years.

A.2. BES4

The time-dependent parameters are given by

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{n}_{0}(t)=68.4{(1.855\times {10}^{5}-t)}^{-0.01}-55.7\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{r}_{0}(t)=-39.0{(1.855\times {10}^{5}-t)}^{-0.01}+38.7\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&a(t)=1.9+0.5\exp (-t/{10}^{5})\end{eqnarray} \tag{ 12 }$$

where n₀ is in ${\mathrm{cm}}^{-3}$, r₀ is in au and t is in years.

A.3. MS

The time-dependent parameters are given by

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{\rho }_{0}(t)=3.54{(5.47-t)}^{-0.15}-2.73\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{r}_{0}(t)=-1.34{(5.47-t)}^{-0.15}+1.47\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&a(t)=2.0-0.5{\left(\displaystyle \frac{t}{5.47}\right)}^{9}\end{eqnarray} \tag{ 15 }$$

with the units determined by the initial central density ρ_c and the sound speed, c_s. ρ₀ is in units of ρ_c, r₀ in units of $\tfrac{{c}_{s}}{\sqrt{2\pi G{\rho }_{c}}}$ and t in units of ${(2\pi G{\rho }_{c})}^{-0.5}$.

A.4. AD

The time-dependent parameters are given by

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{n}_{0}(t)={\mathrm{log}}_{10}(2+1.7\left(\displaystyle \frac{t}{6}-1\right))+3\quad t\lt 6.0\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&5.3{(16.138-t)}^{-0.1}-1.0\quad t\geqslant 6.0\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{r}_{0}(t)=-2.57{(16.138-t)}^{-0.1}+1.85\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&a(t)=2.4-0.2{\left(\displaystyle \frac{t}{16.138}\right)}^{40}\end{eqnarray} \tag{ 19 }$$

where n₀ is in ${\mathrm{cm}}^{-3}$, r₀ is in 0.75 $\mathrm{pc}$ and t is in Myr.

For the MS and AD collapses, the radial velocity profiles were also approximated in order to determine the inwards movement of the gas parcels.

B.1. MS

The radial velocity profile is given by

$$\begin{eqnarray}&&{v}_{r}(r)={v}_{\min }(t)\left[{\left(\displaystyle \frac{r^{\prime} (t)}{{r}_{\min }(t)}\right)}^{2}-1\right]\end{eqnarray} \tag{ 20 }$$

for r < r_min and

$$\begin{eqnarray}&&{v}_{r}(r)={v}_{\min }(t)[\exp (-2a(t)r^{\prime} (t))-2\exp (-a(t)r^{\prime} (t))]\end{eqnarray} \tag{ 21 }$$

for r ≥ r_min, where r'(t) = r − r_min. The time evolution is given by

$$\begin{eqnarray}&&-1.149t+7.2\quad t\lt 4.95\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{r}_{\min }(t)=-9.2\mathrm{log}t+16.25\quad t\lt 5.33\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&-22\mathrm{log}t+37.65\quad t\geqslant 5.33\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&0.0891t\quad t\lt 4.95\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{v}_{\min }(t)=5.5\mathrm{log}t-8.37\quad t\lt 5.33\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&18.9\mathrm{log}t-30.8\quad t\geqslant 5.33\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&0.0101t+0.4\quad t\lt 4.95\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&a(t)=0.695\mathrm{log}t-0.663\quad t\lt 5.33\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&2.69\mathrm{log}t-4\quad t\geqslant 5.33\end{eqnarray} \tag{ 30 }$$

where v_min is in units of c_s, r_min in units of $\tfrac{{c}_{s}}{\sqrt{2\pi G{\rho }_{c}}}$ and t in units of ${(2\pi G{\rho }_{c})}^{-0.5}$.

B.2. AD

The radial velocity profile is given by

$$\begin{eqnarray}&&{v}_{\min }(t)\left[{\left(\displaystyle \frac{r^{\prime} (t)}{{r}_{\min }(t)}\right)}^{2}-1\right]\quad r\lt {r}_{\min }\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}{v}_{r}(r) & = & ({v}_{\min }(t)-{v}_{\mathrm{mid}}(t)){\left(\displaystyle \frac{r^{\prime} (t)}{0.5-{r}_{\min }(t)}\right)}^{0.3}\\ & & -{v}_{\min }(t)\quad r\lt 0.5\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&2{v}_{\mathrm{mid}}(t)(r-1)\quad r\geqslant 0.5\end{eqnarray} \tag{ 33 }$$

The time evolution is given by

$$\begin{eqnarray}&&-0.0039t+0.49\quad t\leqslant 10.2\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{r}_{\min }(t)=-0.0306(t-10.2)+0.45\quad t\leqslant 15.1\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&-0.282(t-15.1)+0.3\quad t\gt 15.1\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{v}_{\min }(t)=3.44{(16.138-t)}^{-0.35}-0.7\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{v}_{\mathrm{mid}}(t)=0.143t\quad t\leqslant 10.2\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&0.217(t-10.2)+1.46\quad t\gt 10.2\end{eqnarray} \tag{ 39 }$$

where r_min is in units of 0.75 $\mathrm{pc}$, v_min and v_mid are in 10⁻² km s⁻¹ and t is in Myr.