Gaussian 16 quantum chemical software was employed for all the ab initio//DFT calculations. (Gaussian 16 et al., 2016). We employed the unrestricted double-hybrid density functional method, specifically the B2PLYP functional, (Grimme, 2006), with augmented triple-zeta basis sets (Dunning, 1989) (aug-cc-pVTZ and aug-cc-pV (T+d)Z, the latter used only for the sulfur system) to optimize the structures of reactants, pre-reactive complex, intermediates, transition states, and products. The B2PLYP functional was supplemented with Grimme’s DFT-D3 zero-damping correction to address long-range London dispersion interactions. (Grimme et al., 2011). The B2PLYP-D3 double-hybrid functional is particularly well-suited for investigating noncovalent interactions in transition states, intermediates, and post-intermediates, especially in systems involving hydrogen bonding, as suggested in the earlier study. (Vazart et al., 2016). Several research groups have successfully used this method to predict the astrochemical formation of nitrogen-containing compounds such as methylamine and formyl cyanide^. (Tonolo et al., 2020; Puzzarini et al., 2020). To confirm the nature of the saddle point, frequency analysis was performed, revealing a single imaginary frequency for the transition states and all positive frequencies for the reactants, intermediates and products. Vibrational analysis was performed for each optimized species to include the zero-point vibrational energy (ZPE). Energy accuracy was improved by performing single-point calculations at the CCSD(T) (Raghavachari et al., 1989)/aug-cc-pVTZ basis set, using geometries optimized at the B2PLYP-D3 level. The combination of CCSD(T) and B2PLYP-D3-optimized methods has been used previously, achieving accuracy within ∼1 kcal/mol. (Ballotta et al., 2021). The zero-point energies (ZPEs) obtained using the B2PLYP-D3 functional are suitable for interstellar medium (ISM) chemistry; however, there may be an uncertainty of 1–2 kcal/mol in the computed energies, which can affect predictions for low-temperature kinetics. B2PLYP-D3 was selected over more conventional hybrid functionals because it offers a consistent and well-established balance between accuracy and computational cost, particularly for systems where both dynamical and static correlation are non-negligible. (Grimme, 2006; Dunning, 1989). The double-hybrid nature of B2PLYP, which includes a perturbative second-order correlation component, improves the description of subtle electronic effects that are critical for accurately modeling low-frequency vibrational modes and tunneling pathways. Furthermore, B2PLYP-D3 has demonstrated superior performance in reproducing vibrational frequencies and barrier heights in similar systems compared to standard hybrid functionals, especially in the presence of non-covalent interactions and anharmonic contributions. (Vazart et al., 2016; Barone et al., 2015). This aligns with the objectives of our study, which emphasize the accurate characterization of low-frequency vibrational modes associated with quantum tunneling. Earlier work by Papamokos and Demetropoulos supports the use of the hybrid functional PW91XC for amide systems and addresses low-frequency vibrational issues, as discussed in their study. (Papamokos and Demetropoulos, 2004a; Papamokos and Demetropoulos, 2004b). The increased accuracy offered by double-hybrid methods such as B2PLYP-D3 makes them particularly suitable for the current investigation, which focuses on low-temperature tunneling-relevant vibrational features and dispersion corrected zero-point energies. Our choice thus reflects a targeted effort to improve the reliability of both geometries and harmonic frequencies in a regime where conventional hybrid functionals may not be sufficiently accurate.

To investigate how the single-reference wave function qualitatively contributes, T1 diagnostic calculations were done at the CCSD(T)/aug-cc-pVTZ level. The resulting T1 diagnostic values were ≤0.03, which is within the accepted range (≤0.04) for a single-reference wave function. (Rienstra-Kiracofe et al. 2000). To evaluate spin contamination, the spin expectation value ⟨S²⟩ was calculated, yielding values in the range of ∼0.75–0.77, indicating that spin contamination is negligible.

All chemical kinetics calculations were performed using the MultiWell program suite. (Barker, 2001; Barker, 2009; Barker et al., 2016). The rate constants were calculated using RRKM/ME theory, where the energy- and angular momentum-dependent unimolecular rate coefficient, k(E, J) (Barker, 2001; Barker, 2009; Barker et al., 2016) is given by Equation 1: k E , J = L ≠ h G ≠ E − E 0 , J , J ρ E , J (1)

