© The Authors 2025

1 Introduction

High-mass stars are important for the energy generation and the life cycle of the interstellar medium (ISM) in the Universe through radiation, outflows, stellar wind, and supernova events (Motte et al. 2018). Thus, it is crucial to track the evolution of high-mass stars throughout their whole lifetime. However, currently, high-mass star formation is still unclear compared with low-mass star formation, because of the difficulties in observation induced by the small number of targets (∼3%; Kroupa 2001; Wang et al. 2023), large distances from Earth (several kiloparsecs; Urquhart et al. 2018; Gieser et al. 2021), and typically short evolutionary timescales (the order of 10⁵ yr; Mottram et al. 2011; Kuiper & Hosokawa 2018; Motte et al. 2018; Gieser et al. 2021). Nevertheless, based on the observed properties, the evolution of high-mass star-forming regions (HMSFRs) can be roughly divided into four distinct stages, including the high- mass starless core (HMSC) stage, the high-mass protostellar object (HMPO) stage, the hot molecular core (HMC) stage, and the ultracompact HII (UCHII) stage, respectively (e.g., Beuther et al. 2007; Zinnecker & Yorke 2007; Gerner et al. 2014, 2015; Jin et al. 2015; Gieser et al. 2021). The HMSCs are molecular cores of cold dense gas and dust with temperatures of 10–20 K (Pillai et al. 2006; Zinnecker & Yorke 2007; Gieser et al. 2021), which are close to isothermal and radiate mostly at millimeter/submillimeter wavelengths; they commonly belong to infrared dark clouds (IRDCs; Sridharan et al. 2005; Rathborne et al. 2006). Subsequently, HMPOs are formed: gravitational collapse generates protostars that are still accreting mass (>8 M_⊙; Beuther et al. 2007; Gerner et al. 2014, 2015), and that can radiate at mid-infrared wavelengths and heat the surrounding environment to become non-isothermal (Beuther et al. 2002; Beltrán & de Wit 2016). Next, a higher temperature (>100 K; Rathborne et al. 2006; Stéphan et al. 2018; Gieser et al. 2021) prompts the formation of HMCs, accompanied by the increased gaseous chemical complexity of various molecules but similar physical conditions compared to HMPOs (Gerner et al. 2014, 2015). Finally, UCHII regions are created by strong radiation from protostars, and are associated with ionized gas, destroyed complex molecules, and free-free emission at centimeter wavelengths (Peters et al. 2010; Sánchez-Monge et al. 2013; Urquhart et al. 2013; Stéphan et al. 2018; Gieser et al. 2021). Due to the smooth evolution, these stages cannot be clearly distinguished, particularly among HMPO, HMC, and UCHII stages (Gerner et al. 2014, 2015; Gieser et al. 2021).

Chemical evolution is a useful tool for further analyzing high-mass star formation. Diverse chemical compositions across various sources can be used to infer physical properties and subsequently constrain different evolutionary conditions (e.g., Herbst & van Dishoeck 2009; Gerner et al. 2014; Jørgensen et al. 2020; McGuire 2022). The column density ratios or integrated intensity ratios between different species are usually used to distinguish different evolutionary stages and constrain lifetimes of HMSFRs, which are defined as chemical clocks (Wakelam et al. 2004; Hoq et al. 2013; Gerner et al. 2015; Yu & Wang 2015; Urquhart et al. 2019; Wang et al. 2023). Such ratios commonly increase or decrease monotonically along with the entire evolution of HMSFRs. For instance, Hoq et al. (2013) derived that the relative abundance ratio of N₂H⁺/HCO⁺ increased slightly throughout the evolution according to observed 333 sources belonging to the Millimeter Astronomy Legacy Team 90 GHz (MALT90) survey. Ohashi et al. (2014) and Tatematsu et al. (2014) proposed that the column density ratios of N(c-C₃H₂)/N(C₂S), N(NH₃ )/N(C₂S), N(NH₃)/N (HC₃N), and N(N₂H⁺)/N(C₂S) in starless regions could be lower than those in star-forming regions based on the observations of six molecular clouds in Orion A. Jin et al. (2015) estimated that the abundance ratio of HCN/HNC increased from 0.97 in quiescent IRDC cores, then to 2.65 in active IRDC cores, next to 4.17 in HMPOs, and finally to 8.96 in UCHII regions among a total of 78 sources. Yu & Wang (2015) suggest that the abundance ratios of N₂H⁺/H¹³CO⁺ and C₂H/H¹³CO⁺ decrease from massive young stellar objects (MYSOs) to UCHII regions toward 31 clumps from the MALT90 survey. For the statistics of larger samples, based on 3246 high-mass clumps from the MALT90 survey, Rathborne et al. (2016) found that the integrated intensity ratios of HCO⁺/HNC, HCN/HNC, HCO⁺/H¹³CO⁺, HNC/HN¹³C, N₂H⁺/H¹³CO⁺, and H¹³CO⁺/SiO increase along with the evolution, while the ratios of HCO⁺/N₂H⁺, HNC/N₂H⁺, and C₂H/N₂ H⁺ decrease from quiescent clumps to protostars then increase into UCHII regions. Similarly, Urquhart et al. (2019) propose 13 integrated intensity ratios as candidate chemical clocks among 570 sources from the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL), and that H¹³CN/N₂H⁺, ¹³CS/N₂H⁺ , HCN/HNC, HCN/N₂H⁺, HC₃N/HN¹³C, H¹³CN/HN¹³C, HC₃N/N₂H⁺, C₂H/N₂H⁺ , HCO⁺/N₂H⁺, and C₂H/c-C₃H₂ increase monotonically, while H¹³CO⁺/N₂H⁺, H¹³CO⁺/HN¹³C, and C₂H/HN¹³C also exhibit a tendency to increase, except in quiescent clumps. Considering the statistical significance of observed samples and the deviation degree of ratios across different stages, the most promising chemical clocks could include HCO⁺/HNC, HCN/HNC, ¹³CS/N₂H⁺, HCN/N₂H⁺, HC₃N/HN¹³C, and HC₃N/N₂H⁺ . Among these chemical clocks, N₂H⁺, HCO⁺, HCN, HNC, and C₂H are primarily involved. The two ions are mainly affected by the intensities of ionization and radiation field (e.g., cosmic-ray ionization and dissociation, photoionization, and photodissociation), while the other three molecules are mainly influenced by the temperature and the density (e.g., neutral-neutral reaction, ion-molecule reaction, accretion, and desorption). Therefore, these species can be utilized to research the evolution of HMSFRs.

