EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 652 (August 2021) A&A, 652 (2021) A71 Full HTML
Free Access
Issue
A&A
Volume 652, August 2021
Article Number A71
Number of page(s) 19
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202140469
Published online 12 August 2021

© ESO 2021

1 Introduction

Although high-mass stars (i.e. M ≳ 8–10 M⊙) are much rarer than their less massive counterparts, they have a critical effect on the physico-chemical characteristics of the interstellar medium (ISM). They also play an important role in the evolution of the host galaxies (e.g. Kennicutt 2005). Massive stars are responsible for the production of large amounts of the α-elements involved in the formation of complex molecules (Woosley & Weaver 1995; Kobayashi et al. 2020) that are created during the final stages of their life-cycle through core-collapse supernovae (e.g. Smartt 2009). They also dominate the energy budget in their immediate surroundings by stirring, heating, and ionising the gas, and they affect the chemical evolution of the ISM as well as the star and planet formation process (e.g. Elmegreen 1998; Bally et al. 2005; Adams 2010). However, observing massive stars is challenging as they evolve quickly, and the initial stages of their formation process take place when their progenitors are still embedded in the parental clump (Zinnecker & Yorke 2007; Motte et al. 2018, for reviews on this topic). In spite of the efforts madein recent years to understand the formation of massive stars, some fundamental questions still remain unanswered.One of these is related to the timescales of the various (evolutionary) phases of their formation process, which would provide crucial information to distinguish among competing star formation theories (e.g. McKee & Tan 2002; Mouschovias et al. 2006; Hartmann et al. 2012; Tigé et al. 2017; Padoan et al. 2020) and to reach a comprehensive view of the chemical evolution of the gas that resides in the parental clumps of massive stars.

Timescales are usually derived by applying statistical approaches (e.g. Wood & Churchwell 1989; Davies et al. 2011; Mottram et al. 2011; Battersby et al. 2017; Tigé et al. 2017; see also Motte et al. 2018). Typical values for the total massive star formation timescale range from ~1 to 5 times 105 yr. The statistical lifetime of a phase is consequently obtained as the fraction of the total time equal to the number of objects in that phase with respect to the total number of objects. These methods depend on the assumption of a total time for the star formation process, however, that is constrained by simulations (e.g. Davies et al. 2011) or derived with respect to the known age of OB stars (e.g. Wood & Churchwell 1989; Motte et al. 2007; Russeil et al. 2010; Mottram et al. 2011; Tigé et al. 2017).

An alternative way to quantify the timescales of the various steps is to study the effects that massive young stellar objects (YSOs) have on their immediate surroundings. The chemical evolution is affected by changes in density and temperatureinduced by the forming stars. During the so-called warm-up phase (e.g. Viti & Williams 1999; Viti et al. 2001; Garrod et al. 2008), molecules formed on the surface of the dust grains are rapidly released into the gas phase, increasing their observed abundance by several orders of magnitude (Viti et al. 2004; Garrod & Widicus Weaver 2013; Choudhury et al. 2015). Thus, the molecular composition of the material around YSOs contains information about the physical conditions of their past (van der Tak 2005). A well-known example are the rich molecular spectral features produced by hot molecular cores (HMCs), which are compact regions around YSOs (≲ 0.1 pc radius; Kurtz et al. 2000; Cesaroni 2005).

Chemical tracers that show a relation between their observed abundances and the different phases of the star formation process are commonly called chemical clocks (e.g. Fontani et al. 2007; Beuther et al. 2009; Hoq et al. 2013; Giannetti et al. 2019; Urquhart et al. 2019; Sabatini et al. 2020, as some recent examples; see also van Dishoeck & Blake 1998 for a review). Through the comparison of their observed abundances obtained from large samples of massive clumps at different evolutionary stages and those predicted by accurate time-dependent chemical models, it is possible to quantify the timescales of each stage. An attempt to adopt this type of analysis was carried out by Gerner et al. (2014) based on the chemical properties of a sample of 59 massive clumps in different evolutionary stages, and modelling the warm-up phase alone. The radial profiles of temperature and density were assumed to be static with time and were described as power laws with the average properties of each evolutionary class (see their Sect. 5.1.1 for more details). The chemistry was then evolved in time until they found the best match with the average observed abundances. They interpreted this time as the typical age of each class. Based on these results, Gieser et al. (2019) studied the chemical history of the HMC VLA 3 in the high-mass star-forming region AFGL 2591 to assess the effects of different initial chemical conditions on the model results and on the derived evolutionary timescales. They showed that the latter are sensitive to the model assumptions. More recently, Gieser et al. (2021) used the same method to investigate 22 cores identified in 18 high-mass star-forming regions with ALMA.

An alternative to these approaches is to adopt a time-dependent analytical function of the temperature to simulate the gradual warming-up of the clumps in a range of typical observed values (i.e. ~10–300 K; e.g. Viti & Williams 1999; Garrod & Herbst 2006). Awad et al. (2010) proposed an improved model of the warm-up phase applied to the low-mass regime, where the temperature evolution was based on the radiative transfer (RT) calculation of Nomura & Millar (2004). More recently, Bonfand et al. (2019) described the evolution of the physical parameters of the high-mass star-forming region Sgr B2 with 3D-RT simulations. They followed the trajectory of a parcel of gas under free-fall collapse and modelled its chemistry by updating the temperature and density of the gas. In this paper, we have developed a model similar to Bonfand et al. (2019) by employing a time-dependent description of the thermal evolution of massive clumps during the warm-up phase, with the aim to propose a new method for deriving the evolutionary timescales of the massive star formation process.

This paper has the following structure: in Sect. 2 we present the reference sample of massive clumps, the evolutionary sequence, and the chemical tracers we employed to derive the duration of each phase. The model and chemical network are described in Sect. 3. In Sect. 4 we present our results, the post-processing procedure, and the derivation of the durations of each evolutionary phase. In Sect. 5 we discuss our estimates and how the selected tracers are reproduced by the models, providing information on the reliability of the chemical clocks. Finally, in Sect. 6 we summarise our conclusions.

2 Reference sample and selected tracers

The survey we used to define the reference sample and the evolutionary stages of the massive star formation process is the APEX1 Telescope Large Area Survey of the Galaxy (ATLASGAL; Schuller et al. 2009). ATLASGAL offers a complete view of the high-mass star-forming regions in the inner Galaxy at 870 μm and provides the ideal basis for detailed studies of a large number of massive clumps in different evolutionary stages, with the unprecedented angular resolution of ~19′′ (Contreras et al. 2013; Csengeri et al. 2014; Urquhart et al. 2013). Estimates of luminosities, masses, kinematic distances, and dust temperature (Urquhart et al. 2014, 2018 and Wienen et al. 2015) have been derived for more than ~104 dust clumps.

Using the unbiased nature of ATLASGAL, the ATLASGAL-TOP100 sample (hereafter TOP100) has been defined as a flux-limited sample of high-column density clumps, selected with additional infrared (IR) criteria to include sources that potentially cover the whole spectrum of ages (see Giannetti et al. 2014). The number of sources included in the TOP100 was slightly refined by König et al. (2017), and to date, it includes 111 clumps. Combining the IR properties of the massive clumps in the TOP100 sample with radio-continuum measurements at 5-9 GHz, the objects were divided into the following four evolutionary classes: (1) 70 μm weak stage (70w; i.e. quiescent), constituted by sources undetected at 24 μm and showing no clear compact emission at 70 μm (or are seen in absorption at this wavelength). This stage represents the earliest phase of massive-star formation and potentially includes starless or prestellar cores. (2) Mid-IR weak stage (IRw; protostellar), composed of compact and visible sources at 70 μm that are still undetected at 24 μm or are associated with a weak IR source (<2.6 Jy fluxes). The clumps in this evolutionary stage are young and likely dominated by cold gas. (3) Mid-IR bright stage(IRb; high-mass YSOs), made by the sources that show a strong mid-IR emission at 8 and 24 μm, which was interpreted to be causedby the progressive dust heating induced by the forming (proto-)star(s). (4) HII regions (HII), where the sources are bright at 70 and 24 μm, and detected in radio continuum at 5-9 GHz.

In the sample, the clump masses range from ~18 to ~5 × 104 M⊙. The bolometric luminosities are between ~60 and ~4 × 106 L⊙ and correspond to a range in luminosity-to-mass ratio (L∕M) of ~0.2–350 L⊙/M⊙. All the TOP100 clumps have the potential of forming high-mass stars and show no bias in terms of distance or mass between the evolutionary classes (König et al. 2017).

Four molecular tracers were studied with the APEX telescope in the TOP100 sample, revealing their potential as chemical clocks: formaldehyde (H2CO; Tang et al. 2018), methyl acetylene, acetonitrile, and methanol (CH3CCH, CH3CN and CH3OH, respectively; Giannetti et al. 2017). The detection rates, abundances, and excitation temperatures derived from each molecular species increase with the L∕M of the clumps, which is known to trace the evolution of the star formation process (Saraceno et al. 1996; Molinari et al. 2008; Urquhart et al. 2018). With respect to the line-of-sight (LOS) and beam-averaged column densities reported by Giannetti et al. (2017) and Tang et al. (2018), we calibrated our method to estimate the duration of the four evolutionary phases defined in the TOP100 sample. We propose a general pipeline of comparisons of models and observations that can be expanded with additional tracers in future follow-ups.

thumbnail Fig. 1

Sketch of the physical model employed in this work. Left panel: collapse phase (see Sect. 3.1.1) solved in a single-zone approximation. The chemical output of phase I is used as initial condition of phase II. Central panel: warm-up phase solved in a 1D approximation from 1 to 105 AU (see Sect. 3.1.2). Right panel: post-processing applied to compare the final chemical outputs of phase II with the column densities observed in the TOP100 (see Sects. 4.1 and 4.2). The physical parameters assumed in phases I and II are summarised in Table 1. The time evolution of the model is indicated by the dotted blue arrows.

3 Methods

