Investigating the hot molecular core, G10.47+0.03: A pit of nitrogen-bearing complex organic molecules

Cesaroni et al. (2010) resolved the hot core G10 and found three distinct HII regions inside the core, which are A (Ultra Compact HII region), B1, and B2 (Hyper compact HII regions). Rolffs et al. (2011b) also observed this source with the submillimeter Array (SMA). They observed the continuum at three different frequencies (201/211 GHz, 345/355 GHz, 681/691 GHz). The beam sizes were insufficient to resolve the continuum emission at 201/211 and 681/691 GHz. However, at 345/355 GHz, the beam size was sufficient ($\sim 0.3^{{}^{\prime\prime}}$) to resolve this source. Here, continuum maps of G10 were observed at four different wavelengths, 1.88 mm (159.45 GHz), 1.94 mm (153.96 GHz), 2.02 mm (148.51 GHz), and 2.3 mm (130.50 GHz), respectively. Summaries of continuum images (synthesized beam, peak flux, integrated flux, etc.) are discussed in detail in Gorai et al. (2020). Gorai et al. (2020) calculated the H₂ column density for each continuum map and obtained an average value of $\sim$ 1.35$\times$10²⁵cm^-2. Gorai et al. (2020) also estimated the dust opacity of each continuum map and found all are optically thin ($\tau_{\nu}$ $\sim$0.135).

