Computational insights into the chemical reaction networks of C₃H₆O₃, C₃H₇O₃ and C₂H₅O₂: implications for the interstellar medium

Abstract

The formation of complex organic molecules (COMs) in the interstellar medium (ISM) is central to astrochemistry and prebiotic chemistry, as these species may act as precursors to biomolecules essential for life. Among COMs, glyceraldehyde (HOCH₂CH(OH)C(O)H, GCA) has attracted attention as a potential building block in early biochemical pathways. Although GCA has not yet been detected in the ISM, the presence of structurally related compounds in various astronomical environments suggests that it may form under interstellar conditions. In this study, we employed the automated reaction discovery tool AutoMeKin to systematically explore the gas-phase chemical reaction networks (CRNs) of C₃H₆O₃ (GCA), C₃H₇O₃ (a hydrogenated analog), and C₂H₅O₂. Reaction pathways were characterized at the ωB97XD/Def2-TZVPP level of theory, and rate coefficients for key processes were computed using the competitive canonical unified statistical (CCUS) model, which accounts for multiple dynamic bottlenecks. Our analysis revealed several barrierless pathways leading to GCA or to GCA and a leaving group. Notably, the reaction between glyoxal (HCOHCO) and the HOCHCH₂OH radical, though neither has yet been detected in the ISM, was found to efficiently produce GCA and a formyl radical, with rate coefficients on the order of 5.4–7.9 × 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹ across the 10–100 K temperature range. However, aside from the aforementioned exception, most GCA formation channels result in highly vibrationally excited intermediates that are more likely to undergo rapid unimolecular decomposition than to be stabilized by radiative emission under typical ISM conditions. These results suggest that while gas-phase GCA formation is chemically feasible, it is likely transient and difficult to detect directly. In contrast, alternative products such as formaldehyde, glycolaldehyde, and (Z)-ethene-1,2-diol dominate many pathways and align better with current astronomical observations. This work provides detailed mechanistic and kinetic insights that enhance astrochemical modeling and advance our understanding of molecular complexity in star-forming environments. Furthermore, it highlights the utility of automated CRN exploration for uncovering viable synthetic routes to prebiotic molecules in space.

---

1 Introduction

The formation of complex organic molecules (COMs) in the interstellar medium (ISM) is a key area of research in astrochemistry, as these molecules are considered precursors to life. Among COMs, glyceraldehyde (HOCH₂CH(OH)C(O)H, GCA), the simplest sugar, has garnered attention due to its potential role in prebiotic chemistry and its relevance as a precursor to biologically relevant molecules.^1,2

The gas-phase decomposition of GCA has been investigated experimentally through pyrolysis,² and the early steps of the Formose mechanism have been studied computationally using density functional theory (DFT).¹ Under pyrolytic conditions, GCA starts decomposing at 400 °C, undergoing retro-aldol fragmentation and dehydration. The retro-aldol reaction predominantly yields glycolaldehyde and formaldehyde.^3,4

Although GCA has not yet been detected in the ISM despite dedicated observational searches,⁵ its structural analogs, such as glycolaldehyde and ethylene glycol, have been observed in star-forming regions and comets,^6–9 suggesting that it may also form under astrochemical conditions. Laboratory experiments and computational studies indicate that GCA may originate from formaldehyde-based reactions in interstellar ice analogues, driven by ultraviolet irradiation, radical recombination, or thermal processes,² or from the reaction of (Z)-ethene-1,2-ediol with hydroxycarbene (HOCH).¹⁰ Understanding the formation pathways of GCA and related species in the ISM provides valuable insights into the chemical processes that drive molecular complexity and, ultimately, the emergence of life's building blocks in space.

This study explores potential gas-phase formation pathways for GCA, as well as other key reactions relevant to astrophysical environments, evaluating their viability under interstellar conditions. The formation of COMs in the ISM can proceed through various reaction mechanisms. Given the extremely low temperatures of interstellar clouds (10–150 K), gas-phase reactions must generally proceed without significant activation barriers.¹¹ In other words, such reactions are expected to follow energetically favorable pathways with submerged barriers,^12,13 often involving at least one radical or ion species. Additionally, due to the extremely low densities of interstellar clouds, product stabilization is inefficient unless radiative cooling provides a competitive mechanism.¹¹

In particular, we investigate two types of GCA gas-phase formation reactions:

R1˙ + R2˙ → GCA + hν (1)

R3˙ + M → GCA + H˙ (2)

Reaction (1), known as radiative association, involves two radicals (R1˙ and R2˙) and requires efficient radiative stabilization of GCA. In contrast, reaction (2) is a radical–molecule reaction that does not rely on radiative stabilization, as it produces GCA and atomic hydrogen as separate products.¹¹ Additionally, we also examined the possibility of reactions analogous to (2), where the leaving group is a small radical rather than H˙.

To explore these reaction mechanisms, we employ AutoMeKin,^14–16 an automated tool for discovering chemical reaction networks (CRNs). AutoMeKin has been successfully used to propose gas-phase formation routes of various interstellar molecules, including benzene,¹⁷ cyanoketene,¹⁸ cyanamide,¹⁹ formaldehyde,²⁰ glycolonitrile,²¹ hydrogen isocyanide,²² indole,¹³ malononitrile,²³ vinyl alcohol,²⁴ and (Z)-1,2-ethenediol.²⁵