In this section we present the model we developed to describe the time evolution of the abundances of the selected molecular tracers in the TOP100 survey (see Sect. 2). We describe the physical and the chemical model separately in Sects. 3.1 and 3.2, respectively.

3.1 Physical model

Our physical model comprises two distinct stages that are sketched in Fig. 1. The first, namely the collapse (left panel of Fig. 1), describes the density evolution of an isothermal clump with a single-zone approximation, exploring different conditions as in Viti & Williams (1999) and Garrod & Herbst (2006), for example. In the second phase, the warm-up (central panel of Fig. 1), the gas and dust temperatures (assumed to be in equilibrium) evolve as a function of time, driven by the luminosity of the central forming protostar. The temperature profiles are computed with the radiative transfer code MOCASSIN (see Ercolano et al. 2003, 2005), and the mass distribution of the clump is described by a static gas radial density profile (e.g. Tafalla et al. 2004).

In the following sections we describe the details of each phase. The specific details of the RT calculations will be presented in a forthcoming paper (Grassi et al., in prep.).

3.1.1 Phase I: Isothermal collapse

The first phase of the model simulates a semi-analytical single-zone isothermal collapse (see Spitzer 1978; Brown et al. 1988; Viti & Williams 1999). The gas number density, nH, at the centre of the clump, evolves with time as (1)

where nH,0 is the initial central gas number density, b is a factor which aims to mimic a slower collapse compared to the ideal free-fall time (i.e. b = 1), , G is the gravitational constant, and mH is the hydrogen mass. Equation (1) is obtained assuming the conservation of mass during an isothermal contraction of a spherical clump (Spitzer 1978).

We set the initial gas number density to the typical values of a clump, nH,0 = n(H) + 2n(H2) = 104 cm−3 (Bergin & Tafalla 2007; Gerner et al. 2014). The collapse assumes a constant temperature of 15 K (for gas and dust temperatures; e.g. König et al. 2017) and a visual extinction Av = 10 mag (e.g. Semenov et al. 2010; Reboussin et al. 2014). The cosmic-ray ionisation rate of hydrogen molecules was set to ζ2 = 5 × 10−17 s−1, as observationally constrained by van der Tak & van Dishoeck (2000) in high-mass star-forming regions, and in agreement with the results reported by Sabatini et al. (2020). The specific density of the dust grains was ρ0 = 3 g cm−3, typical of silicates (e.g. Draine & Lee 1984), the dust-to-gas ratio , and we used a constant grain size μm. The gas mean molecular weight was μ = 2.4.

To simulate various physical conditions for the environment in which the seeds of massive protostars are formed, the collapse was stopped at different final densities nH(r0) (see Fig. 1) that were observationally constrained (Table 1; e.g. Mueller et al. 2002; Sabatini et al. 2019)and were then used as central densities during the warm-up phase.

Table 1

Parameter space explored in our models; i.e. phase I (Sect. 3.1.1) and phase II (Sect. 3.1.2).

3.1.2 Phase II: Warm-up

The second phase of the model, sketched in the central panel of Fig. 1, simulates the warm-up induced by a protostar at the centre of a spherical clump by using a 1D approximation of 100 logarithmic radial steps from 1 to 105 AU (i.e. up to a typical clump size of ~0.5 pc; e.g. Motte et al. 2018).

The mass distribution of the core is described by the gas radial density profile (e.g. Tafalla et al. 2002, 2004), (2)

where nH(r0) is the central number density at r0, at which phase I is stopped (see second row in Table 1). This ensures continuity between the two physical phases of the model. In Eq. (2), Rc is the core radius, that is, the radius at which nH(Rc) = 0.5 nH(r0). The slope 5/2 agrees with the observations in massive star-forming regions (e.g. Mueller et al. 2002; Schneider et al. 2015), while similar profiles were successfully applied to model the H2 distribution in clumps and filamentary structures (e.g. Beuther et al. 2002a; Arzoumanian et al. 2011; André et al. 2016; Sabatini et al. 2019). Alternatively, power-law profiles can be assumed to model the H2 distribution especially on the clumps scale (e.g. Gieser et al. 2021). This provides minor differences in the final clump mass (a factor of ≲2 smaller for a slope of 5/2).

The physical parameters of phase II (i.e. ζ2, , μ, and ) were assumed to be the same as those of the collapse phase, while the visual extinction at each radius r, was Av (r) = N(H, r)∕(2 × 1021 cm−2), where N(H, r) is the column density of the hydrogen nuclei, obtained by integrating Eq. (2) from the edge of the clump (e.g. Tielens 2010; Zhu et al. 2017). We took the final abundances of phase I as initial chemical conditions for each radial grid point of phase II, rescaled by the density profile. The chemistry was evolved in time, assuming that the individual cells are independent throughout the entire time evolution.

To model the thermal evolution of the clump during the warm-up phase, we generated a set of models with different density distributions and protostar luminosities by employing the Monte Carlo radiative transfer code MOCASSIN. The models provide an interpolable look-up table Td (r, t2) = fb[Rc, nH (r0), L* (t2), r], where L* (t2) is the luminosity of the protostar and t2 the time spent since the start of the warm-up phase. The function that describes how L* (t2) evolves in time was derived by Hosokawa & Omukai (2009) and depends on the mass accretion rate, Ṁ, assumed for the protostar. Based on the very few observational estimates (e.g. Beuther et al. 2002b; Herpin et al. 2012; Wang et al. 2013), we explored two typical cases: Ṁ = 10−5 M⊙ yr−1 and Ṁ = 10−3 M⊙ yr−1 (see Figs. 4and 12 in Hosokawa & Omukai 2009). These values also agree with the results of Peters et al. (2011), who simulated the collapse of a magnetised rotating molecular cloud of 1000 M⊙. During phase II, we assumed that the gas radial density profile of the clump and its dust content remained constant over time so that the mass accretion rate of the protostar was completely balanced by the average infall rate observed in a massive clump (e.g. Schneider et al. 2010; Peretto et al. 2013; Wyrowski et al. 2016; Liu et al. 2018).

Analogously to phase I, we assumed that the gas and dust temperatures were in equilibrium; this assumption is less accurate at small radii (i.e. ≲50 AU), where the temperatures are close to the dust sublimation limit, and at large radii (i.e. ~105 AU), where the gas-dust collision term is subdominant with respect to the radiation coupling (Draine 2011). However, neither region is relevant in our analysis and does not affect our findings. An example of how the Td radial profile evolves is shown in Fig. 2 for three (typical) times, Ṁ = 10−3 M⊙ yr−1, and for the following parameters Rc = 104.5 AU and nH(r0) = 107 cm−3. For a qualitative comparison of the Td radial profiles shown in Fig. 2, we refer to the recent results of Gieser et al. (2021), who derived the temperature structure of 18 massive star-forming regions based on two of the four tracers considered in this work (i.e. H2CO and CH3CN; Sect. 2). In this case, the authors assumed a power law (i.e. T(r) ∝ r−q; see their Eq. (1)) to fit the radial temperature in each source, providing an average slope of |q| = 0.4 ± 0.1. Our profiles in Fig. 2 show a similar slope (|q|~ 0.5) at distances larger than the dust sublimation radius (i.e. the flat region in the same figure) and comparable temperatures at their fiducial radius (see Table 3 in Gieser et al. 2021). A more detailed description of the effect of L* (t2) and nH (r) on Td (r, t2) will be discussed in a forthcoming paper (Grassi et al., in prep.).

To compare our results with the observed properties of the TOP100 clumps (Giannetti et al. 2014; König et al. 2017), we evolved one model for each mass in Table 2 (derived from nH (r0) and Rc in Table 1), and for each value of b and Ṁ in Table 1. We considered only clumps with a total mass higher than 15 M⊙, which is the lowest mass associated with a TOP100 source, for a total of 54 models.

thumbnail Fig. 2

Dust temperature radial profile at different times for an arbitrary model with the gas number density profile described by Eq. (2) and Rc = 104.5 AU and nH(r0) = 107 cm−3. With time, Td increases, driven by the protostar luminosity. The flat inner part represents the sublimation of the dust grains. The minimum temperature is set to 15 K to ensure continuity with the collapse phase.
Table 2

Summary of the total masses of each model.

3.2 Chemical model

The chemical evolution of a massive clump (i.e. the collpase and the warm-up phases) is described by employing the publicly available time-dependent code KROME2 (Grassi et al. 2014).

Adsorption and desorption3 processes were included as in Hasegawa et al. (1992) and Hasegawa & Herbst (1993). The surface reactions between two species i and j follow Semenov et al. (2010), where the rate coefficient, in units of cm3 s−1, is (3)

with the diffusion through thermal hopping (; Semenov et al. 2010), Td is the dust temperature, is the characteristic Debye vibration frequency for the adsorbed species, nS = 1.5 × 1015cm−2 is the surface density of binding sites, mi is the mass of the species, and is the binding energy of the ith species on the binding site. In Eq. (3), is the dust number density, where is the dust mass, and the total number of binding sites of a grain assumes an average distance between two contiguous sites of app = 3 Å (Hocuk & Cazaux 2015).

The probability for a reaction to occur is Pij = αij exp(−Ea∕kBTd), where Ea is the activation energy of the reaction and αij is a parameter that depends on the number and on the type of species in the reaction. In the case of exothermic reactions (i.e. Ea = 0), (I) if i≠ j, Pij = 1; (II) if i = j, Pij = αij = 1∕2. For endothermic reactions (i.e. Ea≠0), αij is the inverse of the number of paths in the branching ratios.

The final chemical network was derived from Semenov et al. (2010)4 and contains 654 chemical species and 5869 gas-phase, gas-grain, and grain-surface reactions. The photochemical rate coefficients follow the Av formalism as in Draine (1978) (see KIDA5). In order to ensure that each depleted species was released into the gas phase, we modified the original network by adding 70 missing desorption processes, with the binding energy values updated to the most recent estimates in KIDA (see Table A.1 and Table A.2 for the reactions, and Appendix B for the network benchmark).