All lines are extracted from the 2.0${}^{{}^{\prime\prime}}$ (17200 AU) diameter region centered at G10 ($\alpha$(J2000)=18^h08^m38.232^s, $\delta$(J2000)=-19⁰51^′50.4${}^{{}^{\prime\prime}}$). Line identification has been carried out with CASSIS together with the spectroscopic database, ‘Cologne Database for Molecular Spectroscopy’ (CDMS, Müller et al., 2001, 2005)³³3(https://www.astro.uni-koeln.de/cdms) , and the Jet Propulsion Laboratory (JPL, Pickett et al., 1998)⁴⁴4(http://spec.jpl.nasa.gov/) database. The observed transitions of various nitrogen-bearing molecules together with their quantum numbers, upper state energies (E_u), V_LSR, line parameters such as peak brightness temperature (I_max), line width (full width at half maximum; FWHM), and the integrated intensity ($\rm{\int T_{mb}dV}$) are noted in Table 5. The line parameters (I_max, FWHM, $\rm{\int T_{mb}dV}$) are obtained using a single Gaussian fit to the observed spectral profile of each unblended transition. To see the spatial distribution of the observed species, moment zero maps of each transition are shown in Figures 21, 22, and 23, and 24. The emitting regions of all the transitions are calculated using a two-dimensional Gaussian fitting of their moment maps (see Table 2).

Among all the observed nitrogen-bearing species toward G10, we only find multiple unblended transitions (more than 2) for four molecules (see Table 5), which are vinyl cyanide (C₂H₃CN), ¹³C isotopologue of vinyl cyanide (¹³CH₂CHCN), ethyl cyanide (C₂H₅CN), and ¹³C isotopologue of ethyl cyanide (¹³CH₃CH₂CN). In addition to these four nitrogen-bearing molecules, multiple transitions are observed for methyl acetylene (CH₃CCH). Thus, a rotational diagram analysis was carried out for these five species to measure their excitation temperatures and the total column densities.

For the optically thin transitions, the column density (N${}_{u}^{thin}$) at local upper state can be expressed as

where g_u is the degeneracy of the upper state, k_B is the Boltzmann constant, $\rm{\int T_{mb}dV}$ is the integrated intensity, $\nu$ is the rest frequency, $\mu$ is the electric dipole moment, and S is the transition line strength, (Goldsmith & Langer, 1999). Under the local thermodynamic equilibrium (LTE) condition and for optically thin transitions, the total column density ($N_{total}$) of a molecular species can be prescribed as

where $T_{rot}$ is the rotational temperature, E_u is the upper state energy, $\rm{Q(T_{rot})}$ is the partition function at rotational temperature. This can be rearranged as

There is a linear relationship between the upper state energy and column density at the upper level when all lines are optically thin and all lines are probing the same physical conditions. The column density and rotational temperature are extracted from the rotational diagram. Rotational diagrams of all the five species are shown in Figures 1, 3, and 6.

Though the rotational diagram is used for the species for which multiple transitions are observed, another method, Markov chain Monte Carlo (MCMC), is used to constrain the excitation temperature and column density for the comparison. The MCMC algorithm was initially developed for dealing with the probabilistic behavior of the collection of atomic particles. However, it was challenging to do this analytically. Thus, an iterative process was implemented to solve it through simulation. This stochastic model describes a sequence of possible events in which the probability of each event depends only on the state attained in the previous one. The MCMC method is an interactive process that goes through all line parameters (e.g., molecular column density, excitation temperature, source size, line width) with a random walk and heads into the solution’s space. $\chi^{2}$ minimization gives the final solution.

For MCMC fitting, the FWHM of all the transitions of a species are kept constant at their average measured value, v_LSR is kept constant at $68$ km/s, and source size is kept constant at 2^". Then, the obtained results of these species are compared with those obtained with the rotation diagram method. Obtained results with these two methods are presented in Table 1 and discussed in Section 3. The obtained line parameters of the MCMC fitting of five species are noted in Table 6. MCMC fitted spectra are shown in Figures 18, 19, and 20.

Cyanopolyynes are a group of species with the chemical formula HC_nN $(n=3,5,7..)$. In addition, it contains the cyano group (CN). Previously, Rolffs et al. (2011b) identified HCN, HNC, H¹³CN, and HC₃N in G10. Here, the identifications of cyanoactylene (HC₃N) and its isotopologues and cyanodiacetylene (HC₅N) and its isotopologue are reported.

Six transitions of cyanamide (NH₂CN) are observed toward G10 (Table 5). However, two of them are clearly detected. Our simple LTE model yields a good fit (Figure 7) to the observed spectra with a column density of $4.8\times 10^{15}$ cm^-2. The emitting region (see Figure 23 and Table 2) of the two unblended transitions are compact ($\sim 0.91^{{}^{\prime\prime}}$).

Aminoacetonitrile (H₂NCH₂CN) is a well-known precursor of glycine. It hydrolyses to form glycine. Danger et al. (2011) studied the NH₂CH₂CN formation under the interstellar condition by thermal processing of ices containing methanimine, ammonia, and hydrogen cyanide. Five transitions of aminoacetonitrile are observed toward G10 ( 5), but four of them are blended with other species. Only the transition at 154.51748 GHz with upper state energy of 86 K is successfully identified. Therefore, we consider that this species was tentatively detected in this observation. A column density of $5.5\times 10^{15}$ cm^-2 yields a good fit to the observed spectra (see panel 16 of Figure 7). The emitting diameter of aminoacetonitrile is 1.07${}^{{}^{\prime\prime}}$. The moment zero map of this transition is shown in panel 8 of Figure 23.

Snyder & Buhl (1973) first identified methyl acetylene (CH₃CCH) in Sgr B2. Though methyl acetylene is not a nitrogen-bearing species, being a symmetric top molecule, its identification is reported here because of its utility as a tracer of the physical condition of the star-forming region. Three transitions of methyl acetylene are clearly identified. All the transitions correspond to $J=9-8$. The three transitions have $k_{a}$= 6, 3, and 2, respectively. This is the first identification of CH₃CCH in G10. All three unblended transitions are found to be optically thin. Rotational diagram analysis of the CH₃CCH (Figure 6) depicts a column density of $6.0\times 10^{16}$ cm^-2 and a rotational temperature of $189$ K. Panels 1-3 of Figure 24 shows moment map of CH₃CCH for different $k_{a}$ components of $J=9-8$ transitions. Interestingly, it shows that the emitting regions of CH₃CCH transitions are related to their corresponding $k_{a}$ components. For example, the emitting region of $J=9-8$ transitions with $k_{a}$ = 2 is 1.8${}^{{}^{\prime\prime}}$ ($E_{u}=65.8K$), and for $k_{a}=3$ it is $1.9^{{}^{\prime\prime}}$ ($E_{u}=101.7$ K). Since with our spatial resolution these two transitions are marginally resolved, we did not find significant differences for close-lying energy states. However, for $k_{a}=6$ ($E_{u}=296.9$ K), a compact emission ($1.28^{{}^{\prime\prime}}$) is obtained. With the increase in the $k_{a}$ component, the upper state energy increases, which yields a more compact emitting region. All three detected transitions have $\Delta k_{a}=0$. The ratio of lines between the same J-state but other $k_{a}$-states can be used as temperature probe (Mangum & Wootten, 1993). To estimate the kinetic temperature ($T_{K}$), we chose some selected line ratios ($\frac{J1_{ka}-J2_{ka}}{J3_{ka^{\prime}}-J4_{ka^{\prime}}}$). These ratios are selected based on the following criteria (Mangum & Wootten, 1993; Das et al., 2019): 1) $\Delta J=1$, J1=J3, and J2=J4, 2) $\Delta k_{a}=0$ for both the transitions, 3) ${k_{a}^{\prime}>k_{a}}$, 4) both the transitions are optically thin, 5) The ratio between the difference in excitation temperature belongs to different k_a values.