In this study, we employ AutoMeKin to construct detailed reaction pathways for species with the chemical formulas C₃H₆O₃ and C₃H₇O₃, corresponding to reactions (1) and (2), respectively. C₃H₆O₃ represents GCA, while C₃H₇O₃ corresponds to GCA with an additional hydrogen (GCA–H). Additionally, we investigate the C₂H₅O₂ CRN using AutoMeKin due to its close relationship with the other CRNs. Through the analysis of these networks, we identify key intermediates and reaction channels, and we perform kinetic calculations to determine rate coefficients for the most relevant reactions. The computed rates provide valuable inputs for astrochemical models.

Our findings have the potential to bridge gaps in current astrochemical models by providing accurate rate coefficients and detailed reaction mechanisms. As observational techniques continue to advance, the insights gained from this work will contribute to a deeper understanding of the chemical evolution of star-forming regions and the origins of molecular complexity in space.

2 Computational details

2.1 Automated mapping of reaction mechanisms

The CRNs of C₃H₆O₃, C₃H₇O₃ and C₂H₅O₂ were generated using AutoMeKin,^14–16 an automated software developed in our group. This open-source program enables an efficient exploration of molecular decomposition pathways by integrating reactive MD simulations, graph-theoretic methods, and interactive visualization tools. The workflow begins with a preliminary potential energy surface (PES) exploration using semi-empirical calculations (hereafter referred to as level 1), followed by refinement of stationary points with a higher-level DFT or ab initio method (hereafter level 2). Consistent with previous studies, level 1 corresponds to PM7,²⁶ while level 2 employs ωB97XD/Def2-TZVPP. The latter has been demonstrated to accurately estimate barrier heights.²⁷

Additionally, AutoMeKin now also incorporates an improved algorithm for identifying barrierless reaction pathways.¹³ Since this algorithm relies on level 1 calculations, we further supported the barrierless character of selected pathways by performing relaxed scan optimizations at level 2, corrected with single-point B2PLYP/Def2-TZVPP calculations. As shown in the Supplementary Information (SI), the energy profiles of the selected pathways confirm the absence of barriers.

For a subset of key reactions, energies were further refined using CCSD(T)/CBS calculations (level 3). The CBS limit was estimated via extrapolation using the jun-cc-pVXZ basis set (with X = D, T).²⁸ level 2 and level 3 calculations were performed using Gaussian09.²⁹

Beyond its core capabilities in MD simulations and graph-theoretic analysis, AutoMeKin has recently been extended with a Python-based library, amk_tools, designed for parsing, processing, and analyzing CRNs.¹³ Custom Python scripts further allow for the adaptation of energy profiles generated by amk_tools and facilitate the comparison of the predicted decomposition products with known interstellar molecules.

For the C₃H₆O₃ and C₃H₇O₃ molecular systems, the AutoMeKin workflow was executed iteratively 30 times, with each cycle comprising MD simulations, optimization of minima and transition states (TSs), and CRN generation. Each MD simulation set included 500 trajectories, with a maximum simulation time of 0.5 ps. For the smaller C₂H₅O₂ system, the workflow was run for 20 iterations with 480 trajectories per cycle.

The transition-state filtering criteria from our previous study on indole decomposition¹³ were also applied here. These criteria served to remove structures with very low imaginary frequencies and to eliminate redundant TSs. Additional details on the methodological approach are provided in ref. 13.

To evaluate the convergence of the CRNs obtained through iterative refinement, kinetic Monte Carlo (KMC) simulations were conducted to solve the chemical master equation.³⁰ Specifically, the abundances of decomposition products were determined by exciting 10³ molecules of GCA and GCA–H at 250 kcal mol⁻¹. The reaction rate coefficients were estimated using RRKM theory,³¹ given by:

where, σ represents the reaction coordinate degeneracy, W^‡(E − E₀) is the sum of states of the TS, E₀ is the barrier height, and ρ(E) is the density of states of the reactant. These calculations were performed using AutoMeKin's kinetics module, applying the Beyer–Swinehart algorithm³² for direct counting of sums and densities of states. The product abundances obtained across iterative steps are shown in Fig. S1, indicating convergence after 10–15 iterations. Notably, water, glycolaldehyde, formaldehyde, and their isomers are among the predominant decomposition products of GCA, consistent with previous pyrolysis studies.^3,4

While our study primarily focuses on the barrierless formation pathways of GCA, the complete CRNs are available in Zenodo.³³ The analysis presented here is based on lumped networks, in which conformers are grouped together; however, the full networks, including all individual conformers, are also provided in Zenodo.³³ In these networks, minima, transition states, and products are labeled as MIN, TS, and PROD (or PR), respectively, with numerical identifiers distinguishing different species. The labels for each network are assigned independently. Molecular visualizations and representations were performed using CYLview.³⁴

2.2 Kinetic simulations

As in our previous study on glycolonitrile,²¹ we apply the competitive canonical unified statistical (CCUS)^35,36 model to reactions that exhibit two dynamic bottlenecks. Specifically, the first bottleneck usually corresponds to a free energy barrier for the association of two reactants, A and B, forming a vibrationally excited species (S*):The second bottleneck arises from the competitive formation of the products P₁ and P₂ from S*. These bottlenecks are associated with second-order saddle points: Note that in the above equations the rate coefficients are defined with respect to reactants A and B, rather than S*.The CCUS model enables an approximate of the rate coefficients for the formation of products P₁ and P₂ as follows:The association rate coefficient k_f between two reactants with permanent dipole moments can be approximated by:³⁷