The abundances of chemical species evolve with time, starting from the assumed initial conditions following the recent large-scale simulations of molecular clouds (e.g. Hocuk et al. 2016; Clark et al. 2019), showing that CO is already formed at densities of few × 103 cm−3. Following these findings, we assumed H, C, and O to be in molecular form. In particular, the abundances of H2, H, He, N, O, CO, and N2 were taken from Bovino et al. (2019), while for the other elements, we refer to Garrod & Herbst (2006), as reported in Table 3.

Table 3

Summary of fiducial initial elemental abundances, ni, with respect to the abundance of H-nuclei, nH.

4 Results

An example of the chemical evolution of some observed species during phase I is reported in Appendix C, while here we focus on phase II and the comparison with observations. However, because the outcome of phase I represents the initial chemical conditions for phase II, a few considerations for phase I should be made. The fractional abundances are considerably affected by the dynamics of the collapse, showing a difference up to ~3 orders of magnitudes between the fastest (b = 1) and the slowest (b = 0.1) collapse, even if the global trend remains generally unaltered (more details in Appendix C). Based on the continuity between phase I and phase II, the final effect that b has on the models is to provide different inputs for the warm-up phase. For this reason, we ran phase II for each of the 54 models defined by the parameter space, taking the effect of parameter b on the second phase of the model into account (see also Sect. 3.1.2).

Figure 3 shows the abundances of H2CO, CH3CCH, CH3CN, and CH3OH as a function of time and distance in the gas phase (left panels) and on grains (right panels; labelled ‘(g)’) as a result of phase II. The white contours show the temperature, and the red curves the evaporation fronts, namely the distance from the protostar at which the desorption timescale of a given species (defined as the inverse of the desorption rate, involving the thermal and the cosmic-ray induced desorption; Semenov et al. 2010) is shorter than the dynamical timescale of the model. As a consequence, at distances larger than the evaporation front, the abundances of the tracers on dust grains (Fig. 3) increase rapidly. At the same time, the temperatures at these distances are low enough to enhance two-body reactions on the surface of dust grains. Hydrogenation chains play a major role in the formation of CH3CCH, CH3CN, and H2CO on dust grains, and methanol is also formed through the reaction (g)CH3 + (g)OH →(g)CH3OH.

Each evaporation front grows as a function of time following the temperature evolution and ranges between ~103 and ~105 AU. Our findings agree with those of Choudhury et al. (2015), where the efficiency of the evaporation as a function of the thermal evolution of the clump is responsible for evaporation scales similar to ours (see their Fig. B1).

thumbnail Fig. 3

Temporal and radial evolution of the abundances of the tracers observed in the TOP100 (Giannetti et al. 2017; Tang et al. 2018) during phase II in the gas phase (left panels) and on dust (right panels) for the same reference model as in Fig. 2. White contours indicate the temperature computed with MOCASSIN. Red curves corresponds to the evaporation front of the given tracer.
thumbnail Fig. 4

Panel a: exampleof the LOS gas-phase column density profiles obtained at the end of the post-processing described in Sect. 4.1 for models with b = 1 and Ṁ = 10−3M⊙ yr−1 and for an arbitrarily chosen source at a distance of 3.5 kpc. For the U_W items in the legend, U and W are Rc (colours) and nH(r0) (styles), respectively. Panel b: radial distribution of the gas-phase column densities extracted at the positions A, B, and C in panel a, i.e. the same times as in Fig. 2.

4.1 Post-processing of the model outputs

The chemical structure of each clump was reconstructed assuming that each radial profile from phase II (i.e. each column in the relative abundance maps in Fig. 3) represents a radius of a spherically symmetric clump. For each profile we generated a data cube of 150 × 150 × 150 pixels at a resolution of ~103 AU pixel−1, which corresponds to one-half of the size of the smallest evaporation front radius in Fig. 3, in order to always sample the region at which the gas-phase abundances are enhanced by the desorption of the products of dust-phase chemistry. The column density maps of each tracer were then obtained by integrating the abundances along the LOS, applying for each tracer a convolution to the APEX resolution of the observed transitions6 (see Appendix D), and taking the different distances of the TOP100 objects into account. The final column density profile of a given tracer X, Nmod(X), was calculated along the LOS at the centre of the clump.

With respect to this procedure, CH3OH and CH3CN required an additional constraint to reproduce the observed higher-K lines. In their spectral fit, Giannetti et al. (2017) distinguished the hot and cold components, assuming two domains separated at 100 K and providing the individual column densities. Following the same approach provides additional constraints to reduce the uncertainties and to test the reliability of the temperature of the background model, especially in the innermost part of the clump. After reconstructing the spherical distribution of CH3CN and CH3OH, we applied a 100 K temperature-mask, taking only regions into account in which T >100 K, to generate the averaged column densities profiles of the hot components alone.

As an example, we report the final gas-phase CH3CCH profiles in Fig. 4a, obtained assuming a distance of 3.5 kpc (arbitrarily chosen) to apply the post-processing. We show the results for Ṁ = 10−3 M⊙ yr−1 and for the standard free-fall collapse model (i.e. b = 1), and we note that nH(r0) is the parameter that mostly affects the evolution of Nmod(CH3CCH).

This behaviour is the consequence of two combined effects. First, the initial abundances of the warm-up phase are scaled by the density profile of each model, and thus, models with higher values of nH (r0) have higher initial column densities (see Fig. B.2). Second, the different thermal structures of the clumps have a strong effect on thechemical evolution. At a given time, the average temperature along the LOS is lower in clumps with larger nH (r0) because the radiation of the protostar is more attenuated. Because the luminosity of the protostar increases with time and affects the temperature of the surrounding gas, the evaporation front grows and so does the number of chemical species released into the gas. In Fig. 4b we show the radial profile of gas-phase N(CH3CCH) at three times (A, B, and C in panel a). The column density of this species increases by about one order of magnitude over time, caused by the expanding evaporation front and by the decrease in CH3CCH adsorption rate in the outer parts of the model. The adsorption process is directly proportional to the density, which drops at r > Rc (Fig. 3). This also provides an explanation for the observed increasing trend in the gas-phase CH3CCH abundance with the evolutionary stage of massive clumps (e.g. Molinari et al. 2016; Giannetti et al. 2017).

4.2 Comparison of the modelled column densities with sources from the dTOP100 sample

To infer the ages of the clumps and the relative duration of the evolutionary phases in the TOP100 sample, we compare the column densities derived from the observations (Appendix D and Table D.1) with those from our post-processed models. For a given clump the dataset, D = {L, M, N(X)}, comprises eight (two+six) measurements, that is, the luminosity L, the mass M, and the column densities of the six observed chemical tracers, N(X). At each time step of the warm-up phase, with the free parameters θ = {b, Ṁ, Rc, nH(r0)} given the data D, the likelihood is (4)

where PC(θ|D), PM (θ|D), and PL (θ|D) are probability density distributions of the chemical abundances, the mass, and the bolometric luminosity of each clump, respectively.

The first factor in Eq. (4) is (5)

where for the sake of simplicity, we omitted the arguments of PC, and where (Garrod et al. 2007; Gerner et al. 2014) (6)

with ‘erfc’ the complementary error function. In Eq. (6), Nobs(X) and Nmod(X) are the observed and modelled column densities, respectively. As in Garrod et al. (2007), we set σ = 1, which corresponds to a difference of one order of magnitude between the observed and the modelled column densities. This assumption takes into account the uncertainties on the observed column densities, on the initial chemical conditions assumed for the phase I, and on the shape of the density profile assumed in phase II (see Sect. 3.1.2). For sources without a detection, we assumed a detection limit (see Appendix D) and set PC,X = 1 when Nmod(X) < Nobs(X), while we used Eq. (6) when Nmod(X) > Nobs(X).

For eachclump, PL and PM are Gaussian functions with mean equal to the L and M values derived by König et al. (2017). The standard deviations were computed as follows: We added in quadrature a 50% uncertainty associated with the models based on the approximation of a single clump size of 0.5 pc and a constant accretion rate over time (Sect. 3.1.2), in addition to the errors of 50% on L and 20% on M proposed by Urquhart et al. (2018). This gives a standard deviation for PL of about 60% of L, and for PM of about 50% of M. We note that an additional variation of 10% on these uncertainties does not produce significant variations in our final results.

Figure 5a shows an example of how varies as a function of time when models with the same mass accretion rate (Ṁ = 10−3 M⊙ yr−1) are considered. The models are sorted in ascending order by the difference in mass between the values reported in König et al. (2017) and those of the models (i.e. Table 2), while t2 is the time spent in the warm-up phase defined in Sect. 3.1.2. In both panels, PL defines the range of time spent by the protostar in the warm-up phase (i.e. the range of time in between the vertical cyan lines), which only depends on the uncertainties associated with the bolometric luminosity of each source. Analogously, PM is constrained by the mass range found for a given source (e.g. see discussions in König et al. 2017 and Urquhart et al. 2018), and it is therefore possible to define a lower mass limit that is indicated by the horizontal green line. The combined information given by PL and PM limits the age range of each source. In this area, the models with the same density profile defined by Eq. (2), but different values of b, are degenerate, showing the same as a function of t2 (see Fig. 5b). PC is fundamental to break this degeneracy (see Fig. 5a), providing further information on the dynamics of the collapse that led to the observed chemical properties of each clump, and placing additional constraints on determining the relative duration of the evolutionary classes in the TOP100.

The absolute time is hence the sum of the collapse timescales, t1, and the time spent during the warm-up phase, t2. The collapse time depends on the setup of phase I, and it ranges between (0.4–5) × 106 yr depending on b and nH(r0), while the warm-up time depends on the mass accretion rate of phase II (see Sect. 3.1.2).

The absolute age for the ℓth source is the weighted mean, (7)

where i ranges over the 54 models, j over the 100times used to sample phase II, and we used the log-likelihood as weight.