Recently, Wang et al. (2023) further proposed a candidate chemical clock, N(HC₃N)/N(N₂H⁺) (hereafter HC₃N/N₂H⁺) in HMSFRs, following the integrated intensity ratio of HC₃N/N₂H⁺ proposed by Urquhart et al. (2019). According to their observations concerning 61 UCHII regions via the Insti- tut de Radioastronomie Millimétrique (IRAM) 30 m and the Arizona Radio Observatory (ARO) 12 m telescopes, combined with previous observations toward 17 HMSCs and 28 HMPOs through the Nobeyama 45 m telescope performed by Taniguchi et al. (2019b), the column density ratio increased by a factor of five from ~0.136 at the HMSC stage to ~0.303 at the HMPO stage, and finally to ~0.777 at the UCHII stage (see Table 3). The column density ratio at the HMSC stage is comparable to the integrated intensity ratio (~0.1). However, at the HMPO and UCHII stage, the column density ratio is higher than the integrated intensity ratio by a factor of two or three (~0.13 and ~0.26, respectively; see Fig. 20 in Urquhart et al. 2019). Thus, the column density ratio demonstrates a more distinct increasing trend than the integrated intensity ratio. The properties of these observed sources, including heliocentric distances, temperatures, masses, and the column densities of HC₃N and N₂H⁺ , have been provided by Taniguchi et al. (2019b), Wang et al. (2023), and relevant references, which exclude potential biases. However, some observed ratios between HMSC and HMPO stages or between HMPO and UCHII stages can overlap (see Fig. 10 from Taniguchi et al. 2019b and Fig. 6 from Wang et al. 2023), which implies that a particular source may not exhibit a strictly increasing ratio of HC₃N/N₂H⁺ throughout the evolution. Besides, the observed HC₃N/N₂H⁺ of different sources in each evolutionary stage can fluctuate by over one order of magnitude. This suggests that observed sources in the same stage may span a considerable duration of time, and the value of HC₃N/N₂H⁺ during a certain period in the HMSC, HMPO, and UCHII stage could present an increasing trend. Thus, it is necessary to constrain the potential timescales of each stage by fitting the median value of observed HC₃N/N₂H⁺ via chemical simulations.

Gerner et al. (2014), (2015) constrained the physical conditions and timescales of each stage in a 1D approximation by fitting the abundances of 18 species (including four deuterated species) toward 59 HMSFRs (19 IRDCs, 20 HMPOs, 11 HMCs, and 9 UCHII regions) through the IRAM 30 m, which can be considered as fiducial physical models in chemical simulations for further comparing other observed ratios during each evolutionary stage. Based on such simulations and observations, 15 column density ratios between different species (excluding deuterated species) were also compared. Only CS/SO, HCN/HNCO, HCO⁺/N₂H⁺, HNCO/C₂H, and N₂H⁺/C₂H show monotonic evolutionary trends. However, none of these simulated ratios fit the observed ratios well; indeed, they even exhibit opposite trends, such as HCO⁺/N₂H⁺ and N₂H⁺/C₂H. Additionally, such fiducial models had been further adopted for fitting observed abundances and estimating the lifetime of each HMSFR (e.g., Gieser et al. 2021; Yang et al. 2023). Taniguchi et al. (2019a, 2021, 2023) adopted homogeneous one-point physical models that include a free-fall collapse phase, a warm-up phase, and a hot phase (Garrod & Herbst 2006) to calculate the ratios between several complex organic molecules (COMs; HC₅N/CH₃OH, C₂H/HC₅N, HC₅N/CH₃CN, CH₃OH/CH₃CN, CH₃OCH₃/CH₃CN, HNCO/CH₃CN, CH₂CHCN/CH₃CN, and CH₃CH₂CN/CH₃CN) toward several MYSOs. These modeled ratios can constrain a hot core stage (T > 100 K) but show no clear evolutionary trend, while the one-point physical models were not sufficient to accurately describe HMSFRs. Thus, a systematic theoretical search for chemical clocks via chemical modeling is necessary.

In this paper, we perform chemical simulations to replicate the observed column density ratio of HC₃N/N₂H⁺ and the abundance of these two species across various evolutionary stages of HMSFRs. We focus on this promising chemical clock candidate because it has not been explicitly considered as an evolutionary indicator and analyzed in previous simulations. Simultaneously, we identify the chemical processes responsible for the observed time-dependent trends in these stages. Besides, we analyze other candidate chemical clocks according to numerical simulations. In Sect. 2, we describe adopted chemical models and several physical models for simulating chemical evolution within the range of 10⁵ AU (~0.5pc) in HMSFRs. The 1D models presented by Gerner et al. (2014), (2015) show fluctuations in density and temperature across the entire evolution. To maintain the steady increase in density and temperature over time as predictions of the global hierarchical collapse scenario, we also adjust parameters such as density, temperature, and time spent in each evolutionary stage. In Sect. 3, we compare the modeled average ratio of HC₃N/N₂H⁺ within different radii with observed ratios, and present the evolution of HC₃N/N₂H⁺ at different locations. In Sect. 4, we discuss key mechanisms affecting the evolution of HC₃N/N₂H⁺ , and propose other candidate chemical clocks in HMSFRs according to our simulations. Finally, in Sect. 5 we present our conclusions.

2 Model

2.1 Chemical model

To compare the observations with theoretical predictions in astrochemistry, the evolutions of HC₃N and N₂H⁺ were calculated via the astrochemical code Nautilus (Ruaud et al. 2016). In the simulations, both the two-phase (gas and grain surface) model (Hasegawa et al. 1992) and the three-phase (gas, grain surface, and grain mantle) model (Hasegawa & Herbst 1993) were adopted. Since the results adopting the latter one deviate from the observed HC₃N/N₂H⁺ because most species in the grain mantles are chemically inert when the temperature is low (<50 K), only the simulations using the two-phase model are discussed. The reaction network is identical to the one proposed by Wang et al. (2021) encompassing 9710 reactions and 947 species, which already includes the photodesorption mechanism (Öberg et al. 2009; Chang & Herbst 2014). To exclude potential bugs, this reaction network was benchmarked against the reaction network from Ruaud et al. (2016) in the Nautilus package (see Appendix A), which includes 717 species and over 11 000 reactions. A higher cosmic ray ionization rate of ζ_CR = 1.8 × 10⁻¹⁶ s⁻¹ was utilized according to the survey of 20 HMSFRs in the Galactic disk performed by Indriolo et al. (2015); according to these authors, ζ_CR is supposed to remain constant at a Galactic radius >5 kpc. Given the distances of observed sources from Taniguchi et al. (2019b) and Wang et al. (2023), it is reasonable to adopt this value, which is one order of magnitude higher than the common standard value of 1.3 × 10⁻¹⁷ s⁻¹ for the ISM. This rate is assumed to remain constant throughout the evolution, following the chemical simulations conducted by Gieser et al. (2021). It should be noted that cosmic rays are attenuated inside molecular clouds as a function of column density (Padovani et al. 2022). However, for the inner compact regions of HMSFRs with higher densities and temperatures, several chemical mechanisms become more significant compared to those in the outer diffuse areas. For example, the reaction rates of gaseous reactions and thermal desorption can be several orders of magnitude higher than those in the outer regions. This makes the cosmic ray ionization and dissociation less critical. We compared the results of HC₃N and N₂H⁺, adopting ζ_CR = 1.8 × 10⁻¹⁶s⁻¹ and 1.3 × 10⁻¹⁷ s⁻¹, respectively, using the physical model proposed by Gerner et al. (2014, see our Sect. 2.2). For the inner compact regions, the deviations between these two simulations throughout the entire evolution are fewer than several times. Therefore, for simplicity, a constant ionization rate across HMSFRs was adopted. Similarly, given the potentially stronger FUV radiation field in HMSFRs, χ was assumed to be 10 times higher than the standard value of 1, just like ζ_CR. Actually, considering the high value of visual extinction, A_V (> 10 mag), in the majority of HMSFRs, the impact of such an enhanced radiation field on the results can be ignored. We also adopted the sticking coefficient of adsorption on the grains of P_s = 1, the reactive desorption efficiency of a_RRK = 0.01 (Garrod et al. 2007; Vasyunin & Herbst 2013), the grain surface diffusion barrier of E_b1 = 0.5E_d, and the grain mantle diffusion barrier of E_b2 = 0.7Ed (Chang & Herbst 2016), where E_d is the desorption energy of each species. Other parameters can be consulted in Table 1 presented by Wang et al. (2019). In addition, the initial abundances listed in Table 1 approximate the low-metal environment at the start of the HMSC stage (Gerner et al. 2014).