k_f = Cμ^−1/2(d_Ad_B)^2/3T^−1/6 (7)

where μ is the reduced mass of the two reactants (in amu), T represents the temperature (in K), and d_A and d_B are the dipole moments (in D) of the reactants, calculated via level 2. The constant C has a value of 1.85 × 10⁻⁹, with units such that k_f is expressed in cm³ molecule⁻¹ s⁻¹. If one of the reactants does not have a permanent dipole moment, the association rate coefficient is instead given by C′μ^1/2(d_AQ_B)^1/2, where Q_B is the quadrupole moment (in D Å) and C′ = 3.52 × 10⁻¹⁰, with units ensuring that the rate coefficient is expressed in cm³ molecule⁻¹ s⁻¹.

The rate coefficients k₁ and k₂ are obtained using canonical variational transition state theory (CVT),³⁸ with A and B as the reacting species. For processes involving hydrogen transfer, tunneling corrections are included via the small curvature tunneling (SCT)³⁹ model. The rate coefficients were calculated under the assumption of zero pressure.

When there is no competition between pathways and S* can undergo only a single unimolecular reaction (k₁), the applied model is referred to as the canonical unified statistical (CUS) model, and the rate coefficient is given by k^CUS = k_fk₁/(k_f + k₁).

The vibrationally excited species S* may also undergo a radiative stabilization process:

The rate coefficient for this radiative stabilization was computed using the harmonic approximation:⁴⁰Here, N_m represents the total number of vibrational modes, P_i(n) is the probability of mode i being in level n, I_i is the infrared absorption intensity for the transition from level n = 0 to n = 1 of mode i (in units of km mol⁻¹), and ν_i denotes the frequency of the i-th vibrational mode (in cm⁻¹). The population distribution can be determined using the following equation:⁴⁰In this equation, ρ^N−1,i_vib refers to the density of states for the molecule with the i-th mode absent.

In the calculations above, all vibrational frequencies were scaled by the recommended factor of 0.975.⁴¹

3 Results

3.1 Analysis of the C₃H₆O₃ CRN

The screening of the C₃H₆O₃ CRN to find plausible GCA formation pathways consisted of seeking reactions that follow energetically favorable pathways with submerged barriers. This filtering process revealed several viable pathways for the gas-phase formation of glyceraldehyde. Fig. 1 illustrates the energy profile for the most probable and chemically relevant pathways, described from left to right:

Fig. 1 DFT-computed energy profile for the most probable formation routes of GCA via radiative association. The energies include scaled harmonic zero-point vibrational energy (ZPE), with the zero of energy set for GCA, shown at the bottom right of the figure. The diagram is based on the analysis of the lumped networks.

1. Purple and violet-blue pathways (left of Fig. 1). These routes involve hydroxycarbene (HOCH) reacting with either (Z)-ethene-1,2-diol (HOCHCHOH, violet-blue) or (E)-ethene-1,2-diol (purple). Among these three species, only (Z)-1,2-ethenediol has been detected in the ISM.¹⁰ Notably, the violet-blue route corresponds to the reaction mechanism hypothesized by Rivilla et al.¹⁰ These are the only pathways involving a reaction between two singlet species. Both routes lead to GCA in five elementary steps, with the resulting GCA* possessing an excess energy of approximately 80 kcal mol⁻¹. The first step in both routes involves the association of reactants to form stable van der Waals (vdW) complexes, which then rearrange to produce cyclopropane-1,2,3-triol (at 21.6 kcal mol⁻¹). The rearrangement in the violet-blue route is more favorable than in the purple route.

2. Dark blue pathway (center of Fig. 1). This route begins with the barrierless association of two radical species at 104.4 kcal mol⁻¹: hydroxymethyl (CH₂OH) and HOCCHOH radicals. This pathway converges with the purple one at the (Z)-prop-1-ene-1,2,3-triol intermediate (at 6.5 kcal mol⁻¹), leading to GCA* in just two steps, but with a significant excess energy of 104 kcal mol⁻¹. Additionally, while none of the reactants have been detected in the ISM, protonated formaldehyde (CH₂OH⁺) has been observed.⁴²

From the (Z)-prop-1-ene-1,2,3-triol intermediate, the convergence point of the three aforementioned pathways, the most favorable route proceeds through TS98 to yield 2-hydroxyprop-2-enal (CH₂COHCHO) and water, rather than forming GCA.

3. Light green pathway (center of Fig. 1). This route shares one reactant with the dark blue pathway: hydroxymethyl radical. The other reactant, the hydroxymethylcarbonyl radical (HOCH₂CO), has not been detected in the ISM, but it was identified as a product of the reaction between atomic hydrogen and glycolaldehyde in solid para-hydrogen at 3.2 K.⁴³ As in the previous case, this route leads to GCA in two steps.

4. Red route (right of Fig. 1). In this pathway, the formyl radical (HCO) undergoes a barrierless radical–radical association with HOCHCH₂OH, leading to GCA in a single step. Among the two reagents, only the formyl radical has been detected in the ISM.⁴⁴