We performed the comparison described above only for the sources in the TOP100 sample that agreed within the uncertainties with the masses of the models in Table 2. This reduced the sample to 48 objects, that is, 6 sources in the 70w stage, 16 IRw, 15 IRb, and 11 HII regions. Nevertheless, this limitation does not affect our approach because our model employs a subset of objects that are representative of the whole population, preserving their global features (e.g. Csengeri et al. 2016; Giannetti et al. 2017; König et al. 2017). We summarise the observed properties of the clumps in this subsample in Table D.1. It is worth noting that while the TOP100 sample does not show bias in terms of the source distances for each evolutionary phase (e.g. König et al. 2017), in Table D.1 only the IRw class extends beyond ~10 kpc. To verify whether the distance has an effect on our results, we performed a test considering sources up to a maximum distance of 6 kpc. We did not find significant changes in the final durations reported in Sect. 5.

thumbnail Fig. 5

Panel a: example of the matrix calculated for G353.41-0.36 in all the models with Ṁ = 10−3M⊙ yr−1. Models on the y-axis are sorted from the top to the bottom by the difference between the mass of the models (see Table 2) and the value associated in the TOP100 (see Table D.1). For the sake of clarity, the colour bar is limited to five orders of magnitude and is normalised to maximum likelihood (). Panel b: zoom on the age limits identified by PL and PM showing howthe likelihood would appear if PC in Eq. (4) were not taken into account. Cyan vertical lines indicate the range of time that corresponds to the uncertainties on the observed luminosity in the TOP100, and the horizontal green line shows the same limit associated with the clump mass. Times and are shown in log-scale.

5 Discussion

In this section, we estimate the relative duration of the evolutionary phases identified in the TOP100, following a different approach compared to the statistical lifetimes. Our findings are then compared with the relative number of objects observed in ATLASGAL (Urquhart et al. 2018), and we discuss whether the synthetic column densities of the selected chemical tracers match the observed densities.

5.1 Estimates of the duration of evolutionary phases

To compare the ages of the evolutionary stages defined in the TOP100 with the average absolute time estimated in a sample of objects, , we defined the age factor of the ℓth source as (8)

In this context, the models provide yr starting from t1 = 0 (i.e. the beginning of phase I; see also Fig. 1). The average fA were estimated asthe median of the age-factor distributions in each evolutionary class, computed with Eq. (8) and reported as purple markers in Fig. 6. For the sake of clarity, Fig. 6 shows the corresponding value of yr that is associated with each fA. We quantified the duration of an evolutionaryphase, Δtphase, as the time between the minimum and maximum value of (i.e. the lower and upper limits of the purple shaded areas in Fig. 6). These values were computed as the 5th and 95th percentiles of the fA distributions in each evolutionary phase. The sum of the four Δtphase is the total time of the high-mass star formation process tMSF. The purple percent values in Fig. 6 (M1) indicate the contribution of each phase to the total time tMSF.

We find that following the classification of Giannetti et al. (2014) and König et al. (2017), 12% of the star formation time is spent in the early phase (70w), and the IRw stage is associated with 19% of tMSF. This suggests a fast evolution during the early stage of the massive star formation process, which in total correspond to the ~30% of tMSF. Advanced stages (i.e. IRb sources) have the longest duration (39%), while the remaining 30% of the tMSF is spent in the final stage (HII).

Assuming that the number of objects in an evolutionary stage is also representative of its duration, we compared our findings with the relative number of objects per evolutionary class in ATLASGAL. This is represented as black percentages in Fig. 6 (“O”; from Urquhart et al. 2018). The duration found for the early stages agrees with the observed classification, and the later stages are probably biased by the different definition of IRb and HII in the TOP100 and ATLASGAL. In particular, when radio continuum emission is found at either 4 or 8 GHz within 10′′ of the ATLASGAL peak, the source is classified as HII in the TOP100 (e.g. König et al. 2017). Different criteria have been applied in Urquhart et al. (2018) to classify the advanced stages that contain radio bright HII regions, massive young stellar objects, and sources associated with methanol masers. These surveys have a different coverage than the complete ATLASGAL and also have different sensitivities (see Urquhart et al. 2014 for more details), so that the final number of objects in the advanced stages might be underestimated. Because this limit affects the separation between IRb and HII, we note that if we consider a single phase to describe the more evolved sources (i.e. IRb + HII), the total duration associated with this phase would agree better.

An additional source of uncertainty is a different value of Ṁ throughout the evolutionary sequence. Evidence of an increasing mass accretion rate in the intermediate stages of the massive star formation process is discussed in Beuther et al. (2002b), but when the protostars are close to the main sequence, their radiation pressure might slow down or quench the mass accretion (e.g. Nakano et al. 1995; Stahler & Palla 2005; Klassen et al. 2012). To quantify the effects of Ṁ, we repeated the calculation of the duration of each phase by separating the models by accretion rate and mixing the results of the different phases in all their possible combinations. We find that different Ṁ produces different durations and tMSF. The most accurate solution to interpret the observations in ATLASGAL shows an accretion rate that is initially slow and increases during the two intermediate classes, to decrease again in the more evolved phases (see the orange symbols and percentages in Fig. 6; M2). The latter result allows us to verify how reliable the assumption of a constant H2 radial density profile is that we made for the clumps during phase II (see Sect. 3.1.2). In the most variable scenario where Ṁ = 10−3 M⊙ yr−1 over the t2 identified in our models (and assuming no mass replenishment), the majority of the clumps shows a mass loss of ≲20%. Only two clumps of ~20 M⊙ would suffera relevant mass loss (i.e. G353.07+0.45 and G316.64-0.09 in Table D.1). However, these two sources belong to the IRb and IRw phases, which contain the majority of the clumps, and can be considered outliers in these classes. Conversely, for sources in the 70w and HII evolutionary phases, where Ṁ = 10−5 M⊙ yr−1, the mass loss is negligible. Alongside the few observational estimates available for Ṁ and for the infall rate (e.g. from the parent filament; see Sect. 3.1.2), these results mean that the assumption of a constant-density H2 profile is reasonable on a clump scale.

The individual durations estimated in M2 lead to ~5 × 104 yr for younger objects (i.e. 70w), ~1.2 × 105 yr for the IR-weak, ~2.4 × 105 yr for IR-bright sources, and ~1.1 × 105 yr for HII regions. Additionally, the average fA of each classincreases with the evolution of the sources (circles and stars in Fig. 6), so that the typical object of each class appears statistically older than the average from the previous evolutionary phase. We therefore assume that M2 is the most representative model, with a total massive star formation time yr (1.1 ).

In this context, the chemistry plays a crucial role in breaking the model degeneracy when the likelihood is defined by PM and PL alone and in matching the results from ATLASGAL (see Fig. 6). If PC is not included in (Eq. (4)), the duration of the advanced stages shows a difference up to a factor of ~8 when compared to those predicted by M2. This corresponds to a total time, tMSF, 4.5 times shorter, and the final relative durations differ from the observed relative number of objects in each class by up to a factor of 5.

thumbnail Fig. 6

Summary of the final durations estimated in each evolutionary class defined for the TOP100. The relative number of massive clumps in each phase and observed in ATLASGAL are shown in black (“O”; Urquhart et al. 2018). Purple shaded areas indicate the 5th and 95th percentiles of the age factor, fA, and the corresponding value in yr in our sample (model “M1”), and purple markers (‘–’) represent their medians. The numbers associated with the model “M2 ” (orange), indicate the same quantities obtained by separating the models with respect to Ṁ before calculating the average age of each source (see text). The different markers in M2 show the combination of Ṁ that best match the observed relative number of object in each phase (legend).

5.2 Selected tracers and chemical clocks

Four molecules were selected in this work: formaldehyde, methyl acetylene, acetonitrile, and methanol (see Sect. 2). They all manifest an observed upward trend in their abundances with increasing evolutionary stage, which is indicated by the luminosity-to-mass ratio of the clumps (Giannetti et al. 2017; Tang et al. 2018). Two thermal components were also required to reproduce the observed emission of the higher K-ladders of methanol and acetonitrile lines in Giannetti et al. (2017).

To evaluate to which extend the models are able to reproduce the observed column densities, the same procedure as we used to calculate Δtphase was employed (Eqs. (7) and (8)). Figure 7 summarises the comparison of the modelled and observed column densities, grouped by tracer type (colours) and evolutionary class. Shaded areas represent the range of the observed beam-averaged column densities (Sect. 2 and Appendix D), with an order of magnitude uncertainty (Eq. (6)), not shown in the figure. The median values of these observed column densities are shown as a solid line. The median column densities predicted by our models are plotted as coloured circles and are shown with an uncertainty of one order of magnitude, the same as assumed for the observations.

5.2.1 Methyl acetylene and acetonitrile

The reliability of CH3CCH and CH3CN as thermometers has been widely discussed in the past (e.g. Zhang et al. 1998; Molinari et al. 2016; Giannetti et al. 2017). In particular, Giannetti et al. (2017) detected an increasing trend of the temperature and column density of CH3CCH from the less to the more evolved sources in the TOP100 (i.e. shaded area in panel a of Fig. 7). The same behaviour is found for CH3CN if the contributions of a hot and cold component are taken into account.

In Fig. 7 the observed column densities of methyl acetylene (panel a) and acetonitrile (panels b, and c) derived from the observations increase by a factor of 10 between the 70w and the HII stages. The column densities provided by our models reflect the same trends and reproduce the observed values. The discrepancy between models and observations is within the uncertainties.

The robustness of CH3CCH and CH3CN as chemical clocks was finally tested by removing these two tracers from PC. This had the same effect as removing PC from altogether (see Sect. 5.1), and we conclude that methyl acetylene and acetonitrile are effective chemical clocks for characterising the evolution of massive star-forming clumps and that they contribute to constraining our findings.

5.2.2 Formaldehyde and methanol