2.2 Physical model

We adopted the 1D physical model of HMSFRs including four evolutionary stages derived by Gerner et al. (2014), (2015), and their parameters are presented in Table 2, designated as model 1 and model 2, respectively, serving as benchmark physical models. The densities and the temperatures at different radii were calculated with modified power laws, $$\rho (r)={\rho }_{\text{in}}{\left(r/r_{\text{in}}\right)}^{-p},r\ge r_{\text{in}};$$(1) $$\rho (r)={\rho }_{\text{in}},r<r_{\text{in}};$$(2)

and $$T(r)=T_{\text{in}}{\left(r/r_{\text{in}}\right)}^{-q},r\ge r_{\text{in}};$$(3) $$T(r)=T_{\text{in}},r<r_{\text{in}},$$(4)

which remain constant during each stage. For simplicity, the gas temperature and the grain temperature were assumed to be the same. The UCHII lifetime in model 2 was set to be 1.65 × 10⁴ yr instead of 3.0 × 10³ yr for consistency with model 1 and an easy analysis with t_total = 1.0 × 10⁵ yr. However, simulations adopting these two models cannot fit the observed evolutionary trend of HC₃N/N₂H⁺. To improve the fit between simulations and observations without altering the density distribution or the temperature distribution, the timescales of each stage were adjusted considering an uncertainty of a factor of 2–3 for the best-fit ages proposed by (Gerner et al. 2014, 2015). Consequently, models 3 and 4 only modify the duration of each stage to be 3 × 10⁴ yr compared to model 1 and 2, respectively, with prolonged timescales for the HMSC and UCHII stage, but slightly reduced timescales for the HMPO and HMC stage. Such HMSC and HMPO timescales are consistent with the statistical lifetimes for clumps with masses of ~6.3 × 10³ M_⊙ (see Table 6 from Urquhart et al. 2018), while HMC and UCHII timescales are shorter by a factor of 3–4. Nevertheless, according to the simulated evolutionary trend of HC₃N/N₂H⁺ in the HMC and UCHII stage, if we adopt longer HMC and UCHII timescales (e.g., 10⁵ yr), the ratio could be much lower than observed values in UCHII regions. Therefore, the suggested durations of the HMC and UCHII stage remain reasonable nonetheless. Additionally, ρ_in does not consistently increase throughout the evolution accompanied by the variations in r_in in both models 1 and 2. To maintain the same r_in with reasonable increasing ρ_in and T_in throughout the entire evolution, model 5 combines the parameters from the HMSC and HMC stage of model 3 with those from the HMPO and UCHII stage of model 4. A higher T_in in the HMPO stage was adopted to achieve a better agreement with observations. Finally, to have a consistently increased density and temperature within the outer regions (r ~ r_out) during the late stages, model 6 adopts reduced values of p and q in the HMC and UCHII stage compared to model 5. Similarly, other test models adopting different values of p (1.5–2.0) and q (0.3–0.4) based on model 5 were calculated but are no longer discussed in this paper, since these results cannot fit the observations well. The visual extinction was calculated by using $A_{\text{V}}(r)=A_{{\text{V}}_0}{\left(n_{\text{H}}(r)/n_{{\text{H}}_0}\right)}^{2/3}$ with $A_{{\text{V}}_0}=5\text{\hspace{0.17em}mag}$ and $n_{{\text{H}}_0}={10}^4{\text{cm}}^{-3}$ (Taniguchi et al. 2019a), while n_H(r) = 2ρ(r). For simplicity, the physical parameters undergo an instantaneous variation at the transition point between two stages, as was proposed by (Gerner et al. 2014, 2015). For each physical model, we selected about 45 different locations from the center to the outmost radius of 10⁵ AU (~0.5 pc) for chemical modeling, then the abundances at different radii and the average abundances within different ranges in an HMSFR could be computed. The average number density of each species is equal to the summed molecule number in each shell divided by the total volume of the same shell; at this point, the average abundance with respect to H₂ can be calculated.

Fig. 1 Correlation between the average ratio of HC3N/N2H+ and the average abundance of HC3N with respect to H2 within four different ranges (panel A: 7 × 103 AU, panel B: 104 AU, panel C: 5 × 104 AU, and panel D: 105 AU) around the center by adopting model 4 in chemical simulations compared with observations. The results during the HMSC, HMPO, HMC, and UCHII stages are plotted as squares, circles, diamonds, and triangles, respectively. The evolution of HC3N/N2H+ and HC3N/H2 throughout ttotal during each stage consistently decrease (from right to left and from top to bottom). The observations proposed by Taniguchi et al. (2019b) and Wang et al. (2023, including derived column density of HC3N, N2H+, and H2 simultaneously) are plotted as solid squares, solid circles, and solid triangles, respectively, with error bars.

Fig. 2 Same as Fig. 1, but adopting model 6.

3 Results

3.1 Comparison of chemical simulations and observations

The sources observed by Taniguchi et al. (2019b) and Wang et al. (2023) are mainly located at a distance of ~3–10 kpc. Considering the largest 29″ beam size of the IRAM 30 m, the linear resolution corresponds to ~0.4–1.4pc. Since r_out = 10⁵ AU is comparable to the linear resolution, and also approximately corresponds to the outermost radius of dense parts of HMSFRs (Gerner et al. 2014), the average abundances of HC₃N and N₂H⁺ with respect to H₂ within a radius of 10⁵ AU were calculated, and the corresponding ratio HC₃N/N₂H⁺ was compared with observed targets (Taniguchi et al. 2019b; Wang et al. 2023). To investigate the correlation between the simulated average ratio of HC₃N/N₂H⁺ and HC₃N/H₂ (the abundance of HC₃N with respect to H₂) and compare them with observations, we plot them as scatter diagrams in Figs. 1 and 2, within four different ranges (7 × 10³ AU, 10⁴ AU, 5 × 10⁴ AU, and 10⁵ AU) around the center using models 4 and 6, respectively. These figures show that there is only one radius per model (10⁴ AU for model 4 and 10⁵ AU for model 6) that can reproduce three of the four stages. The detailed comparison between models and observations is discussed below. We note already that due to a lack of observational constraints and observational ambiguities, we did not consider the HMC stage for comparison with our models. Such different ranges were chosen to verify whether the average ratios within certain ranges can fit the observations. The average ratios within 7 × 10³ AU usually severely deviate from the observations. The evolution during the HMSC stage, HMPO stage, HMC stage, and UCHII stage is depicted with different symbols, indicating that the correlation consistently evolves from right to left with a decreasing HC₃N/N₂H⁺ ratio, while it moves from top to bottom with a decreasing HC₃N/H₂ ratio in each panel. A line depicting the evolutionary trend is omitted in each panel for concise pictures, since abrupt jumps are induced by discontinuous changes in the physical parameters between the two stages. However, such abrupt jumps do not affect the comparison between modeled and observed ratios. The average ratios within 10⁴ AU around the center can fit the observations at the final HMSC stage (t_total ~ 2.8–3.0 × 10⁴ yr), the start HMPO stage (t_total ~ 3.0–3.1 × 10⁴ yr), and the final UCHII stage (t_total ~ 1.03–1.2 × 10⁵ yr), respectively (see panel B in Fig. 1). Nevertheless, on the one hand, it is noteworthy that a strictly monotonic increase in HC₃N/N₂H⁺ over the entire evolution cannot be achieved. On the other hand, such a region only accounts for 10% of the radius and 10⁻³ of the volume compared to the entire area of HMSFRs, which seems unlikely to reflect the actual conditions. The simulations within the other three ranges cannot fit the observations, with a lower HC₃N/H₂ ratio in the HMPO stage and a higher HC₃N/N₂H⁺ ratio in the UCHII stage within the range of 7 × 10³ AU (panel A), a higher HC₃N/N₂H⁺ ratio in the HMSC stage and a lower HC₃N/N₂H⁺ ratio in the UCHII stage within the range of 5 × 10⁴ AU (panel C), and a lower HC₃N/N₂H⁺ ratio in the UCHII stage within the range of 10⁵ AU (panel D). Additionally, for the inner region within 7 × 10³ AU, both simulated HC₃N/N₂H⁺ and HC₃N/H₂ ratios can deviate from the observations by over one order of magnitude, because of the high temperature and density. These two simulated ratios can only fit the observations at the late HMSC stage. Thus, such results are no longer plotted.