5. Light blue route (right of Fig. 1). This pathway involves two reactants: the hydroxymethyl radical and the 2-oxoethylideneoxidanium radical (HOCHCHO), the latter of which has been proposed as a precursor to GCA formation in dark regions of the ISM.⁴³ Similar to the previous route, GCA forms through the association of the reactants in a single step.

Although some of the reactions described above appear plausible due to their barrierless nature, GCA* forms with a significant excess energy (62.9 kcal mol⁻¹ in the most favorable case, specifically in the light blue route). In this scenario, GCA* can only dissipate energy through radiative emission of IR light or undergo further unimolecular reactions. An RRKM calculation for the unimolecular decomposition of GCA at 62.9 kcal mol⁻¹ predicts 100% abundance of formaldehyde and ethene-1,2-diol, after surmounting an energy barrier of 36.8 kcal mol⁻¹ (see Fig. 2). The estimated radiative emission rate is 111 s⁻¹, which is orders of magnitude lower than the unimolecular decomposition rate of GCA*, leading to formaldehyde and 1,2-ethenediol, which is approximately 3 × 10⁶ s⁻¹. Given the large disparity between the rates of radiative emission and unimolecular decomposition, we have chosen to neglect radiative stabilization in all kinetic calculations, as the excess energy above the decomposition barrier is consistently high and this process is not expected to compete effectively.

Fig. 2 DFT-computed energy profile for the reactions CH₂OH + HOCHCHO → CH₂O + HOCHCHOH and HCO + HOCHCH₂OH → CH₂O + HOCHCHOH. The energies include scaled harmonic zero-point vibrational energy (ZPE), with the zero of energy set at GCA. The diagram is based on the analysis of the lumped networks.

The following reaction represents the outcome of the radical–radical association described in the light blue route:

CH₂OH + HOCHCHO → CH₂O + HOCHCHOH (11)

As shown in Fig. 2, there are two TSs with relatively similar energies, one leading to the Z isomer and the other to the E isomer, with the former being 2.3 kcal mol⁻¹ lower in energy.

The rate coefficients for reaction (11) evaluated using the CCUS model and presented in Table S1 of the SI, indicate that the abundance of the Z isomer is orders of magnitude higher than that of the E isomer. Notably, only the Z isomer has been detected in the ISM.¹⁰ However, since the reactants involved in reaction (11) have not yet been observed in the interstellar medium (ISM), this reaction cannot be considered a viable formation pathway for 1,2-ethenediol in the ISM.

The rate coefficients for the formation of formaldehyde and (Z)-1,2-ethenediol from reaction (11) in the 10–100 K temperature range vary from 8.7 × 10⁻¹⁰ and 6.0 × 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹.

Fig. 2, along with our kinetic calculations, shows that the predominant decomposition pathway of GCA leads to the formation of formaldehyde and 1,2-ethenediol. Therefore, the light green, red and light blue pathways depicted in Fig. 1 exclusively yield these two products, with GCA serving as a key intermediate. In contrast, the purple, violet-blue, and dark blue pathways result in the formation of water and 2-hydroxyprop-2-enal.

3.2 Analysis of the C₃H₇O₃ CRN

The starting geometry for our analysis with AutoMeKin was a radical formed by GCA with an extra hydrogen atom (GCA–H). The rationale for studying this CRN was to propose that GCA could form according to reaction (2). In this way, the unimolecular decomposition of the system could lead to the formation of GCA, releasing a radical hydrogen, as seen in reaction (2).

As in the CRN described in the previous section, the screening process was similar, i.e., reactions where atomic hydrogen is released and that involved barriers whose energy were below the energy of reactants were sought. Additionally, among the possible reactants, those that could rearrange through an exothermic hydrogen transfer were discarded, as this process would potentially be much faster than reaction (2).

We found that the reactions between the reactant pairs shown in Fig. 3 are particularly relevant.

Fig. 3 Potential reactants that undergo a barrierless association leading to GCA–H. Each is labeled using AutoMeKin's notation (PR followed by a number) to facilitate the identification of the structures within the C₃H₇O₃ CRN. The diagram is based on the analysis of the lumped networks.

The reaction between the fragments depicted in Fig. 3a, although leading to the formation of GCA–H without a barrier in a single step, does not yield GCA and H, as their energy is higher than that of the reactants. However, this reaction is of particular interest because no lower-energy pathway exists for these reactants other than its dissociation back into the original fragments. Consequently, the only possible outcomes of a reactive encounter between these two species are either immediate dissociation or stabilization of the complex via spontaneous radiative emission. As mentioned above, formaldehyde is a compound detected in the ISM,⁴⁵ whereas the larger fragment, the HOCHCH₂OH radical, which also appears in the red route of Fig. 1, has not been observed. However, the reaction between hydroxyacetaldehyde, which has been detected,⁴⁶ with atomic hydrogen to form HOCHCH₂OH has been reported to be exothermic by 24.4 kcal mol⁻¹, albeit with a significant barrier.⁴³

The reactants depicted in Fig. 3b would lead to the formation of GCA + H˙ in four steps, with the reaction being exothermic by 17.5 kcal mol⁻¹. However, neither of these two compounds has been found in the ISM. The two-carbon molecule is the radical derived from ethanol (detected in the ISM)⁴⁷ lacking a hydrogen atom. The other reactant in the high-energy isomer of formic acid (HOCOH), which can undergo hydrogen transfer to form the more stable acid.