Here, L ⁺ is the reaction path degeneracy, h is Planck’s constant, ρ(E,J) is the rovibrational density of states of the reactant, G ⁺ (E - E ₀ ,J, J) is the transition state sum of states, and E ₀ ,J represents the critical energy threshold for the reaction, which is angular momentum-dependent. The reaction path degeneracy was determined from symmetry considerations and optical isomer counts. The DenSum module within MultiWell was used to compute these sums and densities of states as discussed in the MultiWell manual. (Barker, 2001; Barker, 2009; Barker et al., 2016). For efficiency, the two-dimensional terms were integrated over angular momentum to yield one-dimensional forms (Equations 2, 3): G ≠ E − E 0 , 0 = ∑ J = 0 J m a x G ≠ E − E 0 , J , J (2) ρ E , J = ∑ J = 0 J m a x ρ i E , J (3)

The 1-D rate constants k(E) were calculated as given in Equation 4: k E = L ≠ h G ≠ E − E 0 , 0 ρ E (4) where E ₀ , ₀ includes zero-point energy and J = 0 centrifugal corrections. Rotational motion was modeled using the symmetric top approximation (A > B = C), treating the K-rotor as active and the perpendicular 2D-rotor (B = C) as adiabatic. The rotational constants were derived from B2PLYP-D3/aug-cc-pVTZ optimized geometries and vibrational frequencies, and the Stein–Rabinovitch version of the Beyer–Swinehart algorithm was employed for statistical summations. To calculate the pressure-dependent rate constants, N₂ bath gases were used with an approximate value of the energy transfer process < Δ E > d o w n = 200 × T 300 0.85 cm⁻¹. (Ali and Balaganesh, 2023) In the MultiWell software, master equation simulations utilized double arrays containing 500 elements, each representing energy intervals of 10 cm⁻¹. The quasi-continuum region was modeled up to an energy limit of 85,000 cm^-1. Simulations at each specified temperature and pressure began with a chemical activation energy distribution, which is particularly suited for modeling recombination processes. A total of 10⁶ stochastic trials were performed, with each trial simulating a time span equivalent to the average duration of 1,000 molecular collisions. The pressure-dependent forward total rate constants, k _f for CS + NH₂ reactions, have been computed using, k f = Γ   K e q × k ∞ u n i 1 − f PR (5)

The Γ Eckart asymmetric tunneling (EAT) correction was applied as implemented in the MultiWell master equation (ME) code, with the tight transition state explicitly calculated, f is the fraction of the chemical reaction going back to the respective reactive species, and k ∞ u n i is the rate constant at the high-pressure limit.