Similarly, Fig. 2 depicts the comparison between the simulations using model 6 and observations. In panel D, the average ratios within 10⁵ AU from the center can fit observations at the late HMSC stage (t_total ~ 2.0–3.0 × 10⁴ yr), the early HMPO stage (t_total ~ 3.0–3.8 × 10⁴ yr), and the early UCHII stage (t_total ~ 9.4–10.1 × 10⁴ yr), respectively. Therefore, this simulation is more reasonable than the result using model 4, although the strictly monotonic increase in HC₃N/N₂H⁺ throughout the entire evolution remains unreplicated. The results within more internal regions deviate from observations more significantly, especially during the HMPO and UCHII stage (see panels A, B, and C), due to the influence of much higher temperature and density in the innermost regions. In addition, the results adopting the other four models are presented in Appendix B, which cannot fit the observations regardless of the radius of the selected regions. Consequently, we propose that the observed increased HC₃N/N₂H⁺ ratio from the HMSC stage to the HMPO stage and to the UCHII stage in HMSFRs can be fit during a certain period in each stage in chemical simulations, which may constrain the evolution at the late HMSC stage, the early HMPO stage, and the early UCHII stage.

To individually analyze the time-dependent evolution of average ratios of HC₃N/N₂H⁺ , HC₃N/H₂, and N₂H⁺/H₂ in models 4 and 6, we depict the average ratios within the four different aforementioned ranges (7 × 10³ AU, 10⁴ AU, 5 × 10⁴ AU, and 10⁵ AU) around the center across t_total in Fig. 3. In model 4, the average ratio of HC₃N/N₂H⁺ decreases during the HMSC stage and fits the observations at the final HMSC stage, particularly the ratio within the region of 10⁴ AU. Then an abrupt increase occurs at the start HMPO stage within the innermost 10⁴ AU region, mainly induced by the desorption of HC₃N on the grains due to high temperatures. However, the ratio rapidly decreases to be lower than the observations in this region, while within larger regions (5 × 10⁴ AU and 10⁵ AU) the ratio continues to decrease from the HMSC stage. During the HMC stage, the ratio similarly increases at the start because of high temperatures but subsequently decreases within all four regions. The ratio continues to decrease during the UCHII stage, but the values within the innermost 10⁴ AU region can fit the observations, while those within larger regions (5 × 10⁴ AU and 10⁵ AU) can be lower than the observations. For the average abundance of HC₃N, the evolutionary tendency is similar to that of HC₃N/N₂H⁺, but the results within all four regions can roughly agree with the observations. The average abundance of N₂H⁺ increases during each stage across all four regions (also see Gerner et al. 2014, 2015) in contrast to the observed decreased abundance of N₂H⁺, but only the result within 10⁴ AU can fit the observations. In model 6, the average ratio of HC₃N/N₂H⁺ still exhibits a decreasing trend during each stage, with sharp increases at the start of each stage compared with the values at the end of the previous stage. The ratio within 7 × 10³ AU or 10⁴ AU is lower than the observations during the HMSC stage, but much higher than the observations during the HMPO and UCHII stage. The ratio within 5 × 10⁴ AU is higher than the observations during the HMSC stage, but can fit the observations during the HMPO and UCHII stage. Only the ratio within 10⁵ AU can fit the observations during all three stages. For the average ratio of HC₃N/H₂, the results within four regions can agree with the observations during the HMSC and UCHII stage, but only the result within 10⁵ AU can fit the observations during the HMPO stage. For the average ratio of N₂H⁺/H₂, the tendency to increase occurs during each stage accompanied by a dramatic decrease at the start of each stage compared to the values at the end of the previous stage. Still, only the result within 10⁵ AU can fit the observations during the UCHII stage. Obviously, a rigorously monotonically increasing ratio of HC₃N/N₂H⁺ throughout the evolution of HMSFRs cannot be accurately replicated in simulations.

Due to the current lack of observations of HC₃N/N₂H⁺ in HMCs, the evolution of HC₃N/N₂H⁺ in chemical simulations is temporarily not regarded as a primary focus in the analysis. However, there are UCHII regions that have the characteristic spectrum of a hot core, and HMPOs can appear as hot cores, depending on the resolution and sensitivity in observations. These observational uncertainties affect the calculation of column densities and abundances, as well as the clump classification. In fact, a scatter of around one dex in the observed ratios appears within any evolutionary stage, which may induce the ratios during the HMC stage to overlap with those in the HMPO and UCHII stages. Additionally, this ratio during the HMC stage indeed follows the increasing trend, considering the integrated intensity ratio of HC₃N/N₂H⁺ derived by Urquhart et al. (2019). Consequently, it should be higher than that during the HMPO stage and lower than that during the UCHII stage. In the simulations mentioned above, both HC₃N/N₂H⁺ and HC₃N/H₂ in the HMC stage reach higher values compared to those not only in the HMPO stage but also in the UCHII stage, which contradicts the observations. Such a discrepancy is mainly driven by increased temperatures in the larger areas among HMSFRs during the HMC stage, which induces the thermal desorption of HC₃N molecules from the grains and subsequently a rapid increase in HC₃N/N₂H⁺ (see Sect. 4.1). Nevertheless, the average HC₃N/N₂H⁺ ratio within 10⁵ AU at the late HMC stage already approaches the observed values in UCHII regions (see Fig. 3). Combined with the decreasing trend of HC₃N/N₂H⁺ in this stage, adopting a longer duration of this stage (>4 × 10⁴ yr) potentially makes the simulated ratio fit the observed values. In addition, varying temperature distributions across the radius during this stage in physical models could significantly affect the degree of the increase in HC₃N/N₂H⁺, as the scale of the area where the thermal desorption of HC₃N occurs is directly influenced. Such a comparison could be further explored in a future study, when more observations and simulations are available.

To further constrain the evolutionary timescale in each stage, we calculated the confidence level, κ, of HC₃N/N₂H⁺, HC₃N/H₂, and N₂H⁺/H₂ by adopting the “mean confidence level” method (e.g., Garrod et al. 2007; Hassel et al. 2008, 2011; Wang et al. 2019, 2021). For one specific ratio, i, the agreement between the calculated value, X_i, and the observed value, X_obs,i, is defined as the confidence level, κ_i , which is calculated by $${\kappa }_i=\mathrm{erfc}\left(\frac{\left|\lg X_i-\lg X_{\text{obs},i}\right|}{\sqrt{2}\sigma }\right),$$(5)