The reaction between the reactants depicted in Fig. 3c is the most promising and is discussed in more detail below.

Finally, the reaction between the fragments shown in Fig. 3d is the most energetically favorable of all, with an energy difference of 60.6 kcal mol⁻¹ between the reactants and GCA + H˙, leading to the formation of GCA in five steps. Additionally, the small fragment, methylene (CH₂), has been detected in space.⁴⁸ However, the larger fragment is relatively complex, with three oxygen atoms. Nevertheless, the existence of this compound would not be unexpected, as it is structurally similar to ethene-1,2-diol with an additional OH group. Various hypotheses could be proposed regarding its formation.

As shown in Fig. 4, the reaction between hydroxycarbene and HOCHCH₂OH (at 43.5 kcal mol⁻¹), the reactants of Fig. 3c, is barrierless and highly energetically favorable, leading to the formation of glyceraldehyde and atomic hydrogen after three steps. However, the energy profile of Fig. 4 clearly indicates that the most favorable pathway results in the formation of the fragments displayed in Fig. 3a, formaldehyde and the HOCHCH₂OH radical, or to the formation of water and HOCH₂CHCHO. The former represents a peculiar and less efficient pathway for the isomerization of hydroxycarbene into formaldehyde.

Fig. 4 DFT-computed energy profile for the reaction HOCH + HOCHCH₂OH. The energies include scaled harmonic zero-point vibrational energy (ZPE), with the zero of energy set at GCA + H˙, shown at the right of the figure. The diagram is based on the analysis of the lumped networks.

As mentioned above, hydroxycarbene has not been detected in space. One possible reason for this could be that hydrogen transfer between oxygen and carbon, leading to the more stable formaldehyde conformer, occurs relatively easily. This reaction has been previously investigated using higher levels of theory than those employed in the present study.^49,50 Nevertheless, this compound frequently appears in our CRNs and seems to be a highly reactive species, driving several barrierless and highly exothermic reactions. The numerous reaction pathways that the presence of this compound in the ISM would enable motivated us to study its unimolecular kinetics in greater detail.

The HCOH molecule can adopt two possible conformations: syn and anti. The anti conformation can undergo hydrogen transfer via tunneling, leading to the formation of formaldehyde, as shown in Fig. 5. Using both the CCUS and the CVT/SCT models, the rate coefficient for the formation of formaldehyde from hydroxycarbene in the temperature range 10–50 K temperature range varies between 1.3 × 10⁻⁴ and 7.5 × 10⁻³ (see Table S2 in the SI for detailed values). On the other hand, the rate coefficients for the reaction of hydroxycarbene with HOCHCH₂OH radicals are on the order of 7.2–9.4 × 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹. Therefore, the rearrangement of hydroxycarbene depicted in Fig. 5 may potentially compete with its reaction with HOCHCH₂OH, depending on the abundance of these radicals.

Fig. 5 DFT-computed energy profile for the reaction HCOH → CH₂O rearrangement reaction. The energies include scaled harmonic zero-point vibrational energy (ZPE), with the zero of energy set at HCOH. Fig. S3 of the SI shows the corresponding level 3-corrected energies. The diagram is based on the analysis of the lumped networks.

Finally, we also explored other interesting reactions identified within this CRN, beyond those leading to GCA. In this context, we identified a highly efficient formation route for both hydroxymethyl radical and glycolaldehyde (HOCH₂CHO) from hydroxycarbene and the HOCH₂CH₂O radical:

HOCH + HOCH₂CH₂O → CH₂OH + HOCH₂CHO (12)

The rate coefficients for this reaction were computed using the CCUS model and found to range from 7.5 × 10⁻¹⁰ to 1.1 × 10⁻⁹ cm³ molecule⁻¹ s⁻¹ over the 10–100 K temperature interval. See Table S3 and Fig. S4 in the SI for detailed values.

3.3 Analysis of reactions analogous to reaction (2)

The observation that the lowest-energy pathway from GCA–H leads to the formation of formaldehyde and the HOCHCH₂OH radical (vide supra), prompted us to explore other reaction types. Specifically, we looked at reactions analogous to reaction (2):

R3˙ + M → GCA + X˙ (13)

where X˙ is a fragment larger than atomic hydrogen. Eight different reactant combinations were examined, and their reaction energies were computed at the same DFT level employed throughout this study (see Fig. 6). As shown in the figure, the most exothermic and promising reaction occurs between glyoxal (HCOHCO) and the HOCHCH₂OH radical, with an exothermicity of 9.2 kcal mol⁻¹:

HCOHCO + HOCHCH₂OH → GCA + HCO (14)

Fig. 6 DFT-computed reaction energies for various reactions of the type R3˙ + M → GCA + X˙ (reaction (13), see text). In each case, the zero-energy reference corresponds to the products, with the respective reactant pairs indicated in the figure.

Glyoxal can exist in two structural forms: E and Z. Although it has not been detected in the ISM yet, its presence in space has been hypothesized based on the dimerization of HCO radicals.⁵¹ One possible reason for its elusive detection is that the lowest-energy isomer (E) has a vanishing dipole moment, making it undetectable by rotational emission.¹⁰ While we provide rate coefficients for reactions involving glyoxal, we emphasize that these results are intended to support astrochemical modeling efforts and do not constitute evidence that such reactions are responsible for the formation of GCA in the ISM.