With the LTE approximation, the line ratio (R) between the two transitions satisfying the above criteria can be expressed as (Mangum & Wootten, 1993)

where $D=E(J3,k_{a}^{\prime})-E(J1,k_{a})$ and $S_{R}=\frac{S_{J1k_{a}}}{S_{J3k_{a}^{\prime}}}$. Mangum & Wootten (1993) derived this relation for H₂CO, which is also applicable to any slightly asymmetric or symmetric rotor species.

According to the above selection rules, in between the three observed unblended transitions, $\frac{9_{2}-8_{2}}{9_{3}-8_{3}}$, $\frac{9_{2}-8_{2}}{9_{6}-8_{6}}$, and $\frac{9_{3}-8_{3}}{9_{6}-8_{6}}$ are considered for estimating the excitation temperature. The black solid curve in Figure 8(a-c) depicts the expected line ratio with the excitation temperature. The observed line ratio of these three transitions is shown with the solid red line. Interestingly, the observed line ratio intersects the curve between 133-193 K, closer to the estimated hot core temperature.

Additionally, two transitions ($J=9-8$ with $k_{a}=0,3$) of ¹³C isotopologue of methyl acetylene (CH${{}_{3}}^{13}$CCH) are observed toward G10. Our LTE synthetic spectra fit the observed spectra well, with a column density of $8.0\times 10^{15}$ cm^-2. A similar trend of emitting regions is also obtained here. For $k_{a}=0,the$ emitting region is 1.31${}^{{}^{\prime\prime}}$, whereas it is $1.01^{{}^{\prime\prime}}$ for $k_{a}=3$. With the simple LTE calculation, a carbon fractionation ($\frac{{}^{12}C}{{}^{13}C}$ ratio) of $8.8$ is obtained between CH₃CCH and CH${{}_{3}}^{13}$CCH. By considering the column density obtained with the rotational diagram and MCMC methods (for $\rm{CH_{3}CCH}$), a $\frac{{}^{12}C}{{}^{13}C}$ ratio of 7.5-10.9 is obtained.

A similar process for estimating the kinetic temperature from the transitions obtained with the CH${{}_{3}}^{13}$CCH is carried out. Following the selection rule, $\frac{9_{0}-8_{0}}{9_{3}-8_{3}}$ is used to estimate the kinetic temperature. Figure 8d depicts that the observed line ratio intersects the expected line ratio at $102$ K.