Unlike CH3CCH and CH3CN, the model column densities of methanol and formaldehyde (panels d, e, and f in Fig. 7) agree less well with the observed densities. The slightly increasing trend in the observed column densities with evolution is reproduced by the models, but the modelled column densities under- or over-estimate those observed for the TOP100 by one order of magnitude at least, similarly to Gerner et al. (2014). This might be related to the uncertainties in the chemical pathways that determine the formation ofmethanol and its precursors. For temperatures below 20 K, hydrogenation chains are usually invoked to convert CO into H2CO and CH3OH on grains because hydrogenation is made very efficient by the strong CO depletion that occurs at low temperatures (e.g. Caselli et al. 2008, Giannetti et al. 2014 and Sabatini et al. 2019). However, at the same temperatures, thermal desorption is inhibited and methanol is not efficiently released from the dust into the gas phase, suggesting that the surface chemistry alone is not capable of explaining the detection of gas-phase methanol in the cold phases of massive star-forming regions as well (e.g. Cosentino et al. 2018). Alternative formation paths and mechanisms have been proposed to solve this issue (e.g. Viti & Williams 1999; Garrod et al. 2007; Vasyunin & Herbst 2013), but this was found to be relatively inefficient in reproducing the observed abundances of methanol (Geppert et al. 2006), and to no longer accurately predict the abundances of other chemical species (Garrod et al. 2007). The same issue also concerns H2CO, a molecular precursor of methanol. In the advanced evolutionary stages, the amount of methanol is underestimated (although closer to the observed values). The reason might be that the thermal evaporation of methanol in the warmer environments is not well determined, or that further formation pathways are missing in the colder phases. We also note that the absence of the quantum tunneling diffusion in our models might affect the final abundance of methanol, enhancing the efficiency of CO hydrogenation in the colder evolutionary phases, and thus influencing the final abundance of methanol in the warm phases (Vasyunin et al. 2017).

An additional uncertainty that we explored is the effect of different activation energies for reactions that drive the formation of molecules such as formaldehyde and methanol at low temperatures. Recently, a study of CO hydrogenation on water ice by Rimola et al. (2014) reviewed the activation energy of the hydrogenation reaction between H2CO and H producing CH3O, the precursor of methanol. The authors proposed a barrier of 1300 K, which would favour the formation of CH3OH compared to the value commonly employed (~2500 K, as in this work and derived for the gas phase; Woon 2002). Even a lower activation energy of ~500 K has been predicted by Fuchs et al. (2009). However, we have tested these two scenarios without finding significant changes with respect to the modelled column densities of the tracers in Fig. 7. The most substantial change concerns CH3OH when the extreme activation energy proposed by Fuchs et al. (2009) is assumed. In this case, methanol column densities approach the values observed in the two most advanced evolutionary stages, while those of formaldehyde remain always largely overestimated. This suggests that the activation energy has secondary effects in reproducing the observed abundances of methanol and formaldehyde, and other processes or observational strategies, such as different or multiple transitions, should be taken into account.

We finally note that the discrepancies in the abundances of H2CO and CH3OH are not relevant for our results. When these two tracers are removed from PC and , the duration of each phase and the final tMSF remain unaffected, showing a variation of ~10% with respectto the numbers reported in Fig. 6.

thumbnail Fig. 7

Comparison of the median column densities observed in the TOP100 (solid lines) with those predicted by our models (circles) that were obtained by applying the same procedure as was used to quantify Δtphase (see Sect. 4.2). Panels are separated by tracers (colours). The shape of the shaded areas indicates the minimum and maximum observed column density, and the error bars associated with each circle incorporate the uncertainty of one order of magnitude assumed in Eq. (6) for the comparison. The modelled column densities of the hot component of methanol (circles in panel f) are multiplied by a factor 1000. Triangles indicate the upper limits. The column densities are shown in log-scale.

6 Summary and conclusions

We have presented a new and generalised method for deriving the evolutionary timescales of the massive star formation process. Compared to typical statistical approaches (e.g. Russeil et al. 2010; Mottram et al. 2011; Tigé et al. 2017; Motte et al. 2018; Urquhart et al. 2018), our method follows a different procedure in which the tMSF is the final result that is derived as the sum of the individual Δtphase (Sect. 5.1 and Fig. 6), and not an a priori assumption.

We developed a set of 54 models, built to represent the entire population of massive clumps of the ATLASGAL-TOP100 sample (Sect. 1). The models consist of two physical phases: an initial isothermal collapse followed by a warm-up phase induced by a massive protostar at the centre of a spherical clump (Sect. 3). In particular, we assumed a Plummer-like density profile with a slope of 5/2, while the temperature profiles were derived from accurate 3D RT simulations. The recent empirical estimates of the temperature profiles obtained by Gieser et al. (2021) from the same tracers used in this work agree with our computed profiles (see Sect. 3.1.2). We then compared the column densities of formaldehyde (H2CO), methyl acetylene (CH3CCH), acetonitrile (CH3CN), and methanol (CH3OH) as observed in the ATLASGAL-TOP100 sources (Giannetti et al. 2017; Tang et al. 2018) with the modelled densities by post-processing the outputs of the models at the same angular resolution as the observed data (see Appendix D).

The timescales of the evolutionary stages associated with the massive star formation process were derived by considering the physical properties of the clumps, that is, their mass and luminosity, and the observed abundances of each selected molecular tracer. Considering a mass accretion rate depending on the evolutionary phases, we find a total star formation time, tMSF ~ 5.2 × 105 yr, which agrees well with that assumed by the statistical methods. This provides a new and significant validation of the statistical methods. The individual Δtphase that define tMSF are found to be ~5 × 104 yr for the 70w, ~1.2 × 105 yr for the IRw, ~2.4 × 105 yr for IRb sources, and ~1.1 × 105 yr for HII evolutionary stages of the TOP100 sample (Sect. 2). Knowing tMSF and Δtphase, we derived the relative duration of each phase and found an agreement with the relative number of objects in each phase observed in the ATLASGAL survey. This provides further confirmation of the reliability of the method we presented. The chemical constraint included in the likelihood to determine the duration of the different phases is necessary to achieve results that agree with the observed relative number of objects in the ATLASGAL survey. Without this additional constraint, Δtphase becomes shorter by up to a factor of ~8 with respect to those reported above, and the relative durations differ by up to a factor of 5 with the relative number of objects in each class.

Of the selected molecular tracers, CH3CCH and CH3CN are best reproduced by our models, and H2CO and CH3OH differ from the observed values, although they do follow the observed trends. Therefore the final Δtphase reported in this paper are mainly constrained by CH3CCH and CH3CN. Because the chemical constraint is important to identify reliable timescales, we plan to extend the number and complexity of the selected molecular tracers with several chemical clocks proposed in the literature (e.g. HC3N by Taniguchi et al. 2018; CH3OCHO and CH3OCH3 by Coletta et al. 2020; see also Urquhart et al. 2019; Belloche et al. 2020). This might help to better assess the different phases of the massive star formation process and increase the number of constraints on Δtphase reported here.

We also found that the evaporation fronts of the discussed molecular tracers vary between 103 and 105 AU during the warm-up phase (Sect. 4), which in our models marks the regions in which each tracer becomes abundant inthe gas phase. This result suggests that we compare our findings with the observations provided by the new astronomical facilities. In particular, the Atacama Large Millimeter/submillimeter Array (ALMA; Wootten & Thompson 2009) offers the perfect opportunity of achieving the high angular resolutions (e.g. Csengeri et al. 2018; Maud et al. 2019; Sanna et al. 2019; Johnston et al. 2020) needed to sample scales of about the physical sizes of the modelled evaporation fronts. Moreover, ALMA can also reach the sensitivities needed to detect a large number of components in the band. This would allow the removal of possible opacity effects and the accurate definition of the thermal state of the clump that can be compared with the results of this work. The increasing complexity of chemical models and the progressmade in the observational techniques make this paramount goal increasingly achievable in the near future.

To conclude, we have reported relevant and robust results on the high-mass star formation process together with reliable estimates of the duration of the different phases of this complex process. The present pipeline is based on 1D models, which cannot capture the effect of dynamical processes such as magnetic fields and turbulence on the density evolution of the collapsing clumps and the subsequent (proto-)stellar accretion process fully. On the other hand, low-dimensionality models allow distinguishing between the different chemical processes, and building a large number of models to infer statistical properties that can be compared with observations. Although some authors have recently included the evolution of large chemical networks in 3D models (see e.g. Bovino et al. 2019), in the presence of a (proto-)stellar object, modelling the chemistry and the thermal evolution of the gas in 3D still represents a challenge.

Acknowledgements

The authors wish to thank the anonymous referee, for her/his suggestions to improve the manuscript, D. Semenov for fruitful discussions and feedbacks about the chemical model benchmark, and J. Ramsey how developed the KROME python-interface employed in this work (see also Gressel et al. 2020). This paper is based on data acquired with the Atacama Pathfinder EXperiment (APEX). APEX is a collaboration between the Max Planck Institute for Radioastronomy, the European Southern Observatory, and the Onsala Space Observatory. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018; see also http://www.astropy.org), of NASA’s Astrophysics Data System Bibliographic Services (ADS) and of Matplotlib (Hunter 2007). G.S. acknowledges R. Pascale for useful suggestions and comments, and the University of Bologna which partially provided the funds (Marco Polo fellowship) for this project. S.B. acknowledges for funds through BASAL Centro de Astrofisica y Tecnologias Afines (CATA) AFB-17 002.

Appendix A Details of the chemical network

Table A.1 shows the types of chemical reactions that we included in the network. Column (1) contains names of each type of reaction; Col. (2) an example of chemical reaction between two generic reactants A and B, and, Col. (3) summarises the total number of reactions included in the network.

In Table A.2 we provide the complete list of binding energies assumed in this work, updated to the most recent estimates found in KIDA. Where it was not possible to find a recent estimate, we used Semenov et al. (2010) (species in bold). The 35 chemical species in the bottom part of the table are those addedto the original network due to the missing desorption processes (see Sect. 3.2).