The reaction pathways connecting the reactants and products of reaction (14), computed at level 2, are shown in Fig. 7. The corresponding level 3-corrected pathways are presented in Fig. S5.

Fig. 7 DFT-computed energy profile for the reaction HCOHCO + HOCH₂CHOH → GCA + HCO˙. The energies include scaled harmonic zero-point vibrational energy (ZPE), with the zero of energy set at GCA + HCO˙, shown at the right of the figure. Fig. S5 of the SI shows the corresponding level 3-corrected energies. The diagram is based on the analysis of the lumped networks.

The rate coefficients for the formation of GCA and HCO viareaction (14) were computed in this study. The reaction is limited by the initial association step, and the corresponding rate coefficients are provided in Table S4 of the SI. In the 10–100 K temperature range, the rate coefficients for the Z isomer are higher than those for the E isomer and range from 5.4 × 10⁻¹⁰ to 7.9 × 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹.

3.4 Analysis of the C₂H₅O₂ CRN

The radical intermediate HOCHCH₂OH appeared in multiple reaction pathways analyzed above. Consequently, we conducted a more detailed investigation of its corresponding CRN. The energy profile highlighting the most relevant pathways is shown in Fig. 8.

Fig. 8 DFT-computed energy profile for the most relevant pathways of the C₂H₅O₂ network. The energies include scaled harmonic zero-point vibrational energy (ZPE), with the zero of energy set at the HOCHCH₂OH radical intermediate. The diagram is based on the analysis of the lumped networks.

In the C₂H₅O₂ CRN of Fig. 8, we can distinguish two large regions separated by a high barrier (TS93). On the left side, we observe both the E and Z isomers of ethene-1,2-diol and atomic hydrogen (at 35.2 and 39.0 kcal mol⁻¹). The reaction between any of the two isomers and hydrogen has a barrier, which rules it out as a possible pathway for the formation of other molecules. However, ethene-1,2-diol can be formed from the reaction between oxirane and hydroxyl radical (located in the center of the figure at 46.6 kcal mol⁻¹). However, we can assume that this reaction would exclusively lead to the most favorable route, resulting in formaldehyde and the hydroxymethyl radical (at 19.6 kcal mol⁻¹).

Another interesting aspect of the profile is the reaction between ethenol (CH₂CHOH) and hydroxyl radical (at 30.3 kcal mol⁻¹). These two species, both detected in the ISM, react without a barrier to form the HOCHCH₂OH radical. Once this radical is formed, the only two possibilities are either decomposition back into the initial reactants or stabilization via spontaneous radiation. This result supports the consideration of the HOCHCH₂OH radical as a potential reactant, since, even if it has not been detected, its existence in the ISM would be expected.

On the right side of the energy profile, the other isomer of acetic acid, methyl formate (CH₃OCHO), appears along with atomic hydrogen at an energy of 16.2 kcal mol⁻¹. Although its formation would be energetically viable from methoxy radical (CH₃O) and hydroxycarbene (at 78.8 kcal mol⁻¹), it would be expected that formic acid (HCOOH) and methyl radical (at −5.5 kcal mol⁻¹) would be the predominant products.

It should be noted here that a more comprehensive study by Barone's group on the same reaction network indicates that H-abstraction in the reaction between the OH radical and vinyl alcohol primarily yields the CH₂CHO radical and water.²⁵ Our results cannot confirm this conclusion, as they are based on the analysis of a lumped network in which reactions of the type A + B → C + D are not considered.

4 Conclusions

This study presents a comprehensive computational investigation of the chemical reaction networks associated with the gas-phase chemistry of C₃H₆O₃, C₃H₇O₃, and C₂H₅O₂ species in the context of the interstellar medium. Using the automated reaction discovery tool AutoMeKin in combination with DFT calculations and statistical rate theory, we identified multiple energetically favorable pathways relevant to the formation and decomposition of GCA under astrochemically relevant conditions.

Our analysis of the C₃H₆O₃ CRN reveals several plausible formation routes for GCA, including radical–radical and radical–molecule association reactions. Notably, pathways involving hydroxycarbene, HOCHCH₂OH, formyl and hydroxymethyl radicals proceed without barriers and exhibit significant exothermicity. Despite their kinetic feasibility, these reactions typically yield highly energized GCA* molecules, which are prone to rapid unimolecular decomposition before they can be stabilized via radiative cooling. Predicted decomposition products include water, 2-hydroxyprop-2-enal, glycolaldehyde, formaldehyde, and ethene-1,2-diol; all species that have either been observed in the ISM or are considered highly probable.

The C₃H₇O₃ CRN further supports the hypothesis that GCA functions as a transient intermediate. Even within this network, the most accessible pathways favor decomposition over the accumulation of stable GCA. We also explored alternative reactions analogous to GCA formation that involve the departure of larger radicals instead of atomic hydrogen. Among these, the reaction between glyoxal and the HOCHCH₂OH radical is particularly notable for its barrierless and exothermic profile. However, the detection of glyoxal in the ISM remains challenging due to its near-zero dipole moment.

Our complementary analysis of the C₂H₅O₂ network underscores the chemical significance of the HOCHCH₂OH radical, which emerges repeatedly as an intermediate in multiple reaction schemes. Although this species has not yet been detected in the ISM, its formation is supported by barrierless reactions involving known interstellar molecules such as ethenol and hydroxyl radicals.