The chemical kinetics of the pre-reactive complex (PRC) formation CS….NH₂, which is assumed to play an important role in ISM condition, has not been investigated. The rate constants for CS + NH₂ were calculated using a combination of variational transition state theory (VTST) and µVTST. (Ali and Balaganesh, 2023; Ali et al., 2016; Ali et al., 2021). The variational transition state theory (VTST) was used to calculate the rate constants as given in Equation 6 and Equation 7: k G T T , s = L ≠ × k B T h Q T S ≠ T , s Q R T exp − E M E P s k B T (6) k V T S T T = m i n s   k G T T , s = k G T ( T , s V T S T T (7) where k G T T , s and k V T S T T are the rate constants of generalized and variational transition state theory, respectively, E M E P is the zero-point corrected barrier height. Based on VTST calculations, µVTST calculations were carried out as discussed in our previous work. (Ali and Balaganesh, 2023; Ali et al., 2016; Ali et al., 2021) The ktools module within the MultiWell program was used to perform the VTST and µVTST calculations, as described in the MultiWell manual. (Barker, 2001; Barker, 2009; Barker et al., 2016) The equilibrium constant (K _eq ) for the formation of bimolecular complexes were calculated using THERMO code as given in Equation 8: K e q = Q Int Q R exp − E Int − E R k B T . (8)

Barrier heights are important parameters for determining astrochemical reactions in the ISM. To ensure consistency with other theoretical methods, several ab initio and DFT calculations (B2PLYP-D3 (Grimme, 2006; Grimme et al., 2011), B3LYP (Becke, 1993; Lee et al., 1988), M06-2X (Zhao and Truhlar, 2008), and ωB97XD (Chai and Head-Gordon, 2008)) were performed using the CCSD(T)/aug-cc-pVTZ//B2PLYP-D3/aug-cc-pVTZ, CCSD(T)/aug-cc-pVTZ//B3LYP/aug-cc-pVTZ, CCSD(T)/aug-cc-pVTZ//M06-2X/aug-cc-pVTZ, and CCSD(T)/aug-cc-pVTZ//ωB97XD/aug-cc-pVTZ methods for both amide and thioamide schemes. A combination of CCSD(T) and several DFT-based optimization methods has been successfully used in our past work, showing good consistency with experimental findings for similar systems. (Ali and Balaganesh, 2022; Ali et al., 2019; Ali et al., 2018; Ali, 2020; Ali et al., 2022; Ali, 2019; Ali and Saswathy, 2024; Ali and Barker, 2015). The results in Table 1 show that our calculations align with different levels of theory, with a deviation of less than 0.5 kcal/mol for the pre-reactive complex and transition state. The energy analysis in Table 1 clearly indicates that CCSD(T)//B2PLYP-D3 values are accurate compared to the other three DFT methods. As suggested in previous work, (Puzzarini et al., 2020), CCSD(T)//B2PLYP-D3 provides accurate energies compared to the “Cheap” composite method. Therefore, all minimum energy pathway (MEP) and kinetics calculations were performed at the CCSD(T)//B2PLYP-D3 level.

To further validate our results on the formation of PRC in the sulfur system, we employed several quantum composite methods, including CBS-QB3 (Montgomery et al., 2000), G4 (Curtiss et al., 2007), and W1BD (Barnes et al., 2009) methods were employed (see ESI Supplementary Figure SI,13). These calculations indicate that our results are consistent across different levels of theory, with a difference of less than 0.5 kcal/mol for the pre-reactive complex and transition state. Therefore, we believe the results presented in this paper are reasonably accurate for astrochemical implications.

Possible reaction pathways for the formamide (Path 1 and Path 2), urea (Path 3), and protonated urea (Path 4) are shown in Scheme 2.

Figure 1 displays the PES for formamide and urea formation. The optimized structures of reactants, pre-reactive complex, intermediates, transition states and products, along with energy values, are shown in ESI Supplementary Figure S1. The Cartesian coordinates for all these species are tabulated in ESI Supplementary Table S1. When the ^·NH₂ radical attacks the carbon atom of CO, a pre-reactant complex (OC· ··NH₂) is formed (see ESI Supplementary Figure S1). The computed energy for the formation of PRC is nearly −0.2 kcal/mol (at CCSD(T)/aug-cc-pVTZ//B2PLYP-D/aug-cc-pVTZ), which is in good agreement with the energy computed energy at CCSD(T)/aug-cc-pVTZ//ωB97X-D/aug-cc-pVTZ level. This result is attributed to a weak interaction between C and N atoms. To form a C–N bond between N of ^·NH₂ and C of CO, both molecules come closer to each other via a transition state (TS-1), incurring a significant energy penalty. The barrier height for C-N bond formation is approximately ∼6 kcal/mol, leading to the formation of the intermediate Int-1 (^·CONH₂). The energy value is well-supported by calculations performed at the CCSD(T)/aug-cc-pVTZ level using ωB97X-D optimized geometries. The Int-1 (^·CONH₂) can react with hydrogen molecules to form Int-2 (H-H … CONH₂), with an energy barrier of 16.6 kcal/mol, ultimately dissociates via transition state TS-2 to produce formamide (HCONH₂) + H. As shown in path-2, Int-1 can undergo isomerization via H migration, forming Int-3 (HCONH) through a transition state (TS-3), which has a high barrier of 51.5 kcal/mol relative to Int-1 (^⋅CONH₂). Additionally, the transition state (TS-4) involves the addition of hydrogen molecules to form Int-4 (H-H … CONH₂). Due to the high energy barrier and the formation of endothermic intermediates, Path-2 is less feasible compared to Path-1, making such a reaction more relevant under combustion conditions rather than ISM environments.

As shown in Figure 1, the ^·NH₂ radical reacts with CO to form Int-1 (^·CONH₂), which then reacts with H₂ to produce formamide (HCONH₂), following a sequential two-body process. In this study, we consider that in the presence of H₂, CO, and ^·NH₂ molecules, simultaneous three-body collisions are highly unlikely. Instead, the reaction is expected to proceed through the formation of two-body complexes, which subsequently collide with a third species to form three-body complexes and products. A sequential two-body collision mechanism is therefore proposed, as illustrated in ESI Page S18. Such reaction hypotheses for the ISM and EA will be validated through kinetic calculations and vice versa.

To understand the formation of urea (NH₂CONH₂), Int-1 can also react with NH₃ molecules to form CONH₂· ·NH₃ (Int-5), which is more stable than Int-1. This increased stability is due to the formation of two hydrogen bonds from ammonia molecules (see ESI Supplementary Figure SI,1). As shown in Figure 1, the Int-5 involves two pathways, Path 3 and Path 4, to generate NH₂CONH₂ molecules. Path 3 involves the formation of a second C-N bond of Int-5, which then dissociates via transition state TS-5, leading to the production of NH₂CONH₂ and H radical with a barrier of 37.1 kcal/mol. Another possibility involves the simultaneous formation of a C ̶ N bond and protonation of CONH₂, resulting in the formation of the protonated urea (Product-4). The increased stability of this product is due to the formation of three hydrogen bonds. Based on the barrier height calculations, which are relatively high in all four proposed pathways, we believe the formation of NH₂COH and NH₂COH₂ is unlikely to occur in the gas-phase ISM.

Possible reaction pathways for the thioformamide (Path 1 and 2), thiourea (Path 3), and protonated thiourea (Path 4) are shown in Scheme 3.

Figure 2 displays the PES for the formation of thioformamide (NH₂CSH) and thiourea (NH₂CSNH₂), including zero-point energy corrections. Carbon sulfide (CS) is the S-analog of carbon monoxide (CO), with the primary difference between them arising from their dipole moments. The C=S (μ = 1.912 D) has a significantly higher dipole moment than the C=O (μ = 0.113 D), which is reflected in the longer bond length of CS (1.52 Å) than CO 1.12 Å. The formation of a pre-reactant complex (CS····NH₂, PRC_s) is a barrierless process with an energy of −1.7 kcal/mol (w.r.t. Reactants). This stationary point is 0.3 kcal/mol lower than that of the transition state TS-1_s to form a submerged barrier relative to the reactants at infinite separation, facilitating N-S bond formation. This is because carbon sulfide has a higher dipole moment, bond length, and polarizability, leading to more reactivity with NH₂ radicals.

The Int-1_s (^·CSNH₂) can react with hydrogen molecules to form Int-2_s (H₂ … CSNH₂) which then dissociates via transition state TS-2_s to produce thioformamide (HCSNH₂) + H radical, with an energy barrier of 17.7 kcal mol^-1, as seen in Path 1. In Path 2, Int-1 (see ESI Supplementary Figure S2) undergoes isomerization to form Int-3_s (HCSN·H), which then undergoes hydrogenation, followed by breaking the H-H bond, resulting in the formation of thioformamide and a H radical. As shown in ESI Supplementary Figure SI,2, Path 3 and Path-4 lead to two different transition states (TSs) for thiourea (TS-5_s) and protonated thiourea (TS-6_s). In this process, ammonia reacts with C of Int-1_s, forming thiourea and a H radical. Another possible reaction pathway (Path 4) involves Int-5_s reacting with an ammonia molecule, simultaneously leading to C–N bond formation and protonation (H⁺) of Int-1_s. This reaction results in the formation of protonated thiourea, which is more stable than thiourea. As shown in Figure 2, since the first barrier is low, the other barriers are relatively lower, and the exothermic nature of the products, Path 1, Path 3, and Path 4, is likely to be feasible under ISM conditions. On the other hand, Path 2 has a barrier height of ∼50 kcal/mol, making this pathway infeasible under ISM conditions.

The dipole moment is a crucial parameter for detecting molecules in the ISM, as it can be analyzed using microwave spectroscopy (Kroto, 1981) and quantum chemical calculations. Molecules with smaller dipole moments exhibit weak emissions and are difficult to detect, even if they are abundant in the ISM. The dipole moments (in debye) for stable amide species are as follows: CO (0.113), CONH₂ (3.693), HCONH (3.388), HCOONH₂ (3.975), and NH₂CONH₂ (4.384). For stable thioamide species, the values are CS (1.912), Int-1_s (^·CSNH₂) (3.786), Int-3_s (^·HCSNH) (2.209), thioformamide (HCSNH₂) (4.528), and thiourea (NH₂CSNH₂) (5.586). Molecules with higher dipole moments are more likely to be detected in the ISM. Therefore, sulfur-containing compounds such as Int-1_s (^·CSNH₂), thioformamide (HCSNH₂), and thiourea (NH₂CSNH₂) may also be present. However, experimental evidence is needed to confirm this prediction.

To identify the transition state for the CS⋯NH₂ → CS + NH₂ dissociation, the minimum energy path (MEP) was computed by performing constrained optimizations along the bond distance, incorporating zero-point energy corrections for the vibrational modes orthogonal to the forming bond, as illustrated in Figure 3. The MEP was determined using the B2PLYP-D3/aug-cc-pVTZ level of theory, performing constrained optimization at S-N distances ranging from 3.8 Å to 9.0 Å in increments of 0.1 Å.

The reaction path illustrated in Figure 3 was used to calculate VTST rate constants across a temperature range of 10 K–100 K are shown in Figure 4.

Based on the VTST calculation, the µVTST rate constants range from approximately 10^–8 to 10^–11 cm³ molecule⁻¹ s^-1 between 10K and 100K, as shown in Figure 5. The tabulated rate constant values are provided in ESI Supplementary Table SI,7. As highlighted in earlier research, barrierless reactions can exhibit a negative temperature dependence, which supports a capture-type mechanism facilitated by long-range intermolecular forces. These forces grow more influential at lower temperatures as thermal energy diminishes. Our results reflect a similar trend, comparable to the NH₂ + NO reaction. (Klemperer, 2006).

Figure 6 presents the rate constants for the CS + NH₂ reaction over a pressure range from 10^–7 bar–1 bar at various temperatures. The figure reveals that the greatest disparity between the two limits occurs around 90 K, with a difference of nearly three orders of magnitude. In contrast, at 10 K, the difference between the two regimes is approximately a factor of 5. To gain a deeper understanding of the formation of thioformamide and thiourea, the rate constants for these channels were also calculated for temperatures ranging from 10 K to 30 K at a very low pressure of approximately 10^–10 bar. Our calculations show that at low temperatures (<30 K) and very low pressures, Int-1_s forms efficiently. However, at high pressures (>10^–7 bar) and temperatures (>30 K), the reaction tends to revert to the reactants, making the formation of thioformamide and thiourea nearly negligible. Our analysis suggests that thioformamide (NH₂CS) and thiourea (NH₂CSNH₂) can be formed from Int-1_s + H₂ and Int-1_s + NH₃, respectively, under conditions of very low temperature and pressure. We also considered the role of competitive reactions, such as Radiative Association (RA), which are significant in interstellar medium (ISM) conditions. The forward reaction rates at low temperatures and very low pressures are faster than typical RA reactions, which generally exhibit rate constants well below ∼10^–11 cm³ molecules^-1 s^-1. In other words, the reactive intermediate (Int-1_s) will interact with “third body” species such as H₂ or NH₃, leading to the formation of NH₂CSH and NH₂CSNH₂. From the data in Figure 6, the low-pressure analysis (10^–7 bar, typical of experimental conditions) shows that under ISM-like pressures (P < 10^–14 bar, corresponding to a molecular density of approximately 10⁶ molecules/cm³), the reactive intermediate (Int-1_s) can still form effectively. This is because the pressure effect is almost negligible.

In the case of formamide (HCONH₂) and urea (NH₂CONH₂) formation, all the precursor molecules, such as CO, ^·NH₂, H₂, and NH₃, have been detected in ISM. However, all transition states (TS-1 to TS-6) involved in Path 1 and Path 2 (formamide formation), Path 3 (urea formation), and Path 4 (protonated urea formation) are infeasible due to the energy barriers. For thioformamide and thiourea formation, it is important to note that the precursors used in this study, CS, ^·NH₂, H₂, and NH_3, have been observed in the ISM. However, intermediates such as Int-1_s (^·CSNH₂) and Int-3_s (HCSN·H), as well as products like thioformamide (HCSNH₂₎ and urea (NH₂CSNH₂), have not yet been detected. The transition states (TS-1_s to TS-6_s) involved in Path 1 (thioformamide formation), Path-3 (thiourea formation), and Path 4 (Protonated thiourea formation) are feasible under interstellar conditions due to their lower barrier heights and kinetic feasibilities (see Figures 2, 6).

Following the approach suggested in previous studies (Ballotta et al., 2021; Klippenstein et al., 1996), the effective association reaction in the low-pressure limit was calculated as k e f f = k f   k r k b + k r , where k _f is the forward association rate constant (cm³ molecule⁻¹ s^-1), k_r (s^-1) is the radiative stabilization rate constant and k_b (s^-1) is the back-dissociation rate constant. To assess the accuracy of radiative stabilization, the forward rate constant was calculated using the µVTST_RRKM/ME method, as discussed in section 3.4. k_f varies from 10^–8 cm³ molecule⁻¹ s^-1 to 10^–11 cm³ molecule⁻¹ s^-1 from the temperature range of 10–100K. As suggested in previous studies (Dalgarno, 1987; Van Dishoeck, 2014; Douglas et al., 2024), and anticipated in our calculation, the k_r is expected to be higher than the backward reaction (k_b). In that case, the effective association rate for this process is almost expected to be the same as the forward reaction (k_r >> k_b, the k e f f = k f ) . This analysis also agrees with the literature value of similar ISM product formation under low temperature and low-pressure conditions, which is in the range of 10^–9 to 10^–11 cm³ molecule⁻¹ s^-1. (Raghavachari et al., 1989; Klippenstein et al., 1996; Dalgarno, 1987).

The above analysis is further supported by the high dipole moments of each species in the sulfur system. The ISM is characterized by extremely low density (10²–10⁶ atoms, molecules, and ions per cm³) and temperatures ranging from 10K to 100 K). (Klemperer, 2006). Due to these conditions, chemical reactions in the ISM are typically barrierless (or involve submerged TSs) and lead to exothermic products. (Van Dishoeck, 2014; Douglas et al., 2024; Jensen, 2007). On Earth, gas phase reactions typically involve the formation of pre-reactive complexes through weak interactions between reactants, often via two-body collisions, with occasional involvement of a third body. However, in interstellar chemistry, three-body collisions are highly improbable due to the extremely low number density, cold temperatures, and long reaction time scales. (Herbst, 1985). For a dense interstellar cloud with a number density of 10⁵ cm⁻³, the typical chemical reaction timescales are: reactions involving “inert” neutral species ∼10⁵ years (or even longer in some cases), reactions of reactive neutral species ∼100 years, reactions of molecular ions, 1 hour–100 years. This aspect has been well-documented in the literature for several years. (Herbst and Klemperer, 1976).

We estimated the timescale for the formation of Int-1s using an approximate association rate of k∼10^–8 to 10^–13 cm³ s^-1 over a temperature range of 10–100 K, under extremely low-pressure conditions. from 10 K to 100K at extremely low pressure for the formation of Int-1s. These calculations assume particle densities ranging from 10² to 10⁶ cm⁻³, representative of diffuse and molecular clouds.”.

At 10 K, τ∼1/(k⋅n_reactant) = (10^–8 cm³ s^-1*10² cm^-3) = 10⁶ sec = ∼ 12 days.

At 100K τ∼1/(k⋅n_reactant) = (10^–13 cm³ s^-1*10⁶ cm^-3) = 10⁷ sec = 115 days.

Based on the above facts, we believe the reaction mechanism proposed in this work aligns with the current understanding of interstellar gas-phase chemistry.

We explicitly acknowledge that, to date, neither thioformamide nor thiourea has been detected in the interstellar medium (ISM). Although our quantum chemical and kinetic calculations indicate that their formation is feasible under cold gas-phase conditions, there is currently no observational spectroscopic confirmation. This limitation is now clearly emphasized, and our study is framed as a theoretical foundation to motivate targeted astronomical searches—particularly in sulfur-rich interstellar environments. Our proposed reaction mechanisms assume either radiative stabilization or sufficiently long-lived intermediates to permit subsequent reactions under the ultra-low-pressure conditions of the ISM. While we apply RRKM/master equation (RRKM/ME) simulations where applicable, significant uncertainties remain, especially in estimating radiative association rates and the stability of reactive complexes. These assumptions are now explicitly discussed, and we underscore the need for both experimental and theoretical investigations into the dynamics and energetics of such stabilization pathways. We note that ion–molecule and grain-surface reactions are known to play crucial roles in the formation of prebiotic species such as formamide. While our current work focuses exclusively on neutral–neutral gas-phase pathways, we now recognize the potential importance of alternative routes, including ion-mediated reactions and grain-surface processes such as HCO + NH₂ coupling on dust grains. Given the large permanent dipole moments predicted for both thioformamide and thiourea, we strongly recommend laboratory rotational spectroscopy studies to provide reference spectra for radioastronomical detection. Finally, we highlight the need for astrochemical modeling using our computed rate constants to evaluate the viability and relative importance of these proposed formation pathways across a range of ISM environments.