Table A.1

Summary of the reactions in the chemical network.

Table A.2

Complete list of binding energies.

Appendix B Benchmark of the chemical network

We benchmarked our network against Semenov et al. (2010). The initial conditions are summarised in Table B.1. All the elements are initially atomic, with the exception of hydrogen, which was assumed to be completely in molecular form. Except for He, N, and O, all the elements are also ionised, and grains are initially neutral.

TMC1-model: we benchmarked the physical case called TMC1 model. The gas has a constant temperature of 10 K, a constant hydrogen nuclei density nH = 2 × 104 cm−3, and visual extinction Av = 10 mag. The sticking coefficient was S = 1, and the cosmic-ray ionisation rate of hydrogen molecules was set to ζ2 = 1.3 × 10−17 s−1 (Spitzer & Tomasko 1968; Glassgold & Langer 1974). We have assumed an average grain size of m and a bulk density ρ0 = 3 g cm−3, corresponding to silicate grains, as in Sect. 3.1. Figure B.1 shows the perfect agreement with the results reported by Semenov et al. (2010).

Table B.1

Summary of fiducial initial elemental abundances set for the TMC1 model.

thumbnail Fig. B.1

Comparison of the time-dependent variation in the abundances of four arbitrary species (i.e. C, H2 O, CO, and CH3OH) in gas phase(left panels) and on dust (right panels) in the TMC1 model. The green line shows the results obtained with the KROME -package, and the blue crosses are the results of Semenov et al. (2010) with ALCHEMIC.
thumbnail Fig. B.2

Fractional abundance evolution for H2CO, CH3 CN, CH3 CCH, and CH3 OH as a functionof the central density of the clump during the collapse phase. Colours indicate the value of b assumed for the collapse described by Eq. (1). Solid lines show the results employing the canonical binding energies from KIDA, and shaded areas represent the solutions when Tb was varied by ±10%.

Appendix C Chemical evolution and model uncertainties during phase I

Figure B.2 shows the fractional abundance evolution of the relevant tracers (see Sect. 1) as a function of the central density of the clump during the isothermal collapse. We first tested the effect of the collapse velocity (parameter b in Eq. (1)) on the chemical evolution of some of the observed species by delaying the collapse by 50 and the 90%, that is, b = 0.5 and b = 0.1, respectively.Different colours in Fig. B.2 represent different values of b.

Despite the effects induced by varying b, the relative importance of the chemical reactions involving some of the main tracers remains unaltered. In particular, H2CO is mainly formed through CH3 + O →H2CO + H, and it is destroyed by ion-molecule reactions that involve H+, C+, and S+. For densities higher than ~105 cm−3, ion-molecule reactions become predominant and the abundance of H2CO decreases. At higher densities, some relative abundance profiles again increase because the dominant reaction becomes H3CO+ + e-→H2CO + H.

Methyl acetylene is formed by the dissociative recombination (i) , the main gas-phase source of CH3CCH (see Hickson et al. 2016). For densities above ~106 cm−3, the relative abundance is almost constant because CH3CCH reacts with C with a rate coefficient similar to (i).

Conversely, the abundance of acetonitrile is more sensitive to b. The main source of gas-phase CH3CN in our model is (ii) , followed by the dissociative recombination (iii) CH3CNH+ + e−→CH3CN + H. We found that at beginning of the collapse, CH3CN is rapidly formed, in agreement with what was estimated by Loison et al. (2014), and is then slowed down by (iv) . The combined effect of reactions (iii) and (iv) is also sensitive to b because for fast collapse (b = 1), reactions (ii) and (iii) have less time to form CH3CN before reaction (iv) becomes effective. For densities between 105 and 106 cm−3, CH3CN starts to react efficiently with . The effect of this reaction in combination with (iv) produces a slow decrease in the CH3CN profiles at high densities.

The abundance of methanol increases with density, driven by reaction , followed by (see also Viti & Williams 1999).

To determine other sources of uncertainties, we analysed the variation in chemical evolution by changing the binding energies reported in Table A.2 by ± 10%, repeating the collapse with the same physical conditions. The effects on the fractional abundances are represented by the coloured areas in Fig. B.2. The discrepancy is more than one order of magnitude at the highest densities. We note that even larger discrepancies are reported by Penteado et al. (2017), when the binding energies are randomly selected from a normal distribution (see results reported in their Fig. 2).

These tests show that the uncertainties due to binding energies variation are not negligible (see Fig. B.2), but they are not greater than variations caused by variations in b. This suggests that the dynamical state of the clump may play an important role in regulating its chemical evolution (see also the recent results by Kulterer et al. 2020).

Finally, we tested our collapse-phase with a set-up similar to that of Garrod & Herbst (2006), but employing the chemical network andbinding energies of Semenov et al. (2010) because both are publicly available and consequently easy to benchmark. The difference between the abundances found by Garrod & Herbst (2006) and ours is within an order of magnitude. It is caused by the adopted details of the chemical networks.

Appendix D Notes onthe observed column densities

In Sect. 4.1 we discussed that we post-processed the outputs of the models to compare the synthetic column densities with a number of available observations of the TOP100 clumps (right panel of Fig. 1). In this appendix we report some notes about our treatment of the observed column densities before the comparison that we discussed in Sect. 4.2. The observed column densities we used in the comparison are summarised in Table D.1.

D.1 Formaldehyde H2CO

Formaldehyde data were taken from Tang et al. (2018). For simplicity, we assumed a single half power beam width (HPBW) size of 29′′ to convolve the modelled column densities (i.e. the APEX-telescope resolution corresponding to the J = 3−2 transition of the H2CO ortho and the para forms at ~211 and ~218 GHz, respectively). For the sources in which both isomers are observed, the total column density, Nobs(H2CO), is derived as N(o-H2CO) + N(p-H2CO). For sources not detected in o- and/or p-H2CO, we estimated the corresponding detection limit column density using WEEDS (Maret et al. 2011), with a 3σ detection with respect to the average noise level (Tang et al. 2018; 0.06 K main-beam temperature scale) and assuming the H2CO molecular line parameters provided by the CDMS8 database (Müller et al. 2001). For each evolutionary class, we assumed the average values of Tkin (i.e. 52, 73, 81, and 110 K for 70w, IRw, IRb, and HII, respectively) and a full width at half maximum (FWHM) of 5 km s−1 for the line reported by Tang et al. (2018), and that the source size was equal to the APEX beam at the corresponding frequency.

D.2 Methyl acetylene (CH3CCH), methanol (CH3OH), and acetonitrile (CH3CN)

Column densities of methyl acetylene, methanol, and acetonitrile were taken from Giannetti et al. (2017), who observed the J = 20K − 19K, J = 7K − 6K, and J = 19K − 18K bands for CH3CCH, CH3OH, and CH3CN, respectively,with APEX. Because these lines are very close to each other, we assumed the same angular resolution for APEX, that is, 18′′. The final column densities in Giannetti et al. (2017) were corrected for the APEX-beam dilution, , where θS and θbeam are the source and beam HPBW sizes, respectively. We (re-)applied this correction to each source where the detections were not resolved in the APEX-beam angular size (i.e. when θS < θbeam) to obtain the beam-averaged column densities, whereas for sources with θS > θbeam, η = 1. Detection limits were estimated assuming the mean noise of the source spectra, the average line width, and the median excitation temperatures reported in Tables 6 and 7 of Giannetti et al. (2017) for each evolutionary class. We treated the observed hot components for these tracers (see Sect. 4.1) in the same way as the colder ones, considering the appropriate temperatures reported by Giannetti et al. (2017), when we calculated the detection limits.

Table D.1

Summary of the observed properties of the selected TOP100 sources.