with the complementary error function, erfc = 1 − erf, and the standard deviation, σ = 1. A higher value of κ_i (0−1) indicates a better agreement between the calculated and observed values, while κ_i = 0.317 corresponds to the calculated value being one order of magnitude higher or lower than the observed value. Subsequently, the average κ across three ratios, HC₃N/N₂H⁺, HC₃N/H₂, and N₂H⁺/H₂, during each stage was calculated. By identifying the maximum average κ, the best-fit timescale in each stage can be obtained, which can further constrain the corresponding physical parameters.

Table 3 presents the best-fit timescales for the HMSC stage, HMPO stage, and UCHII stage within 10⁴ AU in model 4 and within 10⁵ AU in model 6, respectively, which can fit all three observed ratios. It is crucial to compare these timescales with previous theoretical studies such as (Gerner et al. 2014, 2015, models 1 and 2), which can evaluate the reliability of these timescales and other physical parameters adopted in our modified models. In the HMSC stage, the best fit occurs at t_total ~ 3 × 10⁴ yr, which may imply that a lifetime, t_HMSC = 3 × 10⁴ yr, a factor of 2–3 longer than models 1 and 2 (Gerner et al. 2014, 2015) is required to decrease the ratio of HC₃N/N₂H⁺. A quick evolution within about 10³ yr in the HMPO stage (t_total ~ 3.1 × 10⁴ yr) can fit the observed ratios, which may signify that t_HMPO = 3.0 × 10⁴ yr is sufficient compared to the longer timescale of 6.0 × 10⁴ yr in model 1 and the comparable timescale of 3.2 × 10⁴ yr in model 2. Finally, in the UCHII stage, a duration of about 6.7 × 10³ yr (t_total ~ 9.67 × 10⁴ yr) achieves the best fit, which approaches half of t_UCHII = 1.2 × 10⁴ yr in model 1 and t_UCHII = 1.65 × 10⁴ yr in model 2. Overall, both the best-fit timescale and the total duration of each stage in our models are compatible with those in models 1 and 2 with an uncertainty of a factor of 2–3. The best-fit timescales indicate the order of magnitude of the duration of the process of high-mass star formation, as the modeled chemical composition at these timescales is in broad agreement with statistical estimates. In addition, due to the relatively stable chemical evolutionary trend, certain periods can be confirmed. Different durations of the previous stage mainly influence the abundances at the start of the subsequent stage, with fluctuations by a factor of several times, but do not significantly affect the overall chemical evolutionary trend. Such differences could affect the agreement between simulations and observations at the start of subsequent stages. Overall, our modified physical models incorporating higher densities and temperatures in the outer regions still remain reasonable and provide valuable insights into the evolutionary timescales of HMSFRs, since different HMSFRs could evolve with different physical parameters.

Fig. 3 Time-dependent evolution of the average ratio of HC3N/N2H+ and the average fractional abundance of HC3N and N2H+ with respect to H2 within four different ranges (7 × 103 AU, 104 AU, 5 × 104 AU, and 105 AU) around the center using model 4 (left panels) and model 6 (right panels). The four evolutionary stages are labeled in each panel, and the observations are also plotted, with the dashed lines as the median values and the dotted lines as the maximums or minimums of observations. The evolution before 104 yr is excluded to show the changes in different stages more clearly. The results within 7 × 103 AU and 104 AU during the HMSC stage in model 4 completely overlap because of the large rin > 104 AU.

3.2 Evolution of HC₃N/N₂H⁺ at different locations

To investigate potential physical parameters that mainly dominate the evolution of HC₃N/N₂H⁺, it is necessary to analyze the agreement between observations and simulated ratios across different layers, spanning from the center to the outermost radius, 10⁵ AU, in the best-fit models. Figure 4 presents the timedependent evolution of HC₃N/N₂H⁺, HC₃N/H₂, and N₂H⁺/H₂ at six different radii (10² AU, 10³ AU, 7 × 10³ AU, 10⁴ AU, 5 × 10⁴ AU, and 10⁵ AU) from the center using model 4 and model 6. The evolutionary trends at 7 × 10³ AU, 10⁴ AU, 5 × 10⁴ AU, and 10⁵ AU in both models resemble those plotted in Fig. 3, but with deviations by a factor of 3–5. Nonetheless, the evolutions at 10⁴ AU in model 4 and at 10⁵ AU in model 6 can still roughly fit the observations at similar aforementioned timescales in each stage. In contrast, the evolutions at 10² AU and 10³ AU exhibit greater fluctuations, which can deviate from the observations by several orders of magnitude and can be neglected in our discussion.

Subsequently, we mainly present the densities, the visual extinctions, and the temperatures at large radii of 10⁴ AU and 10⁵ AU in models 4 and 6, which are summarized in Table 4. Firstly, these parameters cannot alter the overall evolutionary trends of HC₃N/N₂H⁺, HC₃N/H₂, and N₂H⁺/H₂ in these two models, especially the results at the same radius (see Fig. 4). Next, according to the chemical kinetics, the density primarily influences the overall evolutionary efficiency in the gas phase, while the visual extinction mainly affects the synthesis of ions via photoionization and photodissociation, especially at low values (e.g., <5 mag). Correspondingly, for the evolution during the HMSC stage with the same temperatures in the same model, a higher density at 10⁴ AU induces a more rapid evolution of HC₃ N/N₂ H⁺ and HC₃N/H₂, with larger fluctuations compared with those at 10⁵ AU, whereas N₂H⁺/H₂ at both locations evolves steadily. In addition, a higher temperature at 10⁴ AU or 10⁵ AU in model 6 induces lower HC₃N/N₂H⁺ and HC₃N/H₂, but higher N₂H⁺/H₂ during the HMSC stage compared with model 4 (20.9 K vs 11.3 K), accompanied by compatible densities and visual extinctions. The discrepancies of these ratios at the same locations in the two models become more distinct during the HMPO stage, as the temperatures in model 6 are higher than those in model 4 by a factor of 3 (36.5 K vs 12.2 K at 10⁴ AU and 14.5 K vs 4.8 K at 10⁵ AU). The abrupt increase in HC₃N/N₂H⁺ and HC₃N/H₂ at the start of the HMPO and HMC stage corresponds to rapidly produced HC₃N, which could be induced by the sublimation of CH₄ from grains at T ~ 25 K followed by gaseous reactions (see Sect. 4.1). This may indicate that the temperature is more significant than the density and the visual extinction to the evolution of HC₃ N/N₂ H⁺ in HMSFRs. For the two best-fit cases, the densities and the visual extinctions at 10⁴ AU in model 4 can be over one order of magnitude higher than those at 10⁵ AU in model 6, but the temperatures are similar except for the HMC stage. The density during each stage probably mainly affects the evolutionary timescale for replicating the observations. The visual extinction mainly influences the evolution of ions via photodissociation in the outer region (several 10⁴ AU – 10⁵ AU) with low A_V ~ 5 during the HMSC stage and the HMPO stage. Consequently, a relatively cold diffuse envelope throughout the entire evolution of HMSFRs is essential to replicate the observed ratio of HC₃ N/N₂ H⁺, while the inner regions characterized by much higher temperatures are not necessary. Such an analysis suggests that the cold envelope could play a crucial role in maintaining specific chemical conditions.

Fig. 4 Same as Fig. 3, but depicting the ratio of HC3N/N2H+ and the fractional abundance of HC3N and N2H+ with respect to H2 at six different radii (102 AU, 103 AU, 7 × 103 AU, 104 AU, 5 × 104 AU, and 105 AU) from the center using model 4 (left panels) and model 6 (right panels). The results at r = 102 AU, 103 AU, 7 × 103 AU, and 104 AU during the HMSC stage in model 4 completely overlap because of the large rin > 104 AU.

