An efficient method for determining the chemical evolution of gravitationally collapsing prestellar cores
aa r X i v : . [ a s t r o - ph . GA ] J un Draft version June 6, 2018
Typeset using L A TEX twocolumn style in AASTeX61
AN EFFICIENT METHOD FOR DETERMINING THE CHEMICAL EVOLUTION OF GRAVITATIONALLYCOLLAPSING PRESTELLAR CORES
F. D. Priestley, S. Viti, and D. A. Williams Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK
ABSTRACTWe develop analytic approximations to the density evolution of prestellar cores, based on the results of hydrodynam-ical simulations. We use these approximations as input for a time-dependent gas-grain chemical code to investigatethe effects of differing modes of collapse on the molecular abundances in the core. We confirm that our methodcan provide reasonable agreement with an exact numerical solution of both the hydrodynamics and chemistry whilebeing significantly less computationally expensive, allowing a large grid of models varying multiple input parametersto be run. We present results using this method to illustrate how the chemistry is affected not only by the collapsemodel adopted, but also by the large number of unknown physical and chemical parameters. Models which areinitially gravitationally unstable predict similar abundances despite differing densities and collapse timescales, whileambipolar diffusion produces more extended inner depleted regions which are not seen in observations of prestellarcores. Molecular observations are capable of discriminating between modes of collapse despite the unknown values ofvarious input parameters. We also investigate the evolution of the ambipolar diffusion timescale for a range of collapsemodes, metallicities and cosmic ray ionization rates, finding that it remains comparable to or larger than the collapsetimescale during the initial stages for all models we consider, but becomes smaller at later evolutionary stages. Thisconfirms that ambipolar diffusion is an important process for diffuse gas, but becomes less significant as cores collapseto higher densities.
Keywords: astrochemistry – stars: formation – ISM: molecules
Corresponding author: F. D. [email protected]
Priestley et al. INTRODUCTIONDense ( ≥ cm − ), cold ( ∼
10 K) regions in molec-ular clouds are the reservoir of gas for low mass starformation. While the chemistry in these regions hasbeen studied for many years and is broadly understood,the chemistry of star-forming cores remains controver-sial. The distribution of core masses within the molecu-lar clouds is similar to the observed stellar initial massfunction (Motte et al. 1998), and young low mass starsare frequently found to be associated with dense cores(Cohen & Kuhi 1979). The process by which a densecore evolves into a star can be partially deduced by ob-servations of molecular line emission, and additionallyby dust observations later on in the evolutionary pro-cess as the gas heats up. Nevertheless the very earlystage of low mass star formation, when collapse begins,is elusive and there is no direct measurement of themode(s) of collapse that the core undergoes. Of course,in the absence of any pressure force, a cloud of gas col-lapses in one free-fall time, t ff = q π Gρ , but in reality,pressure gradients within a cloud will act to inhibit col-lapse. A common solution to the hydrostatic equlibriumequation is the well known Bonnor-Ebert (BE) sphere(Bonnor 1956), which describes the density profile of apressure-confined, self-gravitating isothermal sphere inhydrostatic equilibrium. BE spheres are unstable togravitational collapse if the ratio between central andouter densities exceeds a critical value, which leads to acritical mass for instability. Kandori et al. (2005) wereable to match the density profiles of prestellar cores withBE spheres close to this critical mass, while Kirk et al.(2005) found that for brighter cores, the observed den-sity profiles did not match critical BE spheres. Thissuggests that while cores may pass through a phase sim-ilar to an unstable, thermally supported sphere, it is notsufficient to describe their entire evolution.In fact, a significant number of dense cores are ob-served to be rotating (Goodman et al. 1993), and asangular momentum must be conserved during collapse,the infalling gas would be expected to form a rotat-ing disc rather than collapsing directly into the cen-tre. Simulations of the collapse of rotating gas clouds(Norman et al. 1980; Matsumoto et al. 1997) find thatalthough there is a complex structure of periodic shockwave formation at the centre, the overall picture is ofa disc forming in the plane of rotation, with the cen-tral density and flatness increasing over time. Whilerotation does not seem to change the qualitative na-ture of the collapse, it may be important in determiningthe subsequent evolution of the central object. Rotationalso causes the ”angular momentum problem”, which is that, theoretically, rotating prestellar cores have muchmore angular momentum than a star could contain with-out breaking up (Prentice & ter Haar 1971). Proposedsolutions include the fragmentation of the core into mul-tiple protostellar objects (Matsumoto & Hanawa 2003)and the removal of angular momentum by magneticfields, known as magnetic braking (Basu & Mouschovias1994). Simulations of the fragmentation and collapse ofa magnetised filament (Nakamura et al. 1995; Tomisaka1996) show that, as with rotation, once a cloud has be-gun to collapse dynamically, magnetic forces are not suf-ficient to halt it. However, their conclusions depend onthe coupling between the magnetic field and the gas andthis is determined by the fractional ionization, which it-self is controlled by the chemical evolution during thecollapse and not included in these works.The density structure of prestellar cores can be in-ferred from observations of either the dust continuum,or line emission from molecules. Both methods requirea conversion factor to relate either the dust emissivityor the column density of the observed molecular speciesto the total gas density; these conversion factors intro-duce systematic uncertainties. Additionally, at typicaldistances observations of the central regions are limitedby resolution, while in the low-density outer regions thelower signal-to-noise is a further source of uncertainty.The differences in density structure between many pro-posed models of collapse are often not large enough forobservations to conclusively favour one over the oth-ers. However, the timescale for molecular abundancesin a core to reach equilibrium levels are comparable tothe timescales involved in gravitational collapse, whichmeans that the chemical composition of a cloud shouldbe sensitive to the hydrodynamical situation. Simula-tions of the chemical evolution of a core undergoingcollapse (Rawlings et al. 1992; Bergin & Langer 1997;Aikawa et al. 2001) have found that this is the case,with different hydrodynamical models giving significantdifferences in the abundances of some molecules. Thisraises the possibility that molecular abundances couldbe used as an observational test of theories of star for-mation. This is our motivation.However, there are difficulties in coupling detailedchemical models to a hydrodynamical code, especiallywhen magnetic fields are included, and simpler chemi-cal models cannot necessarily be relied on to give accu-rate molecular abundances. For example, Aikawa et al.(2001) and Rawlings et al. (1992) used analytical solu-tions for isothermal collapse, found by Larson (1969)and Shu (1977) respectively, but these are not necessar-ily applicable to real situations - the Larson (1969) solu-tion only agrees with simulations at small radii (Hunter hemical evolution of prestellar cores O, finding that quasi-equilibrium contrac-tion of a BE sphere was best able to reproduce observa-tions.In this paper we propose a different approach: weparametrize the results of hydrodynamical simulationsof collapsing prestellar cores to describe how the densitybehaves as a function of radius and time for differentmodels, and incorporate these parametrizations into agas-grain time dependent chemical model. Although lessaccurate than a simultaneous solution of both the hydro-dynamics and chemistry, this approach removes the needfor simplifications in either area, while also being muchless computationally expensive, and so enabling the ex-ploration of larger regions of parameter space than hasso far been feasible.This paper is laid out as follows: in Section 2, wedescribe our parametrizations of hydrodynamical simu-lations, and we discuss the chemical model into whichwe incorporate them; in Section 3 we present the resultsof our grid of models, showing how the abundances ofkey molecules are affected; in Section 4 we investigatewhether the inclusion of ambipolar diffusion affects starformation timescales; we briefly discuss our findings andconclude in Section 5. In Appendices A and B we list thefunctions used in our parametrization of the numericalsimulations. METHODOLOGY2.1.
Parametrization of numerical simulations
Empirical models were developed to reproduce theresults from four numerical simulations of collapsing prestellar cores: Aikawa et al. (2005) used a BE sphereas the initial configuration, with the density increasedby a factor of either 1.1 (model BES1) or 4 (BES4) totake the core out of equilibrium and instigate collapse;Nakamura et al. (1995) studied the collapse of a mag-netically supported filament when a density perturba-tion is applied (MS); and Fiedler & Mouschovias (1993)followed the evolution of a core as magnetic support isremoved through ambipolar diffusion (AD).Each of these studies produced the density profile(number density for BES1, BES4 and AD, mass den-sity for MS) of the core as a function of time during thecollapse. We extracted data from published plots of thedensity profiles at different times, as shown in Figure 1with data from Aikawa et al. (2005) as an example.For each density profile, a function, depending on po-sition in the core, reproducing its shape was found. Inthe 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 direc-tion was considered, as the core rapidly collapses into athin disc so that structure along the z -axis is less signif-icant. For the BES1, BES4 and AD models, a functionof the form n ( r ) = n rr ) a (1)with n , r and a free parameters determining the cen-tral density, the width of the central density peak, andthe slope of the profile, provided a good fit. Tafalla et al.(2002) used Eq. 1 to fit the observed density profiles ofprestellar cores. For the MS model, the radial profilescould not be approximated by Eq. 1. The equilibriumdensity of the filament is given by an equation of theform ρ ( r ) = ρ (cid:18) rr (cid:19) ! − a (2)which was adapted to reproduce the data by changingthe outer exponent a , so that the slope at large r is thesame as in the simulated data.The values of the free parameters were chosen to ap-proximate 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 alsoapproximated by a similar method. For example, thecentral density parameter’s time evolution was found tobe well reproduced bylog n ( t ) = A ( t − t ) a + B (3)for all simulations, where t is the simulation’s duration. a was chosen such that log n is approximately linearwith ( t − t ) a , and the coefficients A and B can then by Priestley et al. found by linear regression. Figure 2 shows the variationof the central density parameter, n , with time for thedata from Aikawa et al. (2005), along with the resultingapproximation to the time dependence of this parame-ter. Once all the parameters have been approximated inthis way, the density can be calculated as a function oftime and space, shown in Figure 3 compared with theoriginal data. Figures 4, 5 and 6 show the correspond-ing density profile approximations for the BES4, MSand AD collapses respectively. The equations used toapproximate the time dependence of the density profilesare given in Appendices A and B. The maximum dis-crepancies between the simulations and approximateddensities are 34%, 240%, 66% and 53% for the BES1,BES4, MS and AD cases respectively, while the averagediscrepancies are 10%, 26%, 16% and 16%. The large( > r . For the BES1 and BES4models, we determine the new radius of a parcel at eachtime step by calculating the mass interior to its initialradius at t = 0, M ( < r ), and finding the radius atthe given time which encloses the same mass. The MSand AD approximations are not spherically symmetric,so this approach would require knowledge of the densityat each point. Instead, we use the same methods as forthe density to approximate the radial velocity profiles,and use these to calculate the new parcel radius at eachtimestep. The gas density versus time, for parcels ofdiffering initial radii, in the BES1 collapse is shown inFigure 7. 2.2. Chemical modelling
The four density approximations were used as inputfor the chemical code UCL CHEM. The code is de-scribed in Viti et al. (2004) and references therein. Thiscode has since been made public (https://uclchem.github.io/)and is fully explained in Holdship et al. (2017). Herewe briefly summarize its characteristics: UCL CHEMis a time dependent gas-grain chemical model that cal-culates the abundances of atoms and molecules in thegas and dust in the interstellar medium as a functionof time under chemical and physical conditions set bythe user. The original version of the code uses free-fall collapse to determine the density from the diffusestate to the final density of the gas where the star isborn. Initial atomic elemental abundances are pro- − − − r / pc n / c m − . . . . Figure 1.
Density profiles taken from Aikawa et al. (2005)(solid lines), with the approximate profiles calculated usingEq. 1 (dashed lines). The labels indicate the time since col-lapse in 10 yr. . . . . . . . t / yr n / c m − Figure 2.
Central density n against t for the approxi-mations to the Aikawa et al. (2005) data (solid), with theapproximate fit to the time evolution (dashed). − − − r / pc n / c m − . . . . Figure 3.
As Figure 1, but with the parameters for Eq. 1calculated as a function of time. hemical evolution of prestellar cores -3 -2 -1 r / pc10 n / c m − Figure 4.
As Figure 3, for the BES4 collapse -4 -3 -2 -1 r / pc10 n / c m − Figure 5.
As Figure 3, for the MS collapse -5 -4 -3 -2 -1 r / pc10 n / c m − Figure 6.
As Figure 3, for the AD collapse vided to UCL CHEM which then self-consistently cal-culates gas-phase chemistry, as well as sticking on todust particles with subsequent surface processing. Forthe reaction network we used the UMIST 2012 network(McElroy et al. 2013), and freeze-out and grain surfacereactions as described in Holdship et al. (2017). t / Myr10 n / c m − Figure 7.
Parcel density versus time for the BES1 collapse,for different initial parcel radii.
We chemically model our four parametrized numeri-cal simulations: the collapse of an unstable (BES1) orhighly unstable (BES4) Bonnor-Ebert sphere, collapseagainst magnetic support (MS) and collapse resultingfrom ambipolar diffusion (AD). A grid of models wasrun for each case to investigate the effects of chang-ing other input parameters: the cosmic ray ionisationrate ζ (in units of ζ = 1 . × − s − ), metallicity Z ,and the desorption efficiency parameters ǫ , φ and Y UV ,corresponding to the number of molecules desorbed perH molecule formed, per cosmic ray impact and per UVphoton (produced by cosmic rays) absorbed respectively(Roberts et al. 2007). For our fiducial desorption effi-ciencies, H formation is the dominant desorption mech-anism, as found by Roberts et al. (2007). The valuesadopted for each model are given in Table 1. Each modelwas run once for each of the density approximations. Weassume an external radiation field of 1 Habing, and anexternal extinction at the core boundary of 3 mag, thevalue used by Aikawa et al. (2005). The extinction fromthe core itself is calculated by integrating the densityprofile to the boundary, (0 . . .
75 pc for the ADcase) before the onset of collapse, giving maximum ex-tinctions (at the centremost parcel) of 6.4 (BES1), 15.2(BES4), 7.8 (MS) and 3.3 (AD) mag. We used 13 gasparcels (14 for AD) at initial radii spaced to cover theentire range of the cores, but with an emphasis on themore rapidly evolving central regions.The initial central number densities of the models are n = 2 . × cm − (BES1 and MS), 8 × cm − (BES4)and 300 cm − (AD). The MS equations are in termsof dimensionless variables, which can be converted intophysical values by choosing the intial central density, ρ c ,and the isothermal sound speed, c s = kTµm H . The gas wasassumed to be at a temperature of 10 K and composed Priestley et al.
Table 1.
Model input parametersModel ζ / ζ Z/Z ⊙ ǫ φ Y UV A 1 1 0.01 10 Table 2.
Elemental abundancesElement Abundance Element AbundanceH 1 . . × − He 0 .
085 S 1 . × − C 2 . × − Si 3 . × − O 4 . × − Cl 3 . × − entirely of molecular hydrogen for the purposes of calcu-lating the mean molecular mass, giving c s = 203 ms − ,and we set ρ c = 3 . × − gcm − to give an initialcentral number density equal to the BES1 model. Theinitial number density at each point is then given by therelevant equation for the density profile, and we allowthe chemistry to evolve for 1 Myr before the onset ofcollapse at this density. The models were run until thecentral density reached 10 cm − . The elemental abun-dances relative to H for our standard model, the solarvalues given by Asplund et al. (2009), are listed in Table2 - the abundances for elements other than H and Hein models with varying metallicity are multiplied by thevalue of Z . RESULTS3.1.
Molecular abundances across the different modesof collapse
Figure 8 shows the density profiles of the four collapsemodes at the point when the central number density, n , reaches 2 × cm − . The BES1 and MS profilesdecrease more rapidly with distance than for the BES4and AD approximations, which have similar densitiesup to the end of the BES4 core at 0 . ∼ yr) than for AD ( ∼ yr), so thechemical evolution differs significantly despite the gasdensity being the same. The BES1 and MS modes bothreach n = 2 × cm − after ∼ yr, and as suchthe chemical evolution is similar. The densities at large( & . n = 2 × cm − . In all cases,the CO abundance increases from a central minimum,where freeze-out has depleted most of the molecule ontograin surfaces, to a maximum value of ∼ − . TheBES1, BES4 and MS approximations all behave simi-larly in reaching the maximum, although for BES4 andMS the radius at which this occurs is larger due to thehigher gas densities increasing the effect of freeze-out.The CO abundance in the AD collapse increases muchmore slowly, only reaching 10 − at the edge of the core,whereas in the MS case the abundance begins to fallagain towards the edge. This is due to the higher den-sities in the outer regions for the AD collapse, as mag-netic support can still prevent material here from col-lapsing, unlike in the other three cases. The effect ofthe longer collapse duration is also apparent, as the ADabundances are far lower than the BES4 ones at com-parable radii, despite the densities being very similar.The NH abundance also increases from the centre to amaximum, before falling with radius in all models. Thedecline is much more gradual for the AD case than theother collapse modes, leading to an order of magnitudedifference with the MS model by r = 0 . + shows similar behaviour, while for HCN the AD col-lapse mode produces a nearly constant abundance afterreaching a value of ∼ − , in contrast to the others.3.2. Cosmic ray ionization rate
Raising the cosmic ray ionization rate increases theabundances of molecules in the centre of the core, re-gardless of the collapse mode, as the rate of desorptionfrom grain surfaces is increased. Figure 10 shows theCO abundances for models A, B1 and B2, for the BES1and AD density approximations, at a central density n = 2 × cm − . For the BES1 collapse, the abun-dance at larger radii is mostly unaffected, whereas forAD the A and B2 models show noticeably different be-haviour, with the CO abundance an order of magnitudelower at the edge of the core for the B2 model due tothe higher cosmic ray dissociation rate. Figure 11 showsthe HCO + abundance for the same models. The central hemical evolution of prestellar cores . . . . . . . . . r / pc n / c m − Figure 8.
Density profiles of the BES1 (solid black), BES4(dashed black), MS (blue) and AD (red) approximations, ata central number density of n = 2 × cm − . abundances are again enhanced for models with higherionization rates, but whereas for the AD collapse modethe abundance decreases towards the edge as with CO,for BES1 the abundance is higher throughout the cloud.3.3. Metallicity
Changing the metallicity of the core usually resultsin a corresponding change in the molecular abundances,due to the different availability of atoms to form themolecules. However, some molecules are much less af-fected than others. Figure 12 shows the abundances ofCO and HCN for the BES1 collapse at a central density n = 2 × cm − , for models A, C1 and C2. Whereasthe CO abundance scales nearly linearly with the metal-licity, the HCN abundance is virtually unchanged be-tween models. The main formation and destructionreactions for HCN both involve H + , for which the abun-dance increases with decreasing metallicity (and viceversa), at least partially counteracting the effect of thechanging elemental abundances on the HCN abundance.3.4. Desorption efficiencies
As with the cosmic ray ionization rate, increasingthe desorption efficiences increases the molecular abun-dances, particularly in the denser central regions wheremore freeze-out has taken place. Figure 13 shows theCO abundances for models A, D1 and D2, for the MSand AD density approximations, where the H forma-tion desorption efficiency ǫ has been modified. Whilethe abundances in the central regions are affected sim-ilarly for both collapse modes, at larger radii the ef-fect is negligible for the MS collapse, whereas the COabundance reaches 10 − much more rapidly for AD -for model D2, the abundance profile looks much moresimilar 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 changein the CO abundance beyond 0 . φ has very little effect onthe abundance of any molecule, despite a factor of 100difference between models E1 and E2. The cosmic rayinduced photodesorption efficiency Y UV , however, doesaffect molecular abundances, in particular having a sig-nificant effect on the abundance of NH , which is notgreatly affected by variation of the other parameters in-vestigated. STAR FORMATION EFFICIENCIESOnly one of our modes of collapse, AD, includes the ef-fect of ambipolar diffusion. However, all prestellar coresare expected to be magnetised, and therefore ambipolardiffusion could be important for all collapse models. Anestimate of the timescale on which ambipolar diffusionoccurs is given by t amb = 4 × ( x i / − ) yr (4)(Mouschovias 1979; Hartquist & Williams 1989). If t amb is smaller than the free-fall timescale, magnetic pressureis unlikely to impede gravitational collapse, while if it islarger the impeding effects may be significant.Banerji et al. (2009) showed that the ambipolar diffu-sion timescale becomes very large as the fractional ion-ization increases and the magnetic field is strongly cou-pled to the collapsing core, which, in some cases, mayhalt the collapse and hence the formation of the star.We calculated the ambipolar diffusion timescale at thecentre of the core using Eq. 4 at the beginning of thecollapse, and at central densities of 10 and 10 cm − ,for the BES1, BES4 and MS approximations, for differ-ing metallicities and cosmic ray ionization rates, usingionization fractions calculated in our chemical simula-tions. We compared these timescales with the time ittakes the gas to reach the final density (10 cm − ). Ta-ble 3 shows the collapse duration and t amb for this gridof models.For all models considered, the value of t amb at the fi-nal density, n H = 10 cm − , is significantly lower thanthe collapse duration, while the initial values are com-parable to or larger than the collapse time, particularlyfor the BES4 collapse, and for the models with increasedcosmic ray ionization rates. At 10 cm − , the B1 andB2 models have t amb larger than the collapse time for theBES4 case. Increasing ζ leads to larger values of t amb ,as the ionization, and therefore the coupling to the neu-tral gas, is increased. Metallicity has very little effectat higher densities, but the initial t amb varies with themetallicity, as more readily-ionized atoms such as carbonare present. The BES1 and MS models have lower ini- Priestley et al. . . . . . . . . . r / pc − − − − − − n / n H (a) CO . . . . . . . . . r / pc − − − n / n H (b) NH . . . . . . . . . r / pc − − − − − n / n H (c) HCO + . . . . . . . . . r / pc − − − − − n / n H (d) HCN Figure 9.
Abundances of CO, NH , HCO + and HCN at a central density n = 2 × cm − for model A, using the BES1(solid black), BES4 (dashed black), MS (blue) and AD (red) density approximations. tial values of t amb / t collapse than the BES4 models, sug-gesting that the faster collapse should be impeded morestrongly by the coupling of gas to magnetic fields. Theseresults emphasize that ambipolar diffusion is importantfor diffuse material, where magnetic fields are likely toimpede collapse, but once denser clumps have formedmagnetic support will be removed too rapidly to affectthe subsequent evolution. DISCUSSION & CONCLUSIONSThe results of our chemical modelling show that itis difficult to disentangle the effect on molecular abun-dances of different collapse modes from that of varyingthe input parameters, which are not necessarily knowna priori. Drawing information about the dynamics of star formation from molecular abundances therefore re-quires a full investigation of parameter space, somethingwhich would be extremely time-consuming using com-bined hydrodynamical-chemical modelling. Our gridof models, although not large enough to draw robustconclusions about individual objects, does allow us tocompare results with the general properties of prestellarcores, and determine whether particular collapse mod-els or regions of parameter space are in conflict withobservation.Assuming the values from model A (see Table 1), alldensity approximations predict CO abundances awayfrom the core centre in agreement with observed val-ues of 10 − − − in starless cores (Caselli et al. hemical evolution of prestellar cores -8 -7 -6 -5 -4 -3 n / n H (a) BES1 -8 -7 -6 -5 -4 -3 n / n H (b) AD Figure 10.
Abundance of CO versus radius at a central density n = 2 × cm − for models A (solid line), B1 (dashed line)and B2 (dotted line), using the BES1 (left) and AD (right) density approximations. -11 -10 -9 -8 n / n H (a) BES1 -11 -10 -9 -8 n / n H (b) AD Figure 11.
Abundance of HCO + versus radius at a central density n = 2 × cm − for models A (solid line), B1 (dashedline) and B2 (dotted line), using the BES1 (left) and AD (right) density approximations. (1999); Frau et al. (2012), assuming O / O ≈ − and O / O ≈ − ). However, the AD collapse onlyreaches these values at r & . < . formation and cosmic-ray induced photodesorption(D2 and F2) predict CO abundances of ∼ − at aradius of 0 . + , and the subsequent slow or negligible decline0 Priestley et al. -8 -7 -6 -5 -4 -3 n / n H (a) CO -11 -10 -9 -8 -7 n / n H (b) HCN Figure 12.
Abundance of CO (left) and HCN (right) versus radius at a central density n = 2 × cm − for models A (solidline), C1 (dashed line) and C2 (dotted line), using the BES1 density approximation. -8 -7 -6 -5 -4 -3 n / n H (a) MS -8 -7 -6 -5 -4 -3 n / n H (b) AD Figure 13.
Abundance of CO versus radius at a central density n = 2 × cm − for models A (solid line), D1 (dashed line)and D2 (dotted line), using the MS (left) and AD (right) density approximations. beyond the peak value, are qualitatively different to thesituation with the initially unstable collapse modes, sug-gesting spatially resolved molecular observations couldbe used to discriminate between them.All density approximations predict peak NH abun-dances of ∼ − , consistent with observed values inprestellar cores (Tafalla et al. 2002; Johnstone et al.2010). At n = 2 × cm − , N H + abundances,shown in Figure 14, are intermediate between the valuesfound by Frau et al. (2012) ( ∼ − ) and Tafalla et al.(2002) and Johnstone et al. (2010) ( ∼ − ). However, as with CO, for the AD approximation the abundanceswithin the reported radii of the cores are much lowerthan the observed values. Ambipolar diffusion simula-tions including chemistry by Tassis et al. (2012) showsimilar behaviour, with the inclusion of magnetic effectsleading to a more extended depleted region than in un-magnetised models, suggesting that this is a genuinefeature of collapse under ambipolar diffusion, ratherthan being due to our approximation of the densityevolution. hemical evolution of prestellar cores Table 3.
Collapse duration and ambipolar diffusion timescales at increasing density for BES1, BES4 and MS models withvarying ζ and Z. t amb /10 yrModel Z/Z ⊙ ζ / ζ t collapse /10 yr Initial n H = 10 cm − n H = 10 cm − BES1 A 1 . . .
173 0 .
80 0 .
09 0 . . . .
173 2 .
51 0 .
14 0 . . . .
173 4 .
70 0 .
19 0 . . . .
173 0 .
62 0 .
09 0 . . . .
173 0 .
95 0 .
09 0 . . . .
184 0 .
22 0 .
08 0 . . . .
184 0 .
72 0 .
12 0 . . . .
184 1 .
35 0 .
18 0 . . . .
184 0 .
22 0 .
08 0 . . . .
184 0 .
22 0 .
08 0 . . . .
393 0 .
73 0 .
09 0 . . . .
393 2 .
39 0 .
14 0 . . . .
393 4 .
49 0 .
19 0 . . . .
393 0 .
58 0 .
09 0 . . . .
393 0 .
85 0 .
09 0 . Lee et al. (2004) calculated the chemical evolution ofa quasistatically contracting BE sphere over 10 yr, find-ing similar abundance profiles to our BES1 approxima-tion for CO, NH and HCO + , although for N H + wefind much more depletion in the core centres, and higherHCN abundances overall. Given that our MS and BES4approximations also give similar results to BES1, thissuggests that the collapse timescale, rather than the spe-cific details of the density evolution, are more importantfor the chemical evolution, at least for some observation-ally important molecules such as CO.Our BES1 and BES4 approximations are based onhydrodynamical simulations presented in Aikawa et al.(2005), who modelled the chemical evolution of thesemodels self-consistently, providing a test of the approx-imations’ accuracy. Comparing to their results, we findgood agreement (of the same order of magnitude) be-tween the predicted peak molecular abundances of bothapproaches. However, the variation with radius differs -our approximations generally predict lower abundancesat the core centres, especially for the BES4 collapse,where our results predict lower CO abundances to theBES1 case at the same r , whereas Aikawa et al. (2005)find the opposite. We attribute this to the differing ini-tial conditions - Aikawa et al. (2005) assume the gas isentirely atomic prior to collapse, whereas we allow theabundances to evolve for 10 yr at the initial density be- fore the onset of collapse, leading to significant freeze-out onto grains occuring in the denser central regions.We conclude that despite the dependence of molec-ular abundances on the various parameters mentionedabove, molecular observations can still be useful for dis-criminating between different models of collapse. Whilemodels which begin from an initially gravitationally un-stable state predict similar abundances and radial varia-tions, ambipolar diffusion produces qualitatively differ-ent abundance profiles for many observationally impor-tant molecules, which appear to be in conflict with ob-servations of prestellar cores, although varying the inputparameters may be able to reduce this discrepancy. Amore exhaustive investigation of parameter space, com-bined with observations of multiple species from thesame source, could be used to draw much stronger con-clusions on the nature of core collapse, as well as pro-viding constraints on the values of the input parameterswhich are currently assumed ad hoc. This sort of investi-gation would be extremely time consuming, if not impos-sible, using a coupled hydrodynamical-chemical system,even without the additional complications of magneticfields. The results we have presented here demonstratethat these large grids of models are now feasible usingour method of parametrizing the dynamics, while stillproviding a reasonable level of accuracy compared tofull simulations.2 Priestley et al. . . . . . . . . . r / pc − − − − n / n H Figure 14. N H + abundance at a central density n =2 × cm − for model A, using the BES1 (solid black),BES4 (dashed black), MS (blue) and AD (red) density ap-proximations. FP is supported by the Perren fund and the IMPACTfund. We would like to thank Prof. Fred Adams forhis helpful contribution, and the referee, whose com-ments helped to produce a greatly improved version ofthe original paper.
Software:
UCL CHEM (Viti et al. 2004; Holdship et al.2017)APPENDIX A. DENSITY APPROXIMATIONSThe BES1, BES4 and AD collapse density profiles were approximated with the function n ( r ) = n ( t )1 + ( rr ( t ) ) a ( t ) (A1)whereas the MS collapse required a different functional form, ρ ( r ) = ρ ( t ) (cid:18) rr ( t ) (cid:19) ! − a ( t ) (A2)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 bylog n ( t ) = 61 . (cid:0) . × − t (cid:1) − . − . r ( t ) = − . (cid:0) . × − t (cid:1) − . + 28 .
93 (A4) a ( t ) = 2 . n is in cm − , r is in AU and t is in years. A.2. BES4
The time-dependent parameters are given bylog n ( t ) = 68 . (cid:0) . × − t (cid:1) − . − . r ( t ) = − . (cid:0) . × − t (cid:1) − . + 38 . a ( t ) = 1 . . − t/ ) (A8)where n is in cm − , r is in AU and t is in years. hemical evolution of prestellar cores MS The time-dependent parameters are given bylog ρ ( t ) = 3 .
54 (5 . − t ) − . − .
73 (A9)log r ( t ) = − .
34 (5 . − t ) − . + 1 .
47 (A10) a ( t ) = 2 . − . (cid:18) t . (cid:19) (A11)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 c s √ πGρ c and t in units of (2 πGρ c ) − . . A.4. AD The time-dependent parameters are given bylog n ( t ) = log (2 + 1 . (cid:0) t − (cid:1) ) + 3 t < . . . − t ) − . − . t ≥ . r ( t ) = − .
57 (16 . − t ) − . + 1 .
85 (A14) a ( t ) = 2 . − . (cid:18) t . (cid:19) (A15)where n is in cm − , r is in 0 .
75 pc and t is in Myr. B. VELOCITY PROFILESFor the MS and AD collapses, the radial velocity profiles were also approximated in order to determine the inwardsmovement of the gas parcels. B.1. MS The radial velocity profile is given by v r ( r ) = v min ( t ) "(cid:18) r ′ ( t ) r min ( t ) (cid:19) − (B16)for r < r min and v r ( r ) = v min ( t ) [exp( − a ( t ) r ′ ( t )) − − a ( t ) r ′ ( t ))] (B17)for r ≥ r min , where r ′ ( t ) = r − r min . The time evolution is given by − . t + 7 . t < .
95 (B18) r min ( t ) = − . t + 16 . t < .
33 (B19) −
22 log t + 37 . t ≥ .
33 (B20)0 . t t < .
95 (B21) v min ( t ) = 5 . t − . t < .
33 (B22)18 . t − . t ≥ .
33 (B23)0 . t + 0 . t < .
95 (B24) a ( t ) = 0 .
695 log t − . t < .
33 (B25)2 .
69 log t − t ≥ .
33 (B26)where v min is in units of c s , r min in units of c s √ πGρ c and t in units of (2 πGρ c ) − . .4 Priestley et al.
B.2. AD The radial velocity profile is given by v min ( t ) (cid:20)(cid:16) r ′ ( t ) r min ( t ) (cid:17) − (cid:21) r < r min (B27) v r ( r ) = ( v min ( t ) − v mid ( t )) (cid:16) r ′ ( t )0 . − r min ( t ) (cid:17) . − v min ( t ) r < . v mid ( t )( r − r ≥ . − . t + 0 . t ≤ . r min ( t ) = − . t − .
2) + 0 . t ≤ . − . t − .
1) + 0 . t > . v min ( t ) = 3 . . − t ) − . − . v mid ( t ) = 0 . t t ≤ . . t − .
2) + 1 . t > . r min is in units of 0 .
75 pc, v min and v mid are in 10 − km s − and t is in Myr.REFERENCES Aikawa, Y., Herbst, E., Roberts, H., & Caselli, P. 2005,ApJ, 620, 330Aikawa, Y., Ohashi, N., Inutsuka, S.-i., Herbst, E., &Takakuwa, S. 2001, ApJ, 552, 639Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009,ARA&A, 47, 481Banerji, M., Viti, S., Williams, D. A., & Rawlings, J. M. C.2009, ApJ, 692, 283Basu, S., & Mouschovias, T. C. 1994, ApJ, 432, 720Bergin, E. A., & Langer, W. D. 1997, ApJ, 486, 316Bonnor, W. B. 1956, MNRAS, 116, 351Caselli, P., Walmsley, C. M., Tafalla, M., Dore, L., &Myers, P. C. 1999, ApJL, 523, L165Cohen, M., & Kuhi, L. V. 1979, ApJL, 227, L105Fiedler, R. A., & Mouschovias, T. C. 1993, ApJ, 415, 680Frau, P., Girart, J. M., Beltr´an, M. T., et al. 2012, ApJ,759, 3Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers,P. C. 1993, ApJ, 406, 528Hartquist, T. W., & Williams, D. A. 1989, MNRAS, 241,417Hincelin, U., Commer¸con, B., Wakelam, V., et al. 2016,ApJ, 822, 12 Hincelin, U., Wakelam, V., Commer¸con, B., Hersant, F., &Guilloteau, S. 2013, ApJ, 775, 44Holdship, J., Viti, S., Jim´enez-Serra, I., Makrymallis, A., &Priestley, F. 2017, AJ, 154, 38Hunter, C. 1977, ApJ, 218, 834Johnstone, D., Rosolowsky, E., Tafalla, M., & Kirk, H.2010, ApJ, 711, 655Kandori, R., Nakajima, Y., Tamura, M., et al. 2005, AJ,130, 2166Keto, E., Caselli, P., & Rawlings, J. 2015, MNRAS, 446,3731Kirk, J. M., Ward-Thompson, D., & Andr´e, P. 2005,MNRAS, 360, 1506Larson, R. B. 1969, MNRAS, 145, 271Lee, J.-E., Bergin, E. A., & Evans, II, N. J. 2004, ApJ, 617,360Li, Z.-Y., Shematovich, V. I., Wiebe, D. S., & Shustov,B. M. 2002, ApJ, 569, 792Matsumoto, T., & Hanawa, T. 2003, ApJ, 595, 913Matsumoto, T., Hanawa, T., & Nakamura, F. 1997, ApJ,478, 569McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A,550, A36Motte, F., Andre, P., & Neri, R. 1998, A&A, 336, 150 hemical evolution of prestellar cores15