Taken together, these results suggest that while the gas-phase formation of GCA is chemically feasible under ISM conditions, it is likely outcompeted by pathways leading to other, more stable complex organic molecules (COMs). GCA may exist only fleetingly in interstellar environments, serving as a reactive intermediate en route to more persistent species. The data generated in this study, including structural information, reaction mechanisms, and rate coefficients, provide essential input for astrochemical models and may inform future observational strategies.

Overall, this work highlights the power of automated CRN exploration in astrochemistry and reinforces the notion that molecular complexity in space emerges not only through stepwise buildup but also through dynamic networks of reactive intermediates, shaped by the unique physical conditions of the interstellar medium.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

Data for this article, including the chemical reaction networks of the three studied systems are available at https://doi.org/10.5281/zenodo.15267696.

Supplementary information: convergence of product abundances with iteration number, computed rate coefficients and plots of selected barrierless pathways. See DOI: https://doi.org/10.1039/d5cp01690h.

Acknowledgements

This work was partially supported by the Consellería de Cultura, Educación e Ordenación Universitaria (Centro singular de investigación de Galicia acreditación 2019–2022, ED431G 2019/03) and the European Regional Development Fund (ERDF). ALS thanks Xunta de Galicia for financial support through a predoctoral grant. We also thank the Centro de Supercomputación de Galicia (CESGA) for providing access to their computational facilities.

References