References

  1. Adams, F. C. 2010, ARA&A, 48, 47 [Google Scholar]
  2. André, P., Revéret, V., Könyves, V., et al. 2016, A&A, 592, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Arzoumanian, D., André, P., Didelon, P., et al. 2011, A&A, 529, L6 [CrossRef] [EDP Sciences] [Google Scholar]
  4. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  6. Awad, Z., Viti, S., Collings, M. P., & Williams, D. A. 2010, MNRAS, 407, 2511 [Google Scholar]
  7. Bally, J., Moeckel, N., & Throop, H. 2005, ASP Conf. Ser., 341, 81 [Google Scholar]
  8. Battersby, C., Bally, J., & Svoboda, B. 2017, ApJ, 835, 263 [Google Scholar]
  9. Belloche, A., Maury, A. J., Maret, S., et al. 2020, A&A, 635, A198 [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bergin, E. A., & Tafalla, M. 2007, ARA&A, 45, 339 [NASA ADS] [CrossRef] [Google Scholar]
  11. Beuther, H., Schilke, P., Menten, K. M., et al. 2002a, ApJ, 566, 945 [Google Scholar]
  12. Beuther, H., Schilke, P., Sridharan, T. K., et al. 2002b, A&A, 383, 892 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Beuther, H., Zhang, Q., Bergin, E. A., & Sridharan, T. K. 2009, AJ, 137, 406 [Google Scholar]
  14. Bonfand, M., Belloche, A., Garrod, R. T., et al. 2019, A&A, 628, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Bovino, S., Ferrada-Chamorro, S., Lupi, A., et al. 2019, ApJ, 887, 224 [CrossRef] [Google Scholar]
  16. Brown, P. D., Charnley, S. B., & Millar, T. J. 1988, MNRAS, 231, 409 [Google Scholar]
  17. Caselli, P., Vastel, C., Ceccarelli, C., et al. 2008, A&A, 492, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Cesaroni, R. 2005, IAU Symp., 227, 59 [Google Scholar]
  19. Choudhury, R., Schilke, P., Stéphan, G., et al. 2015, A&A, 575, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Clark, P. C., Glover, S. C. O., Ragan, S. E., & Duarte-Cabral, A. 2019, MNRAS, 486, 4622 [Google Scholar]
  21. Coletta, A., Fontani, F., Rivilla, V. M., et al. 2020, A&A, 641, A54 [CrossRef] [EDP Sciences] [Google Scholar]
  22. Contreras, Y., Schuller, F., Urquhart, J. S., et al. 2013, A&A, 549, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Cosentino, G., Jiménez-Serra, I., Henshaw, J. D., et al. 2018, MNRAS, 474, 3760 [NASA ADS] [Google Scholar]
  24. Csengeri, T., Urquhart, J. S., Schuller, F., et al. 2014, A&A, 565, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Csengeri, T., Leurini, S., Wyrowski, F., et al. 2016, A&A, 586, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Csengeri, T., Bontemps, S., Wyrowski, F., et al. 2018, A&A, 617, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Davies, B., Hoare, M. G., Lumsden, S. L., et al. 2011, MNRAS, 416, 972 [NASA ADS] [CrossRef] [Google Scholar]
  28. Draine, B. T. 1978, ApJS, 36, 595 [NASA ADS] [CrossRef] [Google Scholar]
  29. Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium, ed. Bruce T. Draine (Princeton: Princeton University Press) [Google Scholar]
  30. Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 [Google Scholar]
  31. Elmegreen, B. G. 1998, ASP Conf. Ser., 148, 150 [Google Scholar]
  32. Ercolano, B., Barlow, M. J., Storey, P. J., & Liu, X. W. 2003, MNRAS, 340, 1136 [Google Scholar]
  33. Ercolano, B., Barlow, M. J., & Storey, P. J. 2005, MNRAS, 362, 1038 [Google Scholar]
  34. Fontani, F., Pascucci, I., Caselli, P., et al. 2007, A&A, 470, 639 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Fuchs, G. W., Cuppen, H. M., Ioppolo, S., et al. 2009, A&A, 505, 629 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Garrod, R. T., & Herbst, E. 2006, A&A, 457, 927 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Garrod, R. T., & Widicus Weaver, S. L. 2013, Chem. Rev., 113, 8939 [Google Scholar]
  38. Garrod, R. T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Garrod, R. T., Widicus Weaver, S. L., & Herbst, E. 2008, ApJ, 682, 283 [NASA ADS] [CrossRef] [Google Scholar]
  40. Geppert, W. D., Hamberg, M., Thomas, R. D., et al. 2006, Faraday Discuss., 133, 177 [Google Scholar]
  41. Gerner, T., Beuther, H., Semenov, D., et al. 2014, A&A, 563, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Giannetti, A., Wyrowski, F., Brand, J., et al. 2014, A&A, 570, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Giannetti, A., Leurini, S., Wyrowski, F., et al. 2017, A&A, 603, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Giannetti, A., Bovino, S., Caselli, P., et al. 2019, A&A, 621, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Gieser, C., Semenov, D., Beuther, H., et al. 2019, A&A, 631, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Gieser, C., Beuther, H., Semenov, D., et al. 2021, A&A, 648, A66 [EDP Sciences] [Google Scholar]
  47. Glassgold, A. E., & Langer, W. D. 1974, ApJ, 193, 73 [NASA ADS] [CrossRef] [Google Scholar]
  48. Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  49. Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126 [CrossRef] [Google Scholar]
  50. Güsten, R., Nyman, L. Å., Schilke, P., et al. 2006, A&A, 454, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Hartmann, L., Ballesteros-Paredes, J., & Heitsch, F. 2012, MNRAS, 420, 1457 [Google Scholar]
  52. Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 261, 83 [NASA ADS] [CrossRef] [Google Scholar]
  53. Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167 [NASA ADS] [CrossRef] [Google Scholar]
  54. Herpin, F., Chavarría, L., van der Tak, F., et al. 2012, A&A, 542, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Hickson, K. M., Wakelam, V., & Loison, J.-C. 2016, Mol. Astrophys., 3, 1 [NASA ADS] [CrossRef] [Google Scholar]
  56. Hocuk, S., & Cazaux, S. 2015, A&A, 576, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Hocuk, S., Cazaux, S., Spaans, M., & Caselli, P. 2016, MNRAS, 456, 2586 [Google Scholar]
  58. Hoq, S., Jackson, J. M., Foster, J. B., et al. 2013, ApJ, 777, 157 [Google Scholar]
  59. Hosokawa, T., & Omukai, K. 2009, ApJ, 691, 823 [NASA ADS] [CrossRef] [Google Scholar]
  60. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  61. Johnston, K. G., Hoare, M. G., Beuther, H., et al. 2020, A&A, 634, L11 [CrossRef] [EDP Sciences] [Google Scholar]
  62. Kennicutt, R. C. 2005, in Massive Star Birth: A Crossroads of Astrophysics, eds. R. Cesaroni, M. Felli, E. Churchwell, & M. Walmsley (Cambridge: Cambridge University Press), 227, 3 [Google Scholar]
  63. Klassen, M., Peters, T., & Pudritz, R. E. 2012, ApJ, 758, 137 [Google Scholar]
  64. Kobayashi, C., Karakas, A. I., & Lugaro, M. 2020, ApJ, 900, 179 [CrossRef] [Google Scholar]
  65. König, C., Urquhart, J. S., Csengeri, T., et al. 2017, A&A, 599, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Kulterer, B. M., Drozdovskaya, M. N., Coutens, A., Manigand, S., & Stéphan, G. 2020, MNRAS, 498, 276 [Google Scholar]
  67. Kurtz, S., Cesaroni, R., Churchwell, E., Hofner, P., & Walmsley, C. M. 2000, in Protostars and Planets IV, eds. V. Mannings, A. P. Boss, & S. S. Russell (Tucson: University of Arizona Press), 299 [Google Scholar]
  68. Liu, X.-L., Xu, J.-L., Ning, C.-C., Zhang, C.-P., & Liu, X.-T. 2018, Res. Astron. Astrophys., 18, 004 [Google Scholar]
  69. Loison, J.-C., Wakelam, V., & Hickson, K. M. 2014, MNRAS, 443, 398 [Google Scholar]
  70. Maret, S., Hily-Blant, P., Pety, J., Bardeau, S., & Reynier, E. 2011, A&A, 526, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Maud, L. T., Cesaroni, R., Kumar, M. S. N., et al. 2019, A&A, 627, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. McKee, C. F., & Tan, J. C. 2002, Nature, 416, 59 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  73. Molinari, S., Pezzuto, S., Cesaroni, R., et al. 2008, A&A, 481, 345 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Molinari, S., Merello, M., Elia, D., et al. 2016, ApJ, 826, L8 [Google Scholar]
  75. Motte, F., Bontemps, S., Schilke, P., et al. 2007, A&A, 476, 1243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Motte, F., Bontemps, S., & Louvet, F. 2018, ARA&A, 56, 41 [Google Scholar]
  77. Mottram, J. C., Hoare, M. G., Davies, B., et al. 2011, ApJ, 730, L33 [Google Scholar]
  78. Mouschovias, T. C., Tassis, K., & Kunz, M. W. 2006, ApJ, 646, 1043 [NASA ADS] [CrossRef] [Google Scholar]
  79. Mueller, K. E., Shirley, Y. L., Evans, Neal J., I., & Jacobson, H. R. 2002, ApJS, 143, 469 [Google Scholar]
  80. Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Nakano, T., Hasegawa, T., & Norman, C. 1995, ApJ, 450, 183 [NASA ADS] [CrossRef] [Google Scholar]
  82. Nomura, H., & Millar, T. J. 2004, A&A, 414, 409 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
  83. Padoan, P., Pan, L., Juvela, M., Haugbølle, T., & Nordlund, Å. 2020, ApJ, 900, 82 [Google Scholar]
  84. Penteado, E. M., Walsh, C., & Cuppen, H. M. 2017, ApJ, 844, 71 [NASA ADS] [CrossRef] [Google Scholar]
  85. Peretto, N., Fuller, G. A., Duarte-Cabral, A., et al. 2013, A&A, 555, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Peters, T., Banerjee, R., Klessen, R. S., & Mac Low, M.-M. 2011, ApJ, 729, 72 [Google Scholar]
  87. Reboussin, L., Wakelam, V., Guilloteau, S., & Hersant, F. 2014, MNRAS, 440, 3557 [NASA ADS] [CrossRef] [Google Scholar]
  88. Rimola, A., Taquet, V., Ugliengo, P., Balucani, N., & Ceccarelli, C. 2014, A&A, 572, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Russeil, D., Zavagno, A., Motte, F., et al. 2010, A&A, 515, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Sabatini, G., Giannetti, A., Bovino, S., et al. 2019, MNRAS, 490, 4489 [Google Scholar]
  91. Sabatini, G., Bovino, S., Giannetti, A., et al. 2020, A&A, 644, A34 [EDP Sciences] [Google Scholar]
  92. Sanna, A., Kölligan, A., Moscadelli, L., et al. 2019, A&A, 623, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Saraceno, P., Andre, P., Ceccarelli, C., Griffin, M., & Molinari, S. 1996, A&A, 309, 827 [Google Scholar]
  94. Schneider, N., Csengeri, T., Bontemps, S., et al. 2010, A&A, 520, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Schneider, N., Ossenkopf, V., Csengeri, T., et al. 2015, A&A, 575, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Schuller, F., Menten, K. M., Contreras, Y., et al. 2009, A&A, 504, 415 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Semenov, D., Hersant, F., Wakelam, V., et al. 2010, A&A, 522, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Smartt, S. J. 2009, ARA&A, 47, 63 [NASA ADS] [CrossRef] [Google Scholar]
  99. Spitzer, L. 1978, Physical Processes in the Interstellar (New York: Wiley-) [Google Scholar]
  100. Spitzer, Lyman, J., & Tomasko, M. G. 1968, ApJ, 152, 971 [Google Scholar]
  101. Stahler, S. W., & Palla, F. 2005, The Formation of Stars (Hoboken: John Wiley & Sons, Ltd), 865 [Google Scholar]
  102. Tafalla, M., Myers, P. C., Caselli, P., Walmsley, C. M., & Comito, C. 2002, ApJ, 569, 815 [Google Scholar]
  103. Tafalla, M., Myers, P. C., Caselli, P., & Walmsley, C. M. 2004, Ap&SS, 292, 347 [NASA ADS] [CrossRef] [Google Scholar]
  104. Tang, X. D., Henkel, C., Wyrowski, F., et al. 2018, A&A, 611, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. Taniguchi, K., Saito, M., Sridharan, T. K., & Minamidani, T. 2018, ApJ, 854, 133 [NASA ADS] [CrossRef] [Google Scholar]
  106. Tielens, A. G. G. M. 2010, The Physics and Chemistry of the Interstellar Medium, ed. A. G. G. M. Tielens, (Cambridge, UK: Cambridge University Press) [Google Scholar]
  107. Tigé, J., Motte, F., Russeil, D., et al. 2017, A&A, 602, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Urquhart, J. S., Figura, C. C., Moore, T. J. T., et al. 2013, MNRAS, 437, 1791 [NASA ADS] [CrossRef] [Google Scholar]
  109. Urquhart, J. S., Moore, T. J. T., Csengeri, T., et al. 2014, MNRAS, 443, 1555 [Google Scholar]
  110. Urquhart, J. S., König, C., Giannetti, A., et al. 2018, MNRAS, 473, 1059 [NASA ADS] [CrossRef] [Google Scholar]
  111. Urquhart, J. S., Figura, C., Wyrowski, F., et al. 2019, MNRAS, 484, 4444 [Google Scholar]
  112. van der Tak, F. F. S. 2005, in Massive Star Birth: A Crossroads of Astrophysics, eds. R. Cesaroni, M. Felli, E. Churchwell, & M. Walmsley (Cambridge: Cambridge University Press), 227, 70 [Google Scholar]
  113. van der Tak, F. F. S., & van Dishoeck, E. F. 2000, A&A, 358, L79 [NASA ADS] [Google Scholar]
  114. van Dishoeck, E. F., & Blake, G. A. 1998, ARA&A, 36, 317 [NASA ADS] [CrossRef] [Google Scholar]
  115. Vasyunin, A. I., & Herbst, E. 2013, ApJ, 769, 34 [Google Scholar]
  116. Vasyunin, A. I., Caselli, P., Dulieu, F., & Jiménez-Serra, I. 2017, ApJ, 842, 33 [Google Scholar]
  117. Viti, S., & Williams, D. A. 1999, MNRAS, 310, 517 [NASA ADS] [CrossRef] [Google Scholar]
  118. Viti, S., Caselli, P., Hartquist, T. W., & Williams, D. A. 2001, A&A, 370, 1017 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  119. Viti, S., Collings, M. P., Dever, J. W., McCoustra, M. R. S., & Williams, D. A. 2004, MNRAS, 354, 1141 [NASA ADS] [CrossRef] [Google Scholar]
  120. Wakelam, V., Loison, J. C., Mereau, R., & Ruaud, M. 2017, Mol. Astrophys., 6, 22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  121. Wang, K. S., Bourke, T. L., Hogerheijde, M. R., et al. 2013, A&A, 558, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  122. Wienen, M., Wyrowski, F., Menten, K. M., et al. 2015, A&A, 579, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  123. Wood, D. O. S., & Churchwell, E. 1989, ApJ, 340, 265 [Google Scholar]
  124. Woon, D. E. 2002, ApJ, 569, 541 [Google Scholar]
  125. Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181 [NASA ADS] [CrossRef] [Google Scholar]
  126. Wootten, A., & Thompson, A. R. 2009, IEEE Proceedings, 97, 1463 [NASA ADS] [CrossRef] [Google Scholar]
  127. Wyrowski, F., Güsten, R., Menten, K. M., et al. 2016, A&A, 585, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  128. Zhang, Q., Ho, P. T. P., & Ohashi, N. 1998, ApJ, 494, 636 [NASA ADS] [CrossRef] [Google Scholar]
  129. Zhu, H., Tian,W., Li, A., & Zhang, M. 2017, MNRAS, 471, 3494 [NASA ADS] [CrossRef] [Google Scholar]
  130. Zinnecker, H., & Yorke, H. W. 2007, ARA&A, 45, 481 [NASA ADS] [CrossRef] [Google Scholar]
1

Atacama Pathfinder EXperiment 12 meter submillimeter telescope (Güsten et al. 2006).

3

We include thermal and cosmic-ray induced desorption.

4

Ohio StateUniversity (OSU) chemical network, version March 2008 (http://www.mpia.de/homes/semenov/Chemistry_benchmark/model.html)

6

This step was done using the fft2 function of numpy.fft: https://numpy.org/doc/stable/reference/generated/numpy.fft.fft2.html

7

The gas-phase processes included are ion-molecule, neutral-neutral, charge exchange, radiative association, radiative recombination and dissociative recombination.

All Tables

Table 1

Parameter space explored in our models; i.e. phase I (Sect. 3.1.1) and phase II (Sect. 3.1.2).

In the text
Table 2

Summary of the total masses of each model.

In the text
Table 3

Summary of fiducial initial elemental abundances, ni, with respect to the abundance of H-nuclei, nH.

In the text
Table A.1

Summary of the reactions in the chemical network.

In the text
Table A.2

Complete list of binding energies.

In the text
Table B.1

Summary of fiducial initial elemental abundances set for the TMC1 model.

In the text
Table D.1

Summary of the observed properties of the selected TOP100 sources.

In the text

All Figures

thumbnail Fig. 1

Sketch of the physical model employed in this work. Left panel: collapse phase (see Sect. 3.1.1) solved in a single-zone approximation. The chemical output of phase I is used as initial condition of phase II. Central panel: warm-up phase solved in a 1D approximation from 1 to 105 AU (see Sect. 3.1.2). Right panel: post-processing applied to compare the final chemical outputs of phase II with the column densities observed in the TOP100 (see Sects. 4.1 and 4.2). The physical parameters assumed in phases I and II are summarised in Table 1. The time evolution of the model is indicated by the dotted blue arrows.
In the text
thumbnail Fig. 2

Dust temperature radial profile at different times for an arbitrary model with the gas number density profile described by Eq. (2) and Rc = 104.5 AU and nH(r0) = 107 cm−3. With time, Td increases, driven by the protostar luminosity. The flat inner part represents the sublimation of the dust grains. The minimum temperature is set to 15 K to ensure continuity with the collapse phase.
In the text
thumbnail Fig. 3

Temporal and radial evolution of the abundances of the tracers observed in the TOP100 (Giannetti et al. 2017; Tang et al. 2018) during phase II in the gas phase (left panels) and on dust (right panels) for the same reference model as in Fig. 2. White contours indicate the temperature computed with MOCASSIN. Red curves corresponds to the evaporation front of the given tracer.
In the text
thumbnail Fig. 4

Panel a: exampleof the LOS gas-phase column density profiles obtained at the end of the post-processing described in Sect. 4.1 for models with b = 1 and Ṁ = 10−3M⊙ yr−1 and for an arbitrarily chosen source at a distance of 3.5 kpc. For the U_W items in the legend, U and W are Rc (colours) and nH(r0) (styles), respectively. Panel b: radial distribution of the gas-phase column densities extracted at the positions A, B, and C in panel a, i.e. the same times as in Fig. 2.
In the text
thumbnail Fig. 5

Panel a: example of the matrix calculated for G353.41-0.36 in all the models with Ṁ = 10−3M⊙ yr−1. Models on the y-axis are sorted from the top to the bottom by the difference between the mass of the models (see Table 2) and the value associated in the TOP100 (see Table D.1). For the sake of clarity, the colour bar is limited to five orders of magnitude and is normalised to maximum likelihood (). Panel b: zoom on the age limits identified by PL and PM showing howthe likelihood would appear if PC in Eq. (4) were not taken into account. Cyan vertical lines indicate the range of time that corresponds to the uncertainties on the observed luminosity in the TOP100, and the horizontal green line shows the same limit associated with the clump mass. Times and are shown in log-scale.
In the text
thumbnail Fig. 6

Summary of the final durations estimated in each evolutionary class defined for the TOP100. The relative number of massive clumps in each phase and observed in ATLASGAL are shown in black (“O”; Urquhart et al. 2018). Purple shaded areas indicate the 5th and 95th percentiles of the age factor, fA, and the corresponding value in yr in our sample (model “M1”), and purple markers (‘–’) represent their medians. The numbers associated with the model “M2 ” (orange), indicate the same quantities obtained by separating the models with respect to Ṁ before calculating the average age of each source (see text). The different markers in M2 show the combination of Ṁ that best match the observed relative number of object in each phase (legend).
In the text
thumbnail Fig. 7

Comparison of the median column densities observed in the TOP100 (solid lines) with those predicted by our models (circles) that were obtained by applying the same procedure as was used to quantify Δtphase (see Sect. 4.2). Panels are separated by tracers (colours). The shape of the shaded areas indicates the minimum and maximum observed column density, and the error bars associated with each circle incorporate the uncertainty of one order of magnitude assumed in Eq. (6) for the comparison. The modelled column densities of the hot component of methanol (circles in panel f) are multiplied by a factor 1000. Triangles indicate the upper limits. The column densities are shown in log-scale.
In the text
thumbnail Fig. B.1

Comparison of the time-dependent variation in the abundances of four arbitrary species (i.e. C, H2 O, CO, and CH3OH) in gas phase(left panels) and on dust (right panels) in the TMC1 model. The green line shows the results obtained with the KROME -package, and the blue crosses are the results of Semenov et al. (2010) with ALCHEMIC.
In the text
thumbnail Fig. B.2

Fractional abundance evolution for H2CO, CH3 CN, CH3 CCH, and CH3 OH as a functionof the central density of the clump during the collapse phase. Colours indicate the value of b assumed for the collapse described by Eq. (1). Solid lines show the results employing the canonical binding energies from KIDA, and shaded areas represent the solutions when Tb was varied by ±10%.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.