We identify several transitions of nitrogen-bearing species along with their various ¹³C isotopologues. Only the unblended transitions are used for further analysis. In the case of the ¹³C isotopologue of vinyl cyanide, we identify four transitions of ¹³CH₂CHCN, one of $\rm{C{H_{2}}^{13}CHCN,}$ and one of CH₂CH¹³CN. We obtain two transitions from one of the ¹³C isotopologues of ethyl cyanide ($\rm{CH_{3}C{H_{2}}^{13}CN}$). Two transitions of the ¹³C isotopologue of cyanoacetylene (HC¹³CCN and $\rm{HCC_{13}CN}$) and one transition of the ¹³C isotopologue of cyanodiacetylene (HCC¹³CCCN) are identified. Furthermore, we identify two transitions of a ¹³C isotopologue of methyl acetylene ($\rm{C{H_{3}}^{13}CCH}$). We use the rotational diagram analysis and MCMC method to estimate the column density of these species if more than two transitions were identified. Otherwise, we consider the simple LTE approximation. Despite all these issues, we found that the obtained $\frac{12_{C}}{13_{C}}$ ratio of all these species varies from 0.8-72. In the case of methyl acetylene, vinyl, and ethyl cyanide, it is 7.5-72, which is generally consistent with the elemental $\frac{12_{C}}{13_{C}}$ ratio ($\sim 89$). For the cyanopolyynes (cyanoacetylene and cyanodiacetylene), we obtain a small ¹³C enrichment (0.77-2.93). Since not many transitions were obtained for the cyanopolyynes, we restrict our comments based on the upper limit of the column densities obtained.

We started with the simple 0D simulation. We simulate the source as a single point and do not consider any spatial distribution of physical parameters ($n_{\rm H}$, $A_{v}$, $T_{g}$, $T_{d}$). Following Gorai et al. (2020), we carried out our first simulation with generally accepted parameters for a source such as G10 (initially $n_{\rm H}=10^{3}$ cm^-3 , which goes up to $10^{7}$ cm^-3 during cloud collapse; initial $A_{v}=2,$ which goes up to 250; initial $T_{g}=40$ K, which goes up to 200 during the warm-up stage; and initial $T_{d}$ = 20 K, which goes up to 200 K). Since G10 is massive hot core, we selected a collapse time of $10^{5}$ years and a warm-up time of $5\times 10^{4}$ years.

Knowing that chemistry is very sensitive to certain physical parameters such as temperature and density, we are required to test a number of 0D models by varying the values of most relevant physical parameters. Hence, first we repeat our simulation with a different peak core density of $10^{6}$ cm^-3. We also simulated a number of different models where we used different values of peak gas temperature (400, 300, and 200 K) during the warm-up phase and initial dust temperature (15, 20, and 25 K) during the collapse phase to cover the possible physical evolution history of hot cores such as G10. The warm-up phase is a very important phase in chemical evolution history as the dust grain chemistry is very sensitive to dust temperature. This makes the warm-up duration a very critical parameter. Therefore, we used three different values for the total duration ($T_{w}=10^{4},5\times 10^{4},$ and $10^{5}$ years) of the warm-up phase in order to properly investigate the dependence of simulated abundances of N-bearing species on the warm-up time of the collapsed cloud.

Depending on the values of physical parameters used during the time evolution of these three-stages, a number of 0D models are implemented to study the chemical evolution in G10-like dense cores. In Table 4, we list values of different physical parameters used in different models. More details about these models can be found in Section 5.1, where we discuss results for the 0D model.

To study the spatial distribution of the chemical abundances of G10, we considered the the density and temperature distribution obtained by SED fitting by van der Tak et al. (2013). Very recently, Bhat et al. (2021) constructed a similar model for G31.041+0.31. They divided the cloud envelope into 23 logarithmic spacing shells. Here, for G10, we divided the cloud envelope into ten logarithmic spacing shells. The time evolutions of the physical parameters ($n_{\rm H}$, $A_{v}$, $T_{d}=T_{g}$) for these ten shells are shown in Figure 9.