1. A. Venturini and J. González, ChemPlusChem, 2024, 89, e202300388.
2. Y. Layssac, A. Gutiérrez-Quintanilla, T. Chiavassa and F. Duvernay, Mon. Not. R. Astron. Soc., 2020, 496, 5292–5307.
3. A. Fukutome, H. Kawamoto and S. Saka, ChemSusChem, 2016, 9, 703–712.
4. S. Wu, D. Shen, J. Hu, R. Xiao and H. Zhang, J. Anal. Appl. Pyrolysis, 2013, 101, 79–85.
5. I. Jiménez-Serra, J. Martín-Pintado, V. M. Rivilla, L. Rodríguez-Almeida, E. R. Alonso Alonso, S. Zeng, E. J. Cocinero, S. Martín, M. Requena-Torres, R. Martín-Domenech and L. Testi, Astrobiology, 2020, 20, 1048–1066.
6. J. K. Jørgensen, C. Favre, S. E. Bisschop, T. L. Bourke, E. F. van Dishoeck and M. Schmalzl, Astrophys. J., Lett., 2012, 757, L4.
7. A. Coutens, M. V. Persson, J. K. Jørgensen, S. F. Wampfler and J. M. Lykke, Astron. Astrophys., 2015, 576, A5.
8. F. Goesmann, H. Rosenbauer, J. H. Bredehöft, M. Cabane, P. Ehrenfreund, T. Gautier, C. Giri, H. Krüger, L. Le Roy, A. J. MacDermott, S. McKenna-Lawlor, U. J. Meierhenrich, G. M. M. Caro, F. Raulin, R. Roll, A. Steele, H. Steininger, R. Sternberg, C. Szopa, W. Thiemann and S. Ulamec, Science, 2015, 349, aab0689.
9. M. Meringer, C. Giri and H. J. Cleaves, II, ACS Earth Space Chem., 2018, 2, 1256–1261.
10. V. M. Rivilla, L. Colzi, I. Jiménez-Serra, J. Martín-Pintado, A. Megías, M. Melosso, L. Bizzocchi, Á. López-Gallifa, A. Martínez-Henares, S. Massalkhi, B. Tercero, P. de Vicente, J.-C. Guillemin, J. García de la Concepción, F. Rico-Villas, S. Zeng, S. Martín, M. A. Requena-Torres, F. Tonolo, S. Alessandrini, L. Dore, V. Barone and C. Puzzarini, Astrophys. J., Lett., 2022, 929, L11.
11. C. Puzzarini, Front. Astron. Space Sci., 2022, 8, 811342.
12. C. Robertson, R. Hyland, A. J. D. Lacey, S. Havens and S. Habershon, J. Chem. Theory Comput., 2021, 17, 2307–2322.
13. D. Garay-Ruiz, M. Álvarez-Moreno, C. Bo and E. Martínez-Núñez, ACS Phys. Chem. Au, 2022, 2, 225–236.
14. E. Martínez-Núñez, J. Comput. Chem., 2015, 36, 222–234.
15. E. Martínez-Núñez, Phys. Chem. Chem. Phys., 2015, 17, 14912–14921.
16. E. Martínez-Núñez, G. L. Barnes, D. R. Glowacki, S. Kopec, D. Peláez, A. Rodríguez, R. Rodríguez-Fernández, R. J. Shannon, J. J. P. Stewart, P. G. Tahoces and S. A. Vazquez, J. Comput. Chem., 2021, 42, 2036–2048.
17. M. Castiñeira Reis, E. Martínez Núñez and A. Fernández Ramos, Sci. Adv., 2024, 10, eadq4077.
18. B. Ballotta, T. D. Marforio, S. Rampino, E. Martínez-Núñez, V. Barone, M. Melosso, A. Bottoni and L. Dore, ACS Earth Space Chem., 2023, 7, 1172–1180.
19. J. C. Ramal-Olmedo, C. A. Menor-Salván, A. Miyoshi and R. C. Fortenberry, Astron. Astrophys., 2023, 672, A49.
20. J. C. Ramal-Olmedo, C. A. Menor-Salván and R. C. Fortenberry, Astron. Astrophys., 2021, 656, A148.
21. L. Guerrero-Méndez, A. Lema-Saavedra, E. Jiménez, A. Fernández-Ramos and E. Martínez-Núñez, Phys. Chem. Chem. Phys., 2023, 25, 20988–20996.
22. M. J. Wilhelm, E. Martínez-Núñez, J. González-Vázquez, S. A. Vázquez, J. M. Smith and H.-L. Dai, Astrophys. J., 2017, 849, 15.
23. J. C. Ramal-Olmedo, C. A. Menor-Salván, A. Miyoshi and R. C. Fortenberry, ACS Earth Space Chem., 2023, 7, 1694–1712.
24. B. Ballotta, E. Martínez-Núñez, S. Rampino and V. Barone, Phys. Chem. Chem. Phys., 2023, 25, 22840–22850.
25. B. Ballotta, E. Martínez-Núñez, S. Rampino and V. Barone, ACS Earth Space Chem., 2023, 7, 1467–1477.
26. J. J. P. Stewart, J. Mol. Model., 2013, 19, 1–32.
27. V. K. Prasad, Z. Pei, S. Edelmann, A. Otero-de-la-Roza and G. A. DiLabio, J. Chem. Theory Comput., 2022, 18, 151–166.
28. A. M. Parameswaran, A. Fernández-Ramos and D. G. Truhlar, J. Phys. Chem. A, 2024, 128, 10673–10679.
29. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision B.01, Gaussian, Inc., Wallingford CT, 2010.
30. D. T. Gillespie, J. Comput. Phys., 1976, 22, 403–434.
31. R. G. Gilbert and S. C. Smith, Theory of Unimolecular and Recombination Reactions, Blackwell Science Publishing, Oxford, 1990.
32. T. Beyer and D. Swinehart, Commun. ACM, 1973, 16, 379.
33. E. Martinez-Nunez, A. Lema-Saavedra and A. Fernandez-Ramos, 2025, C₂H₅O₂, C₃H₆O₃ and C₃H₇O₃ chemical reaction networks [Data set], Zenodo DOI:10.5281/zenodo.15267696.
34. C. Y. Legault, CYLview, 1.0b, Université de Sherbrooke, 2009, (https://www.cylview.org).
35. L. G. Gao, J. Zheng, A. Fernández-Ramos, D. G. Truhlar and X. Xu, J. Am. Chem. Soc., 2018, 140, 2906–2918.
36. D. González, A. Lema-Saavedra, S. Espinosa, E. Martínez-Núñez, A. Fernández-Ramos, A. Canosa, B. Ballesteros and E. Jiménez, Phys. Chem. Chem. Phys., 2022, 24, 23593–23601.
37. Y. Georgievskii and S. J. Klippenstein, J. Chem. Phys., 2005, 122, 194103.
38. A. Fernandez-Ramos, B. A. Ellingson, B. C. Garrett and D. G. Truhlar, Reviews in Computational Chemistry, 2007, pp. 125–232.
39. R. T. Skodje, D. G. Truhlar and B. C. Garrett, J. Phys. Chem., 1981, 85, 3019–3023.
40. S. J. Klippenstein, Y. C. Yang, V. Ryzhov and R. C. Dunbar, J. Chem. Phys., 1996, 104, 4502–4516.
41. I. M. Alecu, J. Zheng, Y. Zhao and D. G. Truhlar, J. Chem. Theor. Comput., 2010, 6, 2872–2887.
42. M. Ohishi, S.-i Ishikawa, T. Amano, H. Oka, W. M. Irvine, J. E. Dickens, L. M. Ziurys and A. J. Apponi, Astrophys. J., 1996, 471, L61.
43. P. R. Joshi and Y.-P. Lee, J. Am. Chem. Soc., 2024, 146, 23306–23320.
44. L. E. Snyder, J. M. Hollis and B. L. Ulich, Astrophys. J., 1976, 208, L91–L94.
45. L. E. Snyder, D. Buhl, B. Zuckerman and P. Palmer, Phys. Rev. Lett., 1969, 22, 679–681.
46. J. M. Hollis, F. J. Lovas and P. R. Jewell, Astrophys. J., 2000, 540, L107.
47. B. Zuckerman, B. E. Turner, D. R. Johnson, F. O. Clark, F. J. Lovas, N. Fourikis, P. Palmer, M. Morris, A. E. Lilley, J. A. Ball, C. A. Gottlieb, M. M. Litvak and H. Penfield, Astrophys. J., 1975, 196, L99–L102.
48. J. M. Hollis, P. R. Jewell and F. J. Lovas, Astrophys. J., 1995, 438, 259.
49. V. G. Kiselev, S. Swinnen, V. S. Nguyen, N. P. Gritsan and M. T. Nguyen, J. Phys. Chem. A, 2010, 114, 5573–5579.
50. P. R. Schreiner, H. P. Reisenauer, F. C. Pickard Iv, A. C. Simmonett, W. D. Allen, E. Mátyus and A. G. Császár, Nature, 2008, 453, 906–909.
51. P. M. Woods, B. Slater, Z. Raza, S. Viti, W. A. Brown and D. J. Burke, Astrophys. J., 2013, 777, 90.

---