4 Discussion

4.1 Key mechanisms affecting HC₃N/N₂H⁺

According to derived physical conditions at the outermost radius in model 6, we discuss the key mechanisms for influencing HC₃N/N₂H⁺ in HMSFRs. The major reactions for producing HC₃N are: $${\text{HC}}_3{\text{NH}}^{+}+{\text{e}}^{-}\to {\text{HC}}_3\text{\hspace{0.17em}N}+\text{H},$$(6) $$\text{CN}+{\text{C}}_2{\text{H}}_2\to {\text{HC}}_3\text{\hspace{0.17em}N}+\text{H},$$(7) $$\text{N}+{\text{C}}_3{\text{H}}_3\to {\text{HC}}_3\text{\hspace{0.17em}N}+{\text{H}}_2.$$(8)

Meanwhile, the important reactions for destroying HC₃N are: $${\text{C}}^{+}+{\text{HC}}_3\text{\hspace{0.17em}N}\to {\text{C}}_3{\text{H}}^{+}+\text{CN},$$(9) $${\text{C}}^{+}+{\text{HC}}_3\text{\hspace{0.17em}N}\to {\text{C}}_4{\text{N}}^{+}+\text{H},$$(10) $${\text{H}}_3^{+}+{\text{HC}}_3\text{\hspace{0.17em}N}\to {\text{HC}}_3{\text{NH}}^{+}+{\text{H}}_2,$$(11) $${\text{CH}}_5^{+}+{\text{HC}}_3\text{\hspace{0.17em}N}\to {\text{HC}}_3{\text{NH}}^{+}+{\text{CH}}_4.$$(12)

Therefore, the evolution of HC₃N in HMSFRs can be explained as follows. In the relatively cold HMSC stage with T < 25 K, in the outer diffuse region (several 10⁴ AU – 10⁵ AU), the ion-molecule reactions (9), (10), and (11) mainly destroy HC₃N, while the dissociative recombination reaction (6) mainly converts the ion HC₃NH⁺ back to HC₃N. These reactions can be further enhanced by more ions (especially C⁺) and electrons induced by photodissociation due to the low AV . For the inner dense region (r < 10⁴ AU), the neutral-neutral reactions (7) and (8) are more significant for the production of HC₃N. They are triggered by carbon-chain species, C₂H₂ and C₃H₃ , produced via the reaction of C⁺ with H₂ and subsequent reactions (Sakai & Yamamoto 2013). In addition, reactions (11) and (12) play a major role in the destruction of HC₃N in the dense region, but most products can be converted back to HC₃N via reaction (6). When T > 25 K during the HMPO, HMC, and UCHII stage, reactions (7) and (8) still remain significant, but driven by the thermal desorption of CH₄ and subsequent reactions with C⁺ for synthesizing C₂H₂ and C₃H₃ (Aikawa et al. 2008; Hassel et al. 2008, 2011; Sakai & Yamamoto 2013; Taniguchi et al. 2019a), also known as warm carbon-chain chemistry (WCCC; Sakai et al. 2008, 2009; Sakai & Yamamoto 2013; Wang et al. 2019). For even hotter regions with T > 90 K, the sublimation of HC₃N molecules trapped on grains is the only dominating mechanism for increasing gaseous HC₃N molecules. Thus, at the start of the HMPO and HMC stage, the average abundance of HC₃N increases, as CH₄ and HC₃N molecules evaporate from the grains in specific regions where T becomes higher than their sublimation temperatures. During the HMPO stage, the sublimation of CH₄ and the WCCC mechanism mainly occur within the range of 10³−3 × 10⁴ AU, which is smaller than that during the HMC stage within 3 × 10³−10⁵ AU. In contrast, the sublimation of HC₃N occurs in similar regions during both stages (<10³ AU during the HMPO stage and <3 × 10³AU during the HMC stage, respectively). A larger warm region during the HMC stage results in a higher increase in HC₃N compared to that during the HMPO stage. However, reactions (9)–(12) gradually decrease the abundance of HC₃N over time during each stage, since CH₄ and HC₃ N molecules trapped on the grains are finite.

For N₂H⁺, there is only one major formation reaction, $${\text{H}}_3^{+}+{\text{N}}_2\to {\text{N}}_2{\text{H}}^{+}+{\text{H}}_2,$$(13)

and the major destruction reactions are: $${\text{N}}_2{\text{H}}^{+}+{\text{e}}^{-}\to \text{NH}+\text{N},$$(14) $${\text{N}}_2{\text{H}}^{+}+{\text{e}}^{-}\to {\text{N}}_2+\text{H},$$(15) $${\text{N}}_2{\text{H}}^{+}+\text{CO}\to {\text{HCO}}^{+}+{\text{N}}_2,$$(16) $${\text{N}}_2{\text{H}}^{+}+{\text{CH}}_4\to {\text{CH}}_5^{+}+{\text{N}}_2,$$(17) $${\text{N}}_2{\text{H}}^{+}+{\text{NH}}_3\to {\text{NH}}_4^{+}+{\text{N}}_2,$$(18) $${\text{N}}_2{\text{H}}^{+}+{\text{H}}_2\text{O}\to {\text{H}}_3{\text{O}}^{+}+{\text{N}}_2.$$(19)

During the HMSC stage, T is already higher than the sublimation temperature of N₂ (about 19 K), and sublimated N₂ molecules drive the formation of N₂H⁺ via reaction (13). In contrast, reactions (14) and (15) mainly destroy N₂H⁺, driven by intense cosmic-ray ionization and photodissociation. During the hotter HMPO, HMC, and UCHII stages, reaction (13) remains significant for forming N₂H⁺. But reactions (16), (17), (18), and (19) can sequentially become important for the destruction of N₂H⁺ when the temperature rises above ~21 K, the sublimation temperature of CO, ~25K for CH₄, ~115K for NH₃, and ~120K for H₂ O, respectively. Consequently, at the start of the HMPO, HMC, and UCHII stages, the thermal desorption of CO, CH₄ , NH₃, and H₂O from the grains enhances reactions (16)–(19) and induces the decrease in N₂H⁺ compared to that at the end of the previous stage. However, the efficiency of reaction (13) always overwhelms the total efficiency of reactions (14)-(19) during each stage, which accumulates more N₂H⁺ ions for increasing the abundance during each stage. As a result, combined with the previous analysis of HC₃N, the ratio HC₃N/N₂H⁺ can rapidly increase at the start of each stage, accompanied by a decreasing trend as the stages progress. In summary, the evolution of HC₃N is mainly affected by the WCCC mechanism and its thermal desorption, whereas the evolution of N₂H⁺ is primarily influenced by the thermal desorption of N₂, CO, CH₄, NH₃, and H₂O from the grains, followed by dissociative recombination and ion-molecule reactions.

4.2 Other candidate chemical clocks