For each shell, we consider that, at $t=0$ year, all the shells have $n_{H}=10^{3}$ cm^-3. In the first stage, every shell is allowed to evolve up to its spatial density distribution obtained by van der Tak et al. (2013). For example, for the fifth shell’s (at $\sim 1717$ AU) density may evolve up to $2.08\times 10^{7}$ cm^-3 at the end of $10^{5}$ years. Beyond that (warm-up and post-warm-up stages), the density of each shell is kept constant at its maximum value. Since dust temperature is a critical parameter, we used two profiles (15 and 20 K) for the initial dust temperature. For simplicity, we kept the gas temperature $T_{g}$ equal to the dust temperature $T_{d}$. In the first stage, $T_{d}=T_{g}$ was kept constant at initial value (15 or 20 K depending on model). In the warm-up stage, $T_{d}=T_{g}$ evolves from its initial value to final value respective of the spatial temperature distribution as reported by van der Tak et al. (2013) (for example, $207$ K for the fifth shell). In the post-warm-up stage, gas and dust temperatures are kept constant at their highest values respective to each shell. To calculate $A_{V}$, we used the relation $A_{V}=A_{V,0}\times(n_{\rm{}H}/n_{\rm{}H,0})^{(2/3)}{}$, as given in Sect. 3.3 of Garrod & Pauly (2011), where A_V,0 is the initial value of visual extinction and $n_{\rm{}H,0}$ is the initial total number density. We used A${}_{V,0}=2$ and $n_{\rm{}H,0}=10^{3}$ cm^-3.

To explain the observed abundances of chemical species in G10, a three-phase (gas phase, grain surface, and grain mantle) version of the Nautilus gas-grain chemical code (Ruaud et al., 2016) is used. This code is based on the rate equation approximation (Hasegawa et al., 1992; Hasegawa & Herbst, 1993). In all our simulations, the surface of the dust grain is defined as the top two monolayers of species. All ice species under the surface are treated as mantles. When a species desorbs from the surface, another species from the mantle comes to the top to form the new surface layer, and, similarly, with the accretion of a new species on the surface, another surface species moves to the mantle to form a new mantle layer. In our model, standard physical processes such as accretion, diffusion, recombination, and thermal desorption of species are considered. The cosmic-ray-induced desorption, UV (direct and indirect) photo-desorption, and chemical desorption are also considered. These physical processes are explained in detail in Ruaud et al. (2016); Iqbal et al. (2018); Das et al. (2019); Wakelam et al. (2020). We used a common cosmic-ray ionization rate of $1.3\times 10^{-17}$ s^-1. In Table 3, we list the initial abundances of all the gas and the grain species we used in our all chemical models.

Our gas-phase chemistry is mostly based on the kida.uva.2014 public network (Wakelam et al., 2015). We updated our gas phase chemistry by adding all relevant reactions from KIDA and UMIST databases. Apart from this we also followed Table 8 and Table 9 of Gorai et al. (2020), and Table 5 of Sil et al. (2018) to complete our gas-phase nitrogen chemistry. Our surface chemistry is mostly based on Iqbal et al. (2018). To complete the surface chemistry, we added all relevant formation and destruction reactions in our chemical network from Table 14 of Belloche et al. (2009) and Table 8 of Gorai et al. (2020). We also added all surface reactions, related to the formation and destruction of C₂H₃CN and C₂H₅CN, given in Table B.3 and Table B.4 of Garrod et al. (2017).

Density is one of the fundamental properties of a core that determines the future state of the core. It is challenging to consider the final density of a collapsed core just after its isothermal collapse without having comprehensive knowledge of its collapse. Hence, in our simulations, we adopt two final density profiles for G10 to study their impact on chemical abundances at the end of the simulation.

Here, we show the results for 1D simulations, which are essential to studying the possible spatial distribution of N-bearing species in the source G10. We ran two simulations (based on initial dust temperatures of 15 and 20 K) of the 1D model. In each simulation, we have ten grids representing ten spatial points in the source. Details on these spatial points are given in Section 4.2.