Based on the evolution of HC₃N and N₂H⁺ during different stages of HMSFRs, carbon-chain species are affected by the WCCC mechanism and their thermal desorption just like HC₃N, while ions are mainly influenced by cosmic-ray ionization, photodissociation, and dissociative recombination reactions as N₂H⁺ . Thus, we suggest that the ratio between one carbon- chain species and an ion could be a candidate chemical clock. Other candidate chemical clocks derived from observations (see Sect. 1) can also be further analyzed with our simulations. Therefore, we examined in total 350 ratios involving 27 species, including carbon-chain species (C₂H, C₃H, C₄H, C₅H, C₂H₂, C₃H₂, C₄H₂, C₅H₂, CH₃ CCH, HC₃N, HC₅N, and C₂S), simple molecules (CS, NH₃ , HCN, HNC, CO, SO, H₂CO, and HNCO), COMs (HCOOH, CH₃OH, and CH₃CN), and ions (N₂H⁺, HCO⁺, ${\text{H}}_3^{+}$, and C⁺). Except for the aforementioned HC₃N/N₂H⁺, 178 ratios exhibit increasing or decreasing evolutionary trends at the late HMSC stage, the early HMPO stage, and the early UCHII stage according to our simulations. Figure C.1 depicts the evolution of these average ratios within 10⁵ AU from the center adopting model 6, which could be proposed as candidates of chemical clocks from the perspective of chemical evolution. Among them, 157 ratios are observable and can be verified with observations, whereas ${\text{H}}_3^{+}$ is not observable, and its related ratios cannot be considered as chemical clocks. Thus, such candidate ratios are presented as a foundational database for future observational analysis and simulations. Besides the aforementioned mechanisms affecting carbon-chain species and ions, respectively, for simple molecules, thermal desorption and ion-molecule reactions are significant, and for COMs, two-body reactions on the grains and thermal desorption are important. Nevertheless, the specific reactions affecting the evolution of each ratio can be further discussed in a future study.

According to our results, candidate chemical clocks C₂H/HCO⁺ and HC₃N/HNC show similar evolutionary trends as observations (Yu & Wang 2015; Urquhart et al. 2019), which may be further confirmed as chemical clocks. In addition, CS/SO, HCN/HNCO, C₂H/HNCO, and SO/N₂H⁺ exhibit similar evolutionary trends compared with previous simulations (Gerner et al. 2014) as candidate chemical clocks. In contrast, C₂H/N₂H⁺ and C₂H/C₃H₂ display a decreasing tendency along with the evolution, which contradicts the statistics of observations (Rathborne et al. 2016; Urquhart et al. 2019), although the former one, C₂H/N₂H⁺, is consistent with previous simulations (Gerner et al. 2014). Except for these candidate chemical clocks based on our simulations, we also depict some ratios in Fig. D.1 that have been suggested as chemical clocks in previous studies. The ratio of N₂H⁺ /HCO⁺ is indeed lower during the UCHII stage than that during previous stages and in agreement with observations (Yu & Wang 2015; Urquhart et al. 2019), but the ratio during the HMSC stage and HMPO stage does not show clear differences. For C₃H₂/C₂S, NH₃/C₂S, and N₂H⁺/CS, they increase from the HMSC stage to the HMPO stage and are consistent with observations (Ohashi et al. 2014; Tatematsu et al. 2014), while NH₃/HC₃ N decreases from the HMSC stage to the HMPO stage and contradicts the observations (Ohashi et al. 2014). However, these four ratios apparently do not exhibit monotonic evolutionary trends when one examines their values across the HMSC, HMPO, and UCHII stages. For HCN/HNC, no obvious difference appears during the whole evolution, which is not consistent with the increased ratio that is observed (Jin et al. 2015; Rathborne et al. 2016; Urquhart et al. 2019) or with previous simulations (Gerner et al. 2014). Also, the evolutionary trends of HCO⁺/HNC, HNC/N₂H⁺, HCN/N₂H⁺, and CS/N₂H⁺ contradict observations (Rathborne et al. 2016; Urquhart et al. 2019). Such deviations between simulations and observations may be driven by the simplification of the physical models adopted in the simulations. We also compared CH₃ CN/C₂ H, CH₃OH/CO, CH₃OH/C₂H, CS/C₂H, H₂CO/CO, N₂H⁺/CO, and CH₃ OH/CH₃ CN with simulations presented by Gerner et al. (2014). Among them, CH₃CN/C₂H, CS/C₂H, N₂H⁺/CO, and H₂CO/CO present roughly similar evolutionary trends across four stages, along with deviations in their values during the same stage. CH₃OH/CO, CH₃OH/C₂H, and CH₃OH/CH₃CN exhibit opposite evolutionary trends, accompanied by severe deviations in their values. Since most ratios do not exhibit monotonic evolutionary trends, it is not necessary to systematically analyze the reasons inducing the discrepancies between our simulations and results derived by Gerner et al. (2014). We propose that the discrepancies could be attributed to the different physical parameters between our best-fit model and previous models. Such a detailed analysis could be conducted in a future study.

Given the deviations between observations and simulations, some improvements could be considered in the future analysis for our understanding of HMSFRs. Firstly, the 1D physical model of HMSFRs, with constant parameters during each stage and abrupt transitions between two stages, is apparently not sufficiently accurate, although the chemistry can quickly respond to the discontinuities of physical parameters (Gerner et al. 2014). A more sophisticated physical model that incorporates varying parameters and smoother transitions between stages should be developed. An example model is described by Sabatini et al. (2021), which includes a collapse stage and a warm-up stage with continuity between them. The primary distinction between our optimal model 6 and models proposed by Gerner et al. (2014) or Sabatini et al. (2021) is that our model 6 adopts higher densities and temperatures within the range of 10⁵ AU, along with an enhanced cosmic ray ionization rate. Such a physical model could be further analyzed by adopting a larger range of parameters and supplementing the UCHII stage in a future study. Secondly, the chemical model also needs to be further modified by including more updated reactions and mechanisms, especially those related to the grain surfaces. Finally, when comparing simulations with observations, it is important to recognize that only one ratio between two species may not rigorously constrain the evolutionary stage of HMSFRs. This conclusion is mainly based on the results in panel C of Fig. 2 (also see panel C of Fig. B.3 and panel C of Fig. B.4). The average ratio of HC₃N/N₂H⁺ within 5 × 10⁴ AU from the center can also replicate similar trends and values in the HMSC, HMPO, and UCHII stages, respectively. However, the average abundance ratio of HC₃N/H₂ obviously deviates from the observations. If only one ratio between two species were compared, such a discrepancy could go unnoticed. Also, such a ratio could be replicated in different stages in simulations, accompanied by different magnitudes of the abundances of the involved species. Therefore, to rigorously constrain the evolutionary stage of HMSFRs, not only the column density ratio between two species but also the abundances of these two species should be studied. Alternatively, more than two column density ratios should be simultaneously analyzed when the column density of H₂ is not derived.

5 Conclusions

In this paper, we have performed chemical simulations to reproduce the observed column density ratio of HC₃N/N₂H⁺ and the abundances of these two species across various evolutionary stages in HMSFRs. Observations show an increasing evolutionary trend of HC₃N/N₂H⁺, which is considered a chemical clock. Simultaneously, we have identified the chemical processes responsible for the observed time-dependent trends in these stages. Our simulations adopted the astrochemical code Nautilus and the existing 1D models of HMSFRs proposed by (Gerner et al. 2014, 2015), which include four evolutionary stages with constant densities and temperatures during each stage, as well as abrupt transitions between adjacent stages. However, both density and temperature fluctuate throughout the entire evolution. To maintain a steady increase in density and temperature over time as the global hierarchical collapse scenario predicts (e.g., Vázquez-Semadeni et al. 2019), we adjusted parameters such as the density, temperature, and time spent in each stage. The average ratios ofHC₃N/N₂H⁺ within different ranges from the center were calculated and compared with observations, as were the average abundances of HC₃N and N₂H⁺ with respect to H₂. The main conclusions are listed below:

1. When averaging over large spatial scales (~10⁵AU from the core), the best model produced successfully matches the observed column density ratio of HC₃N/N₂H⁺ and the abundances of the species involved at specific times for each evolutionary stage; that is, the late HMSC stage, the early HMPO stage, and the early UCHII stage, respectively. Assuming that thermal desorption occurs instantaneously, the time required to bring the simulated and observed ratios into agreement reflects the time necessary for chemical adjustments following jumps of physical parameters between stages. The agreement is only achieved when averaging over large spatial scales, and indicates that clusterforming regions cannot be modeled without accounting for the extended envelope.

2. Our simulation fails to demonstrate a continuous and consistent monotonic increase in HC₃N/N₂H⁺ throughout the entire evolution. However, some observed ratios between adjacent stages overlap, which could be induced by observational uncertainties (such as those in deriving column densities and abundances, clump classification, and systematic effects), or indicate that the evolution of HC₃N/N₂H⁺ may not strictly monotonically increase throughout the entire evolution; for example, because a cluster is being formed in the clump, and different parcels of gas can have different physical properties and thus a different chemistry.

3. The evolution of HC₃N is mainly affected by the WCCC mechanism and its own thermal desorption, while the evolution of N₂H⁺ is primarily influenced by the thermal desorption of N₂, CO, CH₄ , NH₃ , and H₂O on the grains, followed by gaseous dissociative recombination and ion-molecule reactions.

4. We also examined 350 other ratios involving 27 species according to our best-fit model. We find that 178 ratios exhibit an increasing or decreasing evolutionary trend at the late HMSC stage, the early HMPO stage, and the early UCHII stage. Among them, 157 ratios are observable and can be compared with observations, indicating that they could be considered as candidate chemical clocks.

Overall, the results obtained from the best-fitting model timescales broadly agree with statistical estimates. This indicates that 1D models with abrupt jumps in physical parameters have reached their limits in terms of the insights they can provide. Thus, more sophisticated models and advanced approaches are necessary in future studies. For instance, physical models incorporating smoother transitions between evolutionary phases might address some of the discontinuities that cause the current models to deviate from the observational data.

Acknowledgements

We thank our referee’s constructive comments to significantly improve the quality of the manuscript. This work is supported by the National Natural Science Foundation of China (NSFC) grant No. 12303032, the Natural Science Foundation of Jiangsu Province (Grant No. BK20221163), and the Excellent Postdoctoral Talent Program of Jiangsu Province (Grant No. 2022ZB475). F.D. is supported by NSFC grant No. 12041305 and the National Key R&D Program of China grant No. 2023YFA1608000. Y.X.W. is a member of the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne.

Appendix A Benchmark of the chemical reaction network

Our reaction network from Wang et al. (2021) is benchmarked against the reaction network from Ruaud et al. (2016), which is provided as an example network in the Nautilus package. The initial abundances are already listed in Table 1. The TMC1 model proposed by Semenov et al. (2010) is utilized in the benchmark, which is assumed as a stable physical model with the temperature T = 10 K, the hydrogen nuclei density n_H = 2 × 10⁴ cm⁻³, the visual extinction A_V = 10 mag, the cosmic ray ionization rate ζ_CR = 1.8 × 10⁻¹⁶ s⁻¹, and the FUV radiation field χ = 1. Other chemical parameters are the same as those described in Sect. 2. Figure A.1 shows the time-dependent evolution of the fractional abundances of several important species with respect to H₂ calculated by the Nautilus code with adopting these two networks, respectively. The comparison shows good agreement between two simulations.

Fig. A.1 Time-dependent evolution of the fractional abundances of several important species (CO, CO2, H2O, CH3 OH, CH4, NH3, HC3N, N2H+, and HCO+) with respect to H2 using the TMC1 model calculated by the Nautilus code. The letter J represents those species on the grains. The black solid lines are the results adopting the reaction network from Wang et al. (2021), and the red crosses represent the results adopting the reaction network from Ruaud et al. (2016), as an example network in the Nautilus package.

Appendix B Additional simulated HC₃N/N₂H⁺

Fig. B.1 Same as Fig. 1, but adopting model 1.

Fig. B.2 Same as Fig. 1, but adopting model 2.

Fig. B.3 Same as Fig. 1, but adopting model 3.

Fig. B.4 Same as Fig. 1, but adopting model 5.

Appendix C Candidate chemical clocks

Fig. C.1 Time-dependent evolution of other average ratios within 105 AU from the center adopting model 6 as candidate chemical clocks.

Appendix D Other average ratios

Fig. D.1 Other average ratios as in Fig. C.1 but do not exhibit monotonic evolutionary trend throughout the entire evolution.

References

All Tables

All Figures

Fig. 1 Correlation between the average ratio of HC3N/N2H+ and the average abundance of HC3N with respect to H2 within four different ranges (panel A: 7 × 103 AU, panel B: 104 AU, panel C: 5 × 104 AU, and panel D: 105 AU) around the center by adopting model 4 in chemical simulations compared with observations. The results during the HMSC, HMPO, HMC, and UCHII stages are plotted as squares, circles, diamonds, and triangles, respectively. The evolution of HC3N/N2H+ and HC3N/H2 throughout ttotal during each stage consistently decrease (from right to left and from top to bottom). The observations proposed by Taniguchi et al. (2019b) and Wang et al. (2023, including derived column density of HC3N, N2H+, and H2 simultaneously) are plotted as solid squares, solid circles, and solid triangles, respectively, with error bars.

In the text

	Fig. 2 Same as Fig. 1, but adopting model 6.
In the text

Fig. 3 Time-dependent evolution of the average ratio of HC3N/N2H+ and the average fractional abundance of HC3N and N2H+ with respect to H2 within four different ranges (7 × 103 AU, 104 AU, 5 × 104 AU, and 105 AU) around the center using model 4 (left panels) and model 6 (right panels). The four evolutionary stages are labeled in each panel, and the observations are also plotted, with the dashed lines as the median values and the dotted lines as the maximums or minimums of observations. The evolution before 104 yr is excluded to show the changes in different stages more clearly. The results within 7 × 103 AU and 104 AU during the HMSC stage in model 4 completely overlap because of the large rin > 104 AU.

In the text

Fig. 4 Same as Fig. 3, but depicting the ratio of HC3N/N2H+ and the fractional abundance of HC3N and N2H+ with respect to H2 at six different radii (102 AU, 103 AU, 7 × 103 AU, 104 AU, 5 × 104 AU, and 105 AU) from the center using model 4 (left panels) and model 6 (right panels). The results at r = 102 AU, 103 AU, 7 × 103 AU, and 104 AU during the HMSC stage in model 4 completely overlap because of the large rin > 104 AU.

In the text

Fig. A.1 Time-dependent evolution of the fractional abundances of several important species (CO, CO2, H2O, CH3 OH, CH4, NH3, HC3N, N2H+, and HCO+) with respect to H2 using the TMC1 model calculated by the Nautilus code. The letter J represents those species on the grains. The black solid lines are the results adopting the reaction network from Wang et al. (2021), and the red crosses represent the results adopting the reaction network from Ruaud et al. (2016), as an example network in the Nautilus package.

In the text

	Fig. B.1 Same as Fig. 1, but adopting model 1.
In the text

	Fig. B.2 Same as Fig. 1, but adopting model 2.
In the text

	Fig. B.3 Same as Fig. 1, but adopting model 3.
In the text

	Fig. B.4 Same as Fig. 1, but adopting model 5.
In the text

	Fig. C.1 Time-dependent evolution of other average ratios within 105 AU from the center adopting model 6 as candidate chemical clocks.
In the text

	Fig. D.1 Other average ratios as in Fig. C.1 but do not exhibit monotonic evolutionary trend throughout the entire evolution.
In the text