The results showing the time evolution of the abundances of selected species for all ten grids are shown in Figure 15. In this simulation, the initial gas and dust temperature was set at 15 K. The same plot for an initial dust and gas temperature of 20 K is given in Appendix E, Figure 29. The observed abundances are marked with horizontal bars. Grid 1 is the innermost grid, whereas Grid 10 is the outermost. Our model shows that the abundances of all species vary by a few orders of magnitude between the inner and outer grids. Some species such as HC₃N and HC₅N are most abundant at inner grids, while others such as NH₂CN and CH₃CCH are most abundant in the outer grids.

For a better comparison with values presented in Rolffs et al. (2011b), we show the spatial distribution of peak abundances of selected species in Figure 16. Again, this plot is for models with initial dust temperatures of 15 K. A similar plot with an initial dust temperature of 20 K is given in Appendix E, Figure 30. The dotted and dashed green vertical lines represent the observed spatial resolutions of Rolffs et al. (2011b), whereas the dashed pink vertical line represents the spatial resolution of our observation. Peak abundances (black curves) are extracted by noting the abundances achieved during warm-up and post-warm-up periods. To avoid the time uncertainty of the peak values, we also show the abundances obtained just after the warm-up phase ($1.5\times 10^{5}$ years) and final time ($2.5\times 10^{5}$ years).

As shown in Figure 16, the spatial resolution of the observation presented in Rolffs et al. (2011b) is much higher than that of the observations in this paper. We also note that the simulated peak abundances of all species is at its maximum toward the inner grids, though not in the same place. This is important to understand why, for most of the species, Rolffs et al. (2011b) recorded a higher maximum abundance compared to us for the source G10. Our ID simulations shows a high abundance of HCN and HC₃N deeper inside the cloud (i.e., around the high-temperature region), which could not be resolved with observation presented in this paper. However, HC₃N is more extended compared to HC₅N. This is reflected in our modeling results (Figure 16 and Figure 30). They show that HC₅N decays in the outer part of the envelope more rapidly compared to HC₃N. The abundance of vinyl cyanide and ethyl cyanide seems to be heavily dependent on the initial ice temperature. In the model with the initial gas and dust temperature of 20 K, shown in Figure 29, abundances of vinyl cyanide and ethyl cyanide drop rapidly compared to the model with the initial gas and dust temperature of 15 K, shown in Figure 15. More interestingly, inside 3000 AU, it depicts a higher abundance of ethyl cyanide than vinyl cyanide. Rolffs et al. (2011b) also obtained a higher abundance of ethyl cyanide ($9\times 10^{17}$ cm^-2) compared to vinyl cyanide ($7\times 10^{17}$ cm^-2). In contrast, we obtained a higher abundance of vinyl cyanide than ethyl cyanide with our resolution. This is also reflected in our model. So, our model at the initial ice phase temperature of 15 K successfully explains a high abundance of ethyl cyanide over vinyl cyanide deep inside, and of vinyl cyanide over ethyl cyanide at the outer part of the envelope. An abundance of CH₃NH₂ seems roughly invariant across the cloud. The abundances of CH₂NH and CH₃NH₂ are overproduced in our 15 K model, whereas they are underproduced in the 20 K model. Results presented here can explain the observed abundances of the peptide-like bond containing species (HNCO, CH₃NCO, and NH₂CHO). With the initial gas and dust temperature of 15 K, CH₃CCH shows a good match with the observational results (see Figure 16).

Simulated (obtained in best-fit 1D model) and observed column densities of the identified species are shown in Figure 17. Various methods are implemented to derive the observed column density of the species. RD (shown in green) represents the column density obtained with the rotational diagram method, and MCMC (shown in blue) illustrates the MCMC Method. LTE1 (shown in red) represents the cases when a fixed source size of $\sim 2^{"}$ is used, whereas LTE2 represents the circumstances for which the average emitting region of species are used (shown in black). Figure 17 shows that all the methods of deriving the observed column density of a species agree well among themselves. The simulated column densities of species are shown in yellow. The grid number (noted in black) represents the peak column density of that grid, for which we have a good match with the observed results. Additionally, for the cases of HC₃N and HC₅N, we include the results obtained for the same grid at $\sim 2.5\times 10^{5}$ years.