The GRETOBAPE gas-phase reaction network: the importance of being exothermic

TL;DR

The paper introduces GRETOBAPE, a revised gas-phase reaction network for astrochemistry, emphasizing the importance of verifying reaction exothermicity through quantum calculations to improve model accuracy.

Contribution

It presents a new astrochemical reaction network, GRETOBAPE, with all reactions checked for exothermicity using quantum mechanics, and provides a reduced version for specific astrochemical applications.

Findings

5% of reactions in the original network are endothermic

The Si-bearing species network is most affected by endothermicity filtering

The updated network and species database are publicly available

Abstract

The gas-phase reaction networks are the backbone of astrochemical models. However, due to their complexity and non-linear impact on the astrochemical modeling, they can be the first source of error in the simulations if incorrect reactions are present. Over time, following the increasing number of species detected, astrochemists have added new reactions, based on laboratory experiments and quantum mechanics (QM) computations as well as reactions inferred by chemical intuition and similarity principle. However, sometimes no verification of their feasibility in the interstellar conditions, namely their exothermicity, was performed. In this work, we present a new gas-phase reaction network, GRETOBAPE, based on the KIDA2014 network and updated with several reactions, cleaned from endothermic reactions not explicitly recognized as such. To this end, we characterized all the species in the GRETOBAPE network with accurate QM calculations. We found that 5% of the reactions in the original network are endothermic although most of them are reported as barrierless. The reaction network of Si-bearing species is the most impacted by the endothermicity cleaning process. We also produced a cleaned reduced network, GRETOBAPE-red, to be used to simulate astrochemical situations where only C-, O-, N- and S- bearing species with less than 6 atoms are needed. Finally, the new GRETOBAPE network, its reduced version, as well as the database with all the molecular properties are made publicly available. The species properties database can be used in the future to test the feasibility of possibly new reactions.

Tables

Table 1. Table 1: ZPE-corrected reaction energy difference, i.e. , Δ ​ ( Δ ​ H ​ ( 0 ) ) Δ Δ H 0 \Delta(\Delta\mathrm{H}(0)) , between previous literature and our data.

Reaction		$\Delta \left(\Delta H\left(0\right)\right)$	Their Level of Theory	Reference
Reactants	Products	[kJ/mol]
C2H5 + O	C2H4 + OH	1.26	CCSD(T)/aug-cc-pVTZ//B2PYLYP/aug-cc-pVTZ	Vazart et al. (2020)
C2H5 + O	CH3CHO + H	1.67
C2H5 + O	H2CO + CH3	1.23
CH + CH3OH	C2H4 + OH	0.36	CCSD(T)/aug-cc-pVTZ//B2PYLYP/aug-cc-pVTZ	Vazart et al. (2020)
CH + CH3OH	CH3CHO + H	0.77
CH + CH3OH	H2CO + CH3	0.34
SiH + S	SiS + H	4.29	CCSD(T)/aug-cc-pVTZ//B3LYP/aug-cc-pVTZ	Rosi et al. (2018)
SiH + S2	SiS + HS	-1.85
Si + SH	SiH + S	-3.92	MRCI/CBS	Mota et al. (2021)

Table 2. Table 2: List of the species removed from the original network after the cleaning process.

Species
SiC4H	SiC6H	SiC8H
SiC2CH3	SiC3H5	H2CSiCH
H2C6N+	H2C8N+	H2C10N+
PNH${}_{}{}^{+}{}_3{}^{}$	C4H5

Table 3. Table 3: Summary of the number of reactions in the original and final GRETOBAPE networks.

	Neutral-Neutral	Ion-Neutral	Cation-Anion	Recombination	Unimolecular	Total
Pre-cleaning	1007	3596	1401	640	596	7240
Endothermic	58	140	108	-	-	306
Domino effect	2	-	-	9	12	23
Post-cleaning	947	3456	1293	631	584	6911

Table 4. Table 4: Summary of the 488 species contained in the cleaned network GRETOBAPE , together with the elements composing them.

Element	Number of species containing n times the element											Total
	1	2	3	4	5	6	7	8	9	10	11
H	123	81	53	38	21	12	7	1	0	0	0	336
He	3	0	0	0	0	0	0	0	0	0	0	336
C	100	66	57	34	21	19	20	17	19	12	2	367
N	115	8	0	0	0	0	0	0	0	0	0	123
O	93	24	0	0	0	0	0	0	0	0	0	117
F	7	0	0	0	0	0	0	0	0	0	0	7
Na	6	0	0	0	0	0	0	0	0	0	0	6
Mg	3	0	0	0	0	0	0	0	0	0	0	3
Si	53	0	0	0	0	0	0	0	0	0	0	53
S	38	7	0	0	0	0	0	0	0	0	0	45
P	38	0	0	0	0	0	0	0	0	0	0	38
Cl	10	0	0	0	0	0	0	0	0	0	0	10
Fe	2	0	0	0	0	0	0	0	0	0	0	2

Table 5. Table 5: List of neutral species with missing reactions with the most abundant cations in the cleaned newtork GRETOBAPE .

Species	H+	He+	H${}_{}{}^{+}{}_3{}^{}$	H3O+	HCO+	Species	H+	He+	H${}_{}{}^{+}{}_3{}^{}$	H3O+	HCO+	Species	H+	He+	H${}_{}{}^{+}{}_3{}^{}$	H3O+	HCO+
NH				$\star$		HPO				$\star$		CH2CHCN				$\star$
OH				$\star$		OCS				$\star$		C6H				$\star$
HF	$\star$			$\star$	$\star$	SiO2	$\star$			$\star$		HC5N				$\star$
NaH				$\star$		SO2	$\star$			$\star$	$\star$	C7				$\star$
C2				$\star$		CH3				$\star$		C6N	$\star$	$\star$		$\star$
MgH				$\star$		H2CN	$\star$	$\star$	$\star$	$\star$	$\star$	C2H6				$\star$
CN					$\star$	SiH3				$\star$	$\star$	C3H5	$\star$	$\star$	$\star$	$\star$	$\star$
CO	$\star$			$\star$	$\star$	HOOH	$\star$	$\star$	$\star$	$\star$	$\star$	CH3OCH2	$\star$		$\star$
N2	$\star$			$\star$		SiCH2				$\star$		CH2CHC2H	$\star$	$\star$	$\star$	$\star$	$\star$
NO				$\star$	$\star$	H2CS				$\star$		CH3C3N				$\star$
PH				$\star$		C3N	$\star$					C7H				$\star$
O2				$\star$		l-SiC3				$\star$		HC6N	$\star$		$\star$	$\star$	$\star$
HS				$\star$		C3P	$\star$					C8				$\star$
HCl				$\star$	$\star$	CH4				$\star$		C7N	$\star$			$\star$
SiN				$\star$	$\star$	CH2NH				$\star$		C7O	$\star$	$\star$		$\star$
SiO				$\bullet$	$\bullet$	CH3O	$\star$			$\star$		CH3CHCH2	$\star$	$\star$	$\star$	$\star$	$\star$
NS				$\star$		SiH4				$\star$		CH3C4H				$\star$
PO				$\star$		H2CCN				$\star$	$\star$	C8H				$\star$
CCl	$\star$		$\star$	$\star$	$\star$	NH2CN				$\star$		HC7N				$\star$
ClO	$\star$		$\star$	$\star$	$\star$	CH2PH				$\bullet$	$\bullet$	C9				$\star$
SO				$\star$		HCOOH				$\star$		C8N	$\star$	$\star$		$\star$
SiS				$\star$		c-HCCHSi				$\star$		C3H7	$\star$	$\star$	$\star$	$\star$	$\star$
S2				$\bullet$		C5				$\star$		CH2CHCHCH2	$\star$	$\star$	$\star$	$\star$	$\star$
H2O				$\star$		C4N	$\star$	$\star$		$\star$		CH3C5N				$\star$
HCO				$\star$		SiC3H				$\star$		C9H				$\star$
HNO				$\star$		SiC4				$\star$		HC8N	$\star$		$\star$	$\star$	$\star$
PH2				$\star$		HCCCHO	$\star$	$\star$				C10				$\star$
O2H	$\star$	$\star$	$\star$	$\star$	$\star$	c-C3H2O	$\star$	$\star$				C9N	$\star$			$\star$
CCN				$\star$		H2CCCO	$\star$	$\star$		$\bullet$		C9O	$\star$	$\star$		$\star$
NaOH				$\star$		C5H				$\star$		C3H8	$\star$	$\star$	$\star$	$\star$	$\star$
CCO				$\bullet$		HC4N	$\star$		$\star$	$\star$	$\star$	CH3C6H				$\star$
OCN	$\star$		$\star$	$\star$	$\star$	C6				$\star$		C10H				$\star$
HNSi				$\bullet$		C5N	$\star$			$\star$		HC9N				$\star$
CO2				$\star$		C5O	$\star$	$\star$		$\star$		C11	$\star$		$\star$	$\star$	$\star$
N2O				$\star$	$\bullet$	C2H5				$\star$		C10N	$\star$	$\star$		$\star$
HCS				$\star$	$\star$	CH3NH2				$\star$		c-C6H6	$\star$
NO2	$\star$	$\star$		$\star$	$\star$	CH2CCH2	$\star$	$\star$	$\star$	$\star$	$\star$	CH3C7N				$\star$

Table 6. Table 6: Species with reported observations in the ISM (McGuire, 2022 ) which are not present in the final GRETOBAPE network.

Carbon, nitrogen and oxygen-bearing species
4-atoms
HCCO	HONO	HNCN	CNCN	HC2N
5-atoms
HNCNH	NH2OH	H2NCO+	NCCNH+
6-atoms
HNCHCN	CH3NC	C5N-
7-atoms
c-C2H4O	c-C3HCCH	CH2CHOH	HC4NC	CH3NCO	HOCH2CN	HC5O
8-atoms
CH2CHCHO	CH3CHNH	NH2CH2CN	NH2CONH2
9-atoms
HC7O	H2CCCHCCH	HCCCHCHCN	H2CCHC3N	CH3NHCHO	CH3CONH2
10-atoms
HOCH2CH2OH	CH3CH2CHO	CH3CHCH2O	CH3OCH2OH
11-atoms
C2H5OCHO	CH3COOCH3	CH3COCH2OH	NH2CH2CH2OH
12-atoms
n-C3H7CN	i-C3H7CN	1-C5H5CN	2-C5H5CN
$>$13-atoms
HC11N	c-C6H5CN	C9H8	1-C10H7CN	2-C10H7CN
Sulfur-bearing species
HSC	NCS	HNCS	HCCS	CHOSH	HCSCN	HCSCCH
H2CCS	H2CCCS	CH3SH	C5S	CH3CH2SH
Silicon-bearing species
SiCN	SiCSi	SiC3	SiH3CN	CH3SiH3

Table 7. Table 7: Results of the modeling of the cold molecular cloud (§  4.1.1 ). For each species, listed in the first column, columns 2-5 list the minimum ( min ) and maximum ( max ) values of the network post- ( GRETOBAPE ) over pre-cleaning predicted abundance ratios in the time range 1 × 10 4 absent superscript 10 4 \times 10^{4} –1 × 10 5 absent superscript 10 5 \times 10^{5} and 1 × 10 5 absent superscript 10 5 \times 10^{5} –1 × 10 6 absent superscript 10 6 \times 10^{6} yr, respectively. Column 6-7 report the ratio and the abundance at steady-state (SS). The last two columns list the species largest reached abundance and at what time. The species are ordered following the maximum ratio variation during the entire simulation. Finally, boldfaced species have been detected in the ISM.

Interval	1$\times {10}^4$–1$\times {10}^5$ yr		1$\times {10}^5$–1$\times {10}^6$ yr		SS		max(Abd/H)
Species	min(ratio)	max(ratio)	min(ratio)	max(ratio)	ratio	Abd/H	Abd/H	time[yr]
HCSi	4.2$\times {10}^{-2}$	6.2$\times {10}^{-1}$	6.7$\times {10}^{-2}$	8.6$\times {10}^{-1}$	5.2$\times {10}^{-2}$	1.1$\times {10}^{-14}$	1.1$\times {10}^{-12}$	7.5$\times {10}^4$
SiO+	3.8$\times {10}^2$	1.6$\times {10}^3$	1.3$\times {10}^3$	2.1$\times {10}^3$	2.7$\times {10}^3$	8.2$\times {10}^{-12}$	8.2$\times {10}^{-12}$	2.4$\times {10}^6$
CHSi+	3.1$\times {10}^3$	1.0$\times {10}^5$	4.7$\times {10}^4$	2.1$\times {10}^5$	9.2$\times {10}^4$	1.0$\times {10}^{-13}$	2.7$\times {10}^{-12}$	6.1$\times {10}^4$
HCCSi	1.2$\times {10}^{-1}$	5.8$\times {10}^{-1}$	8.1$\times {10}^{-3}$	1.1$\times {10}^{-1}$	2.9$\times {10}^{-1}$	5.2$\times {10}^{-13}$	6.5$\times {10}^{-12}$	1.3$\times {10}^5$
SiH	1.2	6.7	4.5	1.3$\times {10}^1$	9.5	2.1$\times {10}^{-11}$	2.2$\times {10}^{-11}$	2.6$\times {10}^3$
HNSi	8.4$\times {10}^{-1}$	2.0	2.2$\times {10}^{-1}$	3.4	2.2$\times {10}^{-2}$	1.4$\times {10}^{-11}$	1.3$\times {10}^{-10}$	7.1$\times {10}^5$
HSiNH+	8.6$\times {10}^{-1}$	2.5	1.8$\times {10}^{-1}$	3.2	1.8$\times {10}^{-2}$	2.5$\times {10}^{-14}$	2.9$\times {10}^{-13}$	6.4$\times {10}^5$
Si+	1.4	1.0$\times {10}^1$	4.5	2.4$\times {10}^1$	9.4	1.8$\times {10}^{-9}$	7.8$\times {10}^{-09}$	1.0$\times {10}^2$
SiF+	1.4	1.0$\times {10}^1$	4.3	2.3$\times {10}^1$	8.9	6.6$\times {10}^{-13}$	8.0$\times {10}^{-13}$	9.6$\times {10}^5$
SiH${}_{}{}^{+}{}_2{}^{}$	1.4	1.0$\times {10}^1$	4.4	2.3$\times {10}^1$	9.2	3.7$\times {10}^{-12}$	4.7$\times {10}^{-12}$	8.7$\times {10}^5$
H2CS	1.8$\times {10}^{-2}$	3.7$\times {10}^{-2}$	2.8$\times {10}^{-2}$	6.7$\times {10}^{-2}$	1.2$\times {10}^{-1}$	3.4$\times {10}^{-12}$	2.9$\times {10}^{-11}$	1.7$\times {10}^5$
SiC	1.9	5.4	6.8$\times {10}^{-1}$	5.1	7.1	6.0$\times {10}^{-13}$	2.5$\times {10}^{-11}$	2.9$\times {10}^3$
SiH+	1.3	3.1	2.0	4.9	7.7	1.9$\times {10}^{-12}$	2.3$\times {10}^{-12}$	8.7$\times {10}^5$
HSiO+	1.5$\times {10}^{-1}$	1.9$\times {10}^{-1}$	1.3$\times {10}^{-1}$	1.9$\times {10}^{-1}$	1.5$\times {10}^{-1}$	7.4$\times {10}^{-12}$	7.4$\times {10}^{-12}$	1.0$\times {10}^7$
SiN	1.3	5.4	1.8	5.8	3.1	2.5$\times {10}^{-12}$	1.8$\times {10}^{-11}$	6.4$\times {10}^5$
SiNC	8.8$\times {10}^{-1}$	2.8	1.6	5.8	3.4	2.5$\times {10}^{-15}$	1.2$\times {10}^{-13}$	6.8$\times {10}^4$
l-SiC3	1.1	5.5	3.6	7.7	6.5	1.0$\times {10}^{-14}$	9.0$\times {10}^{-12}$	3.5$\times {10}^5$
SiC3H	5.8$\times {10}^{-1}$	8.8$\times {10}^{-1}$	5.8$\times {10}^{-1}$	1.1	4.5	1.8$\times {10}^{-15}$	1.5$\times {10}^{-11}$	2.1$\times {10}^5$
SiCH3	8.8$\times {10}^{-1}$	9.6$\times {10}^{-1}$	9.9$\times {10}^{-1}$	3.3	6.4	5.9$\times {10}^{-14}$	7.6$\times {10}^{-13}$	1.5$\times {10}^5$
l-C3H${}_{}{}^{+}{}_2{}^{}$	5.3$\times {10}^{-1}$	9.0$\times {10}^{-1}$	1.0$\times {10}^{-1}$	4.7$\times {10}^{-1}$	2.7$\times {10}^{-1}$	2.5$\times {10}^{-15}$	3.7$\times {10}^{-11}$	1.0$\times {10}^5$

Table 8. Table 8: Results of the modeling of the warm molecular outflow shock (§  4.1.2 ). For each species, listed in the first column, columns 2-5 list the minimum ( min ) and maximum ( max ) values of the network post- ( GRETOBAPE ) over pre-cleaning predicted abundance ratios in the time range 1 × 10 2 absent superscript 10 2 \times 10^{2} –1 × 10 3 absent superscript 10 3 \times 10^{3} and 1 × 10 3 absent superscript 10 3 \times 10^{3} –5 × 10 3 absent superscript 10 3 \times 10^{3} yr, respectively. Column 6-7 report the ratio and the abundance at 5 × 10 3 5 superscript 10 3 5\times 10^{3} yr. The last two column list the species largest reached abundance and at what time. The species are ordered following the value of the ratio at 5 × 10 3 absent superscript 10 3 \times 10^{3} yr, which is believed to approximately be the L1157-B1 age (§  4.1.2 ). Finally, boldfaced species have been detected in the ISM.

Interval	1$\times {10}^2$–1$\times {10}^3$ yr		1$\times {10}^3$–5$\times {10}^3$ yr		5$\times {10}^3$ yr		max(Abd/H)
Species	min(ratio)	max(ratio)	min(ratio)	max(ratio)	ratio	Abd/H	Abd/H	time[yr]
SiS+	3.9$\times {10}^7$	4.2$\times {10}^7$	4.2$\times {10}^7$	4.5$\times {10}^7$	4.2$\times {10}^7$	3.9$\times {10}^{-11}$	6.3$\times {10}^{-11}$	1.0$\times {10}^2$
CHSi+	1.5$\times {10}^6$	1.6$\times {10}^6$	1.4$\times {10}^6$	1.9$\times {10}^6$	1.4$\times {10}^6$	4.3$\times {10}^{-13}$	4.3$\times {10}^{-13}$	4.9$\times {10}^3$
SiO+	1.3$\times {10}^3$	1.5$\times {10}^3$	1.6$\times {10}^3$	7.9$\times {10}^3$	7.9$\times {10}^3$	1.1$\times {10}^{-11}$	1.1$\times {10}^{-11}$	5$\times {10}^3$
SiN	1.2$\times {10}^{-2}$	1.9$\times {10}^{-2}$	2.8$\times {10}^{-3}$	7.7$\times {10}^{-3}$	5.5$\times {10}^{-3}$	3.6$\times {10}^{-13}$	3.6$\times {10}^{-13}$	5$\times {10}^3$
HSiO+	1.1$\times {10}^{-1}$	1.1$\times {10}^{-1}$	8.3$\times {10}^{-2}$	1.0$\times {10}^{-1}$	8.2$\times {10}^{-2}$	3.5$\times {10}^{-10}$	3.5$\times {10}^{-10}$	3.9$\times {10}^3$
CH3CHCH2	5.0$\times {10}^{-2}$	5.2$\times {10}^{-2}$	5.3$\times {10}^{-2}$	6.5$\times {10}^{-2}$	6.5$\times {10}^{-2}$	3.6$\times {10}^{-13}$	4.2$\times {10}^{-13}$	1.0$\times {10}^2$
SiH${}_{}{}^{+}{}_2{}^{}$	4.0$\times {10}^2$	4.3$\times {10}^2$	4.4$\times {10}^2$	5.0$\times {10}^2$	4.9$\times {10}^2$	1.5$\times {10}^{-12}$	1.5$\times {10}^{-12}$	5$\times {10}^3$
Si+	4.0$\times {10}^2$	4.3$\times {10}^2$	4.4$\times {10}^2$	4.9$\times {10}^2$	4.8$\times {10}^2$	1.5$\times {10}^{-10}$	1.6$\times {10}^{-10}$	3.9$\times {10}^3$
SiH	7.1$\times {10}^1$	1.4$\times {10}^2$	1.6$\times {10}^2$	3.7$\times {10}^2$	3.7$\times {10}^2$	3.9$\times {10}^{-12}$	3.1$\times {10}^{-12}$	4.3$\times {10}^3$
HNSi	1.6$\times {10}^{-2}$	1.9$\times {10}^{-2}$	1.1$\times {10}^{-2}$	1.5$\times {10}^{-2}$	1.1$\times {10}^{-2}$	1.3$\times {10}^{-11}$	1.4$\times {10}^{-11}$	8.0$\times {10}^2$
SiNC+	5.5$\times {10}^1$	7.7$\times {10}^1$	8.4$\times {10}^1$	1.2$\times {10}^2$	1.1$\times {10}^2$	2.5$\times {10}^{-12}$	2.4$\times {10}^{-12}$	4.7$\times {10}^3$
HSS	9.6$\times {10}^{-1}$	1.0	2.3$\times {10}^{-1}$	9.2$\times {10}^{-1}$	2.3$\times {10}^{-1}$	4.7$\times {10}^{-13}$	4.7$\times {10}^{-13}$	5$\times {10}^3$
c-HCCHSi	4.3$\times {10}^{-1}$	1.2	2.3$\times {10}^{-1}$	3.7$\times {10}^{-1}$	2.3$\times {10}^{-1}$	6.2$\times {10}^{-12}$	6.1$\times {10}^{-12}$	4.9$\times {10}^3$
c-SiC2	4.3$\times {10}^{-2}$	5.5$\times {10}^{-2}$	6.3$\times {10}^{-2}$	2.1$\times {10}^{-1}$	2.1$\times {10}^{-1}$	7.2$\times {10}^{-13}$	7.3$\times {10}^{-13}$	4.7$\times {10}^3$
SiS	1.7$\times {10}^{-1}$	2.8$\times {10}^{-1}$	1.5$\times {10}^{-1}$	1.7$\times {10}^{-1}$	1.5$\times {10}^{-1}$	8.4$\times {10}^{-12}$	8.5$\times {10}^{-12}$	4.9$\times {10}^3$
SO${}_{}{}^{+}{}_2{}^{}$	1.3$\times {10}^1$	1.4$\times {10}^1$	1.4$\times {10}^1$	1.5$\times {10}^1$	1.5$\times {10}^1$	1.4$\times {10}^{-13}$	1.4$\times {10}^{-13}$	5$\times {10}^3$
HCCSi	2.4$\times {10}^{-2}$	3.8$\times {10}^{-2}$	4.4$\times {10}^{-2}$	1.3$\times {10}^{-1}$	1.3$\times {10}^{-1}$	4.8$\times {10}^{-12}$	4.7$\times {10}^{-12}$	4.9$\times {10}^3$
Si	1.6$\times {10}^{-1}$	1.7$\times {10}^{-1}$	1.2$\times {10}^{-1}$	1.6$\times {10}^{-1}$	1.2$\times {10}^{-1}$	3.3$\times {10}^{-10}$	2.1$\times {10}^{-09}$	1.0$\times {10}^2$
HSSH	8.0$\times {10}^{-1}$	1.0	1.0$\times {10}^{-1}$	6.7$\times {10}^{-1}$	1.0$\times {10}^{-1}$	1.6$\times {10}^{-13}$	1.6$\times {10}^{-13}$	5$\times {10}^3$
SiH+	2.4	3.3	3.7	5.7	5.1	5.3$\times {10}^{-13}$	5.4$\times {10}^{-13}$	4.9$\times {10}^3$
SiC	1.2	1.2	1.2	3.0	3.0	2.3$\times {10}^{-11}$	2.3$\times {10}^{-11}$	5$\times {10}^3$
C3O	9.8	2.1$\times {10}^1$	2.0	1.6$\times {10}^1$	2.0	2.4$\times {10}^{-12}$	7.9$\times {10}^{-12}$	1.0$\times {10}^2$

Table 9. Table 9: Summary of the files containing the reaction networks and species chemical data provided in the Supporting Material and a brief description of their content.

File Name	Description
ReadMe.txt	File with the full description of all the files in the Supporting Material.
GRETOBAPE-pre.dat	Original reaction network, before the cleaning process, described in § 3.1.
GRETOBAPE.dat	New reaction network, after the cleaning process, described in § 3.2.
GRETOBAPE-red.dat	Reduced network from GRETOBAPE.dat, described in § 3.3.
GRETOBAPE-endo.dat	List of reactions deleted because of the endothermicity criteria, described in § 2.3.1.
GRETOBAPE-endo-0-10kJmol.dat	List of reactions with an endothermicity lower than 10 kJ/mol: they are present in GRETOBAPE.dat.
GRETOBAPE-domino.dat	List of reactions deleted because of the domino effect, described in § 2.3.2.
GRETOBAPE-domino-endo-Si.dat	List of reactions involving Si-bearing species and deleted because of the endothermicity and domino
	effect criteria, discussed in § 3.2.2.
Database_Molecules.csv	Extracted and organized chemical information of all species in GRETOBAPE.dat, described in 2.1.
species.zip	Zip file containing the .XYZ files of all species studied in this work, described in 2.1.

Equations

$$k(\mathrm{T})=\alpha {\left(\frac{\mathrm{T}}{300\mathrm{K}}\right)}^{\beta }\exp \left(-\frac{\gamma }{\mathrm{T}}\right)$$

Full text

astrochemistry — ISM: Reaction network; astrochemistry; cold molecular cloud; molecules

††journal: ApJS††software: NetworkX (Hagberg et al., 2008), Gaussian16 (Frisch et al., 2016), krome (Grassi et al., 2014).

1 Introduction

Since the discovery of the first molecules in the interstellar medium (ISM) in the visible spectrum (Swings & Rosenfeld, 1937; Dunham Jr, 1937; McKellar, 1940), the astrochemical community has been trying to understand their origin and how they can survive the harsh interstellar conditions. With the evolution of astronomical facilities and, more specifically, with the advent of radio telescopes, new and more complex molecules have been detected in various astronomical objects with an almost constant discovery rate of $\sim 4-6$ new molecules per year (McGuire, 2022). To date, more than 270 molecules, composed of 19 different elements, have been detected in the interstellar and circumstellar medium, and this number is steadily increasing. The list of detected interstellar species can be found in the Cologne Database for Molecular Spectroscopy (CDMS) (Endres et al., 2016)111https://cdms.astro.uni-koeln.de/classic/molecules, in the McGuire (2022) census (and the related Python package222https://github.com/bmcguir2/astromol), and in The Astrochymist website333http://www.astrochymist.org/astrochymist_ism.html.

Understanding how interstellar molecules are formed, destroyed and connected, in addition to be interesting per se, can help inferring the properties, history, and evolution of the astronomical objects in which they have been detected. Thus, in parallel with the identification of molecules and the measurements of their abundances in different astronomical objects, soon a community dedicated to developing astrochemical models began to play an important role (Herbst & Klemperer, 1973; Prasad & Huntress Jr, 1980). A crucial element of astrochemical models is the gas-phase reaction network, namely the lists of the gas-phase chemical reactions with their rate constants, that describe how efficient the reactions are as a function of the temperature, and product branching ratios.

The increasing number of detected species has prompted modelers and chemists to add new reactions to the existing networks in order to reproduce the observed species and their abundances. In this way, the number of reactions, reactants and products have been increasing each time a new molecule was detected, which also led to an increase in the network size and complexity. Nowadays, the publicly available astrochemical gas-phase networks, KIDA444http://kida.obs.u-bordeaux1.fr/ (Wakelam et al., 2015) and UMIST555http://udfa.ajmarkwick.net/ (McElroy et al., 2013), contains about 8000 reactions involving more than 500 species.

Unfortunately, the vast majority ($\geq 80$%) of the rates and branching ratios of the reported gas-phase reactions has not been measured or computed and is often based on approximate estimates (e.g. using the Capture Theory: Su & Chesnavich, 1982; Herbst, 2006; Woon & Herbst, 2009; Loison et al., 2013), educated guesses on similarity principles or simple chemical intuition. In addition, even when some experimental measurements are available, the estimated rate constants may have substantial uncertainties as they are often based on experiments at room temperature (e.g. Anicich, 2003). The uncertainty on the reaction rates is a well-known problem in the astrochemistry community. In parallel with the efforts of studying reactions via experiments and theoretical computations, sensitivity analysis studies have provided the impact of these uncertainties on the abundances predicted by astrochemical models (e.g. Wakelam et al., 2006; Penteado et al., 2017).

One possible and mandatory first step to improve the reliability of the gas-phase reaction networks, that has never been performed systematically on the astrochemical networks is to estimate the exo/endothermicity of the gas-phase reactions, using accurate physico-chemical quantum mechanical (QM) data for each species. Indeed, given the ISM conditions (temperatures usually less than $\sim$100-200 K), strongly endothermic reactions are inhibited and, therefore, can be excluded from the astrochemical networks without complex and time-consuming reaction transition states studies or experiments (Smith, 2006, 2011). This verification is now possible, thanks to the availability of two theoretical studies that together provide reliable QM data (electronic state, electron spin multiplicity, geometry, harmonic frequencies, absolute electronic energy, and dipole moment) for the totality of the species in the KIDA network. The first one by Woon & Herbst (2009) computed the QM data for a large number of neutral molecules involved in the KIDA reaction network. The second study by Tinacci et al. (2021) complemented the first one and provided QM data for all the ions and the remaining neutral species of the KIDA network. Please note that the two studies have been carried out at the same electronic level of theory, so that they can be reliably used for computing the exo/endothermicity of all the reactions. Also, Tinacci et al. (2021) created a database666https://aco-itn.oapd.inaf.it/aco-public-datasets/theoretical-chemistry-calculations with the above QM data for all the cations.

The KIDA and UMIST networks contain such a large number of reactions because their goal is to reproduce the observed abundances of large molecules, which have usually low (between $10^{-8}$ and $10^{-12}$) abundances. However, in some cases, such as astrochemical models coupled with hydrodynamical simulations of cloud collapse or protoplanetary disks formation and evolution, the goal is to reliably reproduce only the most abundant species. In these cases, the use of large reactions networks is not only useless but also detrimental, because it demands spareable computing time. Indeed, several authors have proposed reaction network reduction techniques in the last few years, each using different criteria and focusing on specific goals (Oppenheimer & Dalgarno, 1974; Nelson & Langer, 1999; Glover et al., 2010; Grassi et al., 2013, 2021).

The goal of the present article is twofold. The first goal is to produce a gas-phase reaction network ”cleaned” from endothermic reactions not reported as such. To this purpose, we perform a systematic study of the endo/exo thermicity of all reactions present in a new network, GRETOBAPE, which is based on the KIDA2014 one and updated through the years by our group. In the process, we also verify that each species present in the network is not a sink and that the obvious reactions with the most abundant molecular cations are present. The second goal is to produce a reduced network from the cleaned one, which only includes H-, C-, N-, O-, and S-bearing species with not more than six atoms, to be employed when the use of a complete network is unnecessary.

The article is organized as follows. Section 2 details the adopted tools to deal with the reaction network, the cleaning procedure, and the computational methodology used to characterize all the chemical species present in the network. In Section 3, we describe the new reaction networks, GRETOBAPE and GRETOBAPE-red, respectively resulting from the cleaning and reduction processes, as well as a list of the obvious reactions possibly missing in the network. In Section 4, we discuss the obtained results and their impact on astrochemical models of cold molecular clouds and warm molecular outflow shocks simulations. Finally, in Section 5, we present our conclusions, and in Section 6 we provide the hyperlinks to the publicly available online databases and reaction networks produced in this work.

2 Methodology

2.1 Computational details for QM calculation

All the species present in the astrochemical reaction network were characterized via QM calculations to obtain their optimized structure, absolute electronic energy, and Zero Point Energy (ZPE) correction. Most of the neutral species structures and their corresponding electronic spin multiplicity ($m_{s}$) in the ground state, a fundamental piece of information for the electronic QM calculations, are taken from Woon & Herbst (2009) while for the cations, the information is taken from our previous work (Tinacci et al., 2021). For the chemical species that have not been theoretically characterized yet, we inferred here the structure and $m_{s}$ using the same approach introduced in Tinacci et al. (2021). The list of those species with their chemical data are reported in the Supporting Material, in the file species.zip.

In order to have all species characterized at the same level of theory, i.e., so to have a consistent method, we used the DFT M06-2X (Zhao & Truhlar, 2008) coupled with the triple-$\zeta$ Dunning’s correlation consistent basis set (cc-pVTZ) (Kendall Jr, 1992; Woon & Dunning Jr, 1993) for geometry optimization and harmonic frequency calculations. The electronic absolute energy was refined at a CCSD(T) (Knowles et al., 1993, for closed shell species) or RO-CCSD(T) (Watts et al., 1993, for open shell species) level of theory, in conjunction with an augmented triple-$\zeta$ correlation consistent basis set (aug-cc-pVTZ; Kendall Jr, 1992). All the molecular structures and properties computed with the adopted methodology were benchmarked in our previous work (Tinacci et al., 2021). All calculations were carried out with the Gaussian16 program (Frisch et al., 2016) and we kept the default values set up in the program.

Due to the importance to have accurate reaction energies for our reactions network cleaning methodology, we report a benchmark of our resulting ZPE-corrected reaction energy with respect to other published studies in Table 1. All differences are smaller than 4 kJ/mol, which corroborates the accuracy of present results. Please note that the differences are mostly due to a different treatment of the spin contamination, which we improved in the present study.

We decided not to compute the reaction energies using the standard enthalpy of formation ($\Delta$H${}_{f}^{\circ}$(T), usually given at 298.15 K or 0 K), but instead we used the absolute electronic energy. This choice minimizes the propagation error in the QM computations, which would be larger if carried out using the formation enthalpy. One way to overcome this limitation would be to compute the $\Delta$H${}_{f}^{\circ}$(T) by the Isodesmic and Homodesmic cycles (Cramer, 2013; Wheeler et al., 2009), which compensates the errors. The problem with this technique is the difficulty of automatize it. In the past, some attempts were carried out but they cannot be generalized to all chemical species (e.g. Cavallotti, 2022).

It is worth noticing that many astrochemical species, radicals and exotic species with respect to the terrestrial chemistry, lack laboratory measured $\Delta$H${}_{f}^{\circ}$(T). Therefore, limiting the control on the exo/endo thermicity of the reaction to only those having laboratory experimental data would greatly reduce our capacity to ”clean” the astrochemical networks.

Overall, we characterized 542 species. The whole list is publicly available at the ACO (Astro-Chemical Origin) project website777https://aco-itn.oapd.inaf.it/home. In addition, a data frame (in .csv format) with all the properties of the studied species is available in the Supplementary Material of this article (see Table 9.

2.2 Graph theory for reaction networks

In the astrochemical community, the reaction network is usually encoded and manipulated with a matrix approach. In this work, we adopt an approach based on the graph theory, so that the network manipulation can be easily obtained via the NetworkX (Hagberg et al., 2008) Python package (e.g., node connectivity, centrality, sinks, and other characteristics). Instead of using a multi-di-graph approach as done by Grassi et al. (2013), the mathematical object we used to describe the reaction network is the colored directed graph (also known as colored di-graph). Here, each node has a color attribute and the connections are directed (Barabási, 2013). A colored node can have different attributes depending on the color: in our case, it can be a species-node or a rate constant-node. In addition, in a directed graph, the edges have a direction, i.e. the adjacency (or connectivity) matrix is not symmetrical, so that the connection between nodes is orientated (e.g. Barabási, 2013; Newman, 2018). In Figure 1, we report the codification of a generic reaction network given in the above-mentioned three ways: the simple reaction list, the reaction network visualization and, finally, the corresponding adjacency matrix of the colored directed graph. The Python scripts used to encode and characterize the network properties are publicly available in a GitHub repository888https://github.com/TinacciL/GreToBaPe_Cleaning.

2.3 Network cleaning procedure

2.3.1 Identification of endothermic reactions

As explained in the Introduction, the first goal of the present work is to produce a network ”cleaned” from endothermic reactions that have not been reported as such (see below), since they cannot occur in the cold ISM. To this end, we computed the reaction enthalpy $\Delta\mathrm{H}$ of the neutral-neutral, ion-neutral, and ion-ion reactions present in the reaction network. Note that we did not consider photo- and electron- induced reactions (involving electrons, cosmic-ray particles and UV radiation) because they entail strongly energetic processes, whose products are usually not in their fundamental electronic state.

In the ISM, molecules primarily reside in their ground electronic state and lowest vibrational level, with only the rotational levels populated, depending on the gas conditions. Since the rotational energy has a negligible contribution to the total reaction energy, we used the electronic energy and the vibrational energy at 0 K (also called Zero-point Energy correction, i.e., ZPE) to evaluate the reaction enthalpy $\Delta\mathrm{H}(0)$ of each reaction present in the reaction network. Then, an endothermic reaction has $\Delta\mathrm{H}(0)$ larger than zero.

In addition, in the astrochemical networks the rate constants of the neutral-neutral and ion-neutral (but not ion-polar) reactions are reported following the modified Arrhenius equation (Kooij, 1893; Laidler, 1996):

[TABLE]

where T is the gas temperature in K, and $\alpha$, $\beta$, and $\gamma$ are often derived as fitted parameters of either experimental or theoretical values of $k(\mathrm{T})$. It should be emphasized that $\gamma$ does not always have an actual physical meaning. In other words, $\gamma$ is often not a real activation barrier. That said, it can be safely assumed that, if the reaction has a barrier, it is lower or equal to $\gamma$. Therefore, in order to consider that a reaction has been incorrectly inserted in the reaction network, we apply the additional criterion that $\Delta\mathrm{H}(0)>\gamma~{}R$, where $\gamma$ is multiplied by the constant $R$ to have it in kJ/mol.

Finally, our computations have an intrinsic uncertainty, i.e. the accuracy of the electronic structure QM calculations, which we conservatively evaluate to be 10 kJ/mol. Therefore, when retaining or excluding a reaction from the network we consider a 10 kJ/mol threshold regardless the actual meaning of $\gamma$ either a real activation barrier or the results of an analytical fit.

In summary, a reaction is considered correctly included in the network if it satisfies the following criteria:

if the reaction is not encoded with the modified Arrhenius equation (i.e. the ion-pol formula for barrierless reactions), $\Delta\mathrm{H}(0)\leq 10$ kJ/mol;

• 2)

if the reaction is encoded with the modified Arrhenius equation, $\Delta\mathrm{H}(0)\leq\gamma~{}R+10$ kJ/mol.

Reactions that do not satisfy the two above criteria are excluded.

2.3.2 The domino effect of removing endothermic reactions

In order for the network to be meaningful, each species must have at least one destruction route (i.e. loss) and one formation route (i.e. production).

Removing reactions from the network can affect other species if they are not formed or destroyed through other reactions than those removed. In other words, if a species becomes a source (species with no formation routes) or a sink (species with no destruction routes) it is also removed from the network. As a consequence, the reactions involving the removed species are also removed.

Therefore, the removal of one or more initial reactions triggers a domino effect, causing the removal of other reactions and species. Our graph theory approach, described in § 2.2, allows us to easily follow the domino effect and identify the reactions and species to be removed.

2.4 Network reduction procedure

As stated in the Introduction, the second goal of this work is to produce a reduced network from the cleaned one, which only includes H-, C-, N-, O- and S- bearing species with at most six atoms, to be used when the larger network is unnecessary. To this end, we used the following criteria and steps:

We only consider C-, N-, O- and S-bearing species, plus Fe and Mg atoms, because their ions can be important positive charge carriers. 2. 2.

We only consider molecules with six or less atoms (either reactant or product). 3. 3.

The removal of larger molecules may cause the domino effect, explained in § 2.3.2, which can lead to the removal of other reactions and species also with less than six atoms. 4. 4.

Methanol is the only exception to rules 2 and 3. In the gas phase, it is formed by reactions involving molecules with more than 6 atoms, so it should be removed, based on rule 3. However, methanol is an important component of the interstellar grain mantles, a major carbon carrier, believed to be mainly formed on the grain surfaces and released into the gas phase by thermal and non-thermal effects. In practice, if the model contains surface reactions that form methanol, it could be injected into the gas-phase by e.g. the sublimation or sputtering of the grain mantles. Once in the gas-phase, methanol can react with species having less than 6 atoms (e.g. HCO+, H${}_{3}^{+}$, OH etc..). For this reason, we keep methanol in the reduced network. If no methanol is produced by surface chemistry, then its abundance would be zero because no gas-phase reaction would form it so that it is not a source nor a sink. 5. 5.

We remove all anions. We had two reasons for that. The first one is that, usually, they are not crucial species in studies that necessitate a reduced network. The second one is that there is still a debate on their rate of formation (e.g. Lara-Moreno et al., 2019b, a).

3 Results

3.1 Original astrochemical network

The reaction network from which we start is based on the KIDA 2014 network (Wakelam et al., 2015), updated following studies from our and other groups, as follows: Hoyermann et al. (1981); DeFrees et al. (1985); Meot-Ner et al. (1986); Dóbé et al. (1991); Anicich (2003); Hamberg, M. et al. (2010); Cheikh (2012); Lawson et al. (2012); Fournier (2014); Loison et al. (2014); Balucani et al. (2015); Neufeld et al. (2015); Vazart et al. (2016); Loison et al. (2016); Urso et al. (2016); Codella et al. (2017); Skouteris et al. (2017); Fontani et al. (2017); Balucani et al. (2018); Gao et al. (2018); Rosi et al. (2018); Skouteris et al. (2018); Sleiman et al. (2018); Ascenzi et al. (2019); Ayouz et al. (2019); Lara-Moreno et al. (2019a); Skouteris et al. (2019); Urso et al. (2019); Balucani (2020); Codella et al. (2020); Vazart et al. (2020); Blázquez et al. (2020); Mancini & Aragão (2021), and data, not yet published, courteously provided by Luca Mancini on some P-bearing reactions. We also add the recombination reactions of F+ and P+ with free electrons, which are not present in the KIDA astrochemical reaction network, although the second one it is in the UMIST database. These data are taken from the database of Badnell (2006) and given in the usual Arrhenius-Kooij format (Eq. 1).

In summary, the pre-cleaned original network comprises 499 species and 7240 reactions. The list is reported in the Supporting Material (§ 7: file GRETOBAPE-pre.dat).

Note that only 5768 of the 7240 reactions have thermodynamic data. The remaining 1472 reactions do not either because they contain species for which the electronic structures are not available (236 reactions999Reactions involving the species Fe, Fe+, C2H${}_{7}^{+}$, C9H3N+ or bimolecular reactions having as a product an electron or photon.) or they are recombination (640 reactions), or cosmic-ray particles and UV radiation (596 reactions).

3.2 New cleaned network: GRETOBAPE

In this section, we introduce the new reaction network obtained applying the criteria of § 2.3 to the network just described in the previous section. It is important to emphasize that we only tested the endothermicity of the original network, which is based on the products listed in KIDA and the works cited above. We did not try to correct reactions a posteriori comparing the products in other astrochemical networks, notably UMIST (but see § 3.2.3), where sometimes exothermic products of the reactions exist, for the reasons that will become clearer at the end of this section.

3.2.1 Overview of the removed species and reactions

Among the 5768 reactions that have thermodynamic data (see §  3.1) and for which we could carry out the endothermicity verification, 306 do not satisfy the two criteria described in § 2.3, namely $\sim 5$% of the studied reactions101010 Please note that the reaction H + O+ $\rightarrow$ H+ + O is endothermic by  17kJ/mol according our computations. However, the reaction has been studied by other authors (e.g. Stancil et al., 1999) and found to be exothermic. The reason for our erroneous result is that the ionization energy of O is a critical case that pushes the computational calculations to their limits and a difficult case for QM calculations due to the correlation errors even if we are using one of the best available methods.. Figure 2 shows the distribution of their endothermic energy. The majority of these reactions have endothermicities lower than about 150 kJ/mol while about 140 reactions even larger.

The list of the removed reactions based on the endothermicity criteria is reported in the Supporting Material (§ 7: file GRETOBAPE-endo.dat). In addition, we supply the list of 47 reactions which are endothermic but with a reaction enthalpy lower than 10 kJ/mol and for which more accurate calculations should be carried out before removing them for the network (§ 7: fileGRETOBAPE-endo-0-10kJmol.dat). Finally, because of the domino effect, described in § 2.3.2, 11 species and 23 additional reactions are removed, leading to a total of 329 removed reactions (§ 7: files GRETOBAPE-endo.dat plus GRETOBAPE-domino.dat).

The full list of deleted species is reported in Tab. 2. Table 3 provides an overview of the removed reactions, grouped according to each reaction class. The ion-neutral reactions are the most affected by the cleaning process in number (140), while the cation-anion reactions are the most affected in percentage ($7.7$%). Table 4 lists the number of species in the final cleaned network classified with their elements.

In summary, the new cleaned network GRETOBAPE contains 6911 reactions and 488 species.

3.2.2 Description of the removed reactions

The 329 removed reactions have the following properties:

• •

65 reactions are not encoded with the modified Arrhenius equation: they either do not satisfy the criterion 1 listed in § 2.3.1, namely $\Delta\mathrm{H}(0)\leq 10$ kJ/mol, or they are removed because of the domino effect.

• •

219 reactions are encoded via the Arrhenius-Kooji equation and have $\gamma=0$: they either do not satisfy the criterion 2 listed in § 2.3.1, namely $\Delta\mathrm{H}(0)\leq 10$ kJ/mol, or they are removed because of the domino effect.

• •

45 reactions are encoded via the Arrhenius-Kooji equation and have $\gamma\neq 0$: they either do not satisfy the criterion 2 listed in § 2.3.1, namely $\Delta\mathrm{H}(0)\leq\gamma~{}R+10$ kJ/mol, or they are removed because of the domino effect.

Overall, the cleaning process strongly affects the Silicon chemistry, as it will be shown by the modeling simulations of § 4.1. Indeed, out of a total number of 329 removed reactions, 48 involve Si-bearing species (GRETOBAPE-domino-endo-Si). In practice, about 15% of the total removed reactions regard Si-bearing species, and about 9% of the original reactions involving Si-bearing species (526) are removed. Likewise, 6 out of 11 removed species are Si-bearing species (see Tab. 2) and represents the $\sim 10$% of the total original Si-bearing species (61).

Figure 3 shows a graph representation of the reactions involving Si-bearing species and present in the original network, together with the removed and detected species. Figure 4 presents the same reaction network showing the deleted reactions only. First, the 6 removed Si-bearing species (all containing more than 6 atoms) are only formed and deleted by one reaction, so that their removal did not introduce a domino effect on other species. As can be seen from Figs. 3 and 4, these species are in fact at the border of the reaction network, meaning that they have only a minor role in the overall Si reaction network. Second, all the removed reactions involving Si-bearing species have endothermic formation routes.

Digging into the KIDA database to understand their literature origin, we find that almost all these reactions are introduced not because of specific experiments or calculations, but based on educated guesses. Specifically, their rate constants are estimated via the capture theory approach described in Herbst (2006) and Woon & Herbst (2009) for ion-neutral reactions. For example, the formation route of SiC4H, SiC6H and SiC8H is assumed to be Si + CxH2 $\rightarrow$ SiCxH + H, and their constant rates are estimated.

Furthermore, looking in detail at the deleted reactions involving Si-bearing species (Fig. 4) most of them are reactions that increase the complexity of the molecules, namely they increase their number of atoms. This is a general characteristic of the cleaning process: the most affected class of reactions are those that increase the number of atoms of at least one reactant rather than those that reduce it. This result can be explained by the fact that the excess of energy of ”destructive” reactions is generally much easier to guess than that of formation reactions.

Another source of endothermic reactions comes from the cation-anion reactions which involve carbon chain anions, whose data were predominantly taken from Harada & Herbst (2008).

Final remarks

It is of paramount importance to emphasize that the network GRETOBAPE is the result of the removal of endothermic reaction channels. Consequently, some loss/production channels have disappeared from the network, some of which are essential for the production of detected interstellar species. The list of removed reactions needs, therefore, a dedicated systematic (and time-consuming) revision as will be shown in § 3.2.3.

3.2.3 Comparison with the UMIST database

As already emphasized at the beginning of § 3.2, the original network on which we tested the endothermicity of the reactions is based on the KIDA 2014 network and, therefore, the 329 reaction channels removed refer to the KIDA database. In this subsection, we discuss some deleted reactions in comparison with the data reported in UMIST, in which sometimes a more detailed bibliography is reported. It is also worth mentioning that our scope here is not to substitute the precious work done by databases maintainers, but rather to highlight some identified problems.

Reaction Si + C2H2 $\rightarrow$ SiC2H + H

UMIST reports two products, SiC2 + H2 and SiC2H + H, with the second one three times faster than the first. They assumed this branching ratio based on the study at 15 K by Canosa et al. (2001) of the global rate, and based on Kaiser & Gu (2009) for the products. However, the analysis by Kaiser & Gu (2009) shows that the SiC2H + H is indeed endothermic.

Reaction S+ + CH4 $\rightarrow$ H3CS+ + H

This reaction is crucial for the formation of the organo-sulphur products, and we deleted it since endothermic, despite it is considered as barrierless in KIDA and UMIST (with the same parameters). UMIST reports the reference of the experimental data by Smith et al. (1981). These authors assumed that H3CS+ is formed, which is impossible because S+ is a quadruplet while H3CS+ is a triplet and it would be impossible energetically for the spin conservation rules. Indeed, Yu et al. (2020) found that the corrected product seen by Smith et al. (1981) is H2CSH+ (singlet) and not of H3CS+, as corroborated by our calculations that found the channel producing H2CSH+ to be exothermic. Please note that to re-introduce the reaction with the correct reaction channel (H2CSH+) would require to revise all the reactions involving H3CS+ and H2CSH+, which is beyond the scope of the present work.

Reaction CN+ + CO2 $\rightarrow$ NO + C2O+

Experiments sometimes provided wrong products, as in the case of the channel CN+ + CO2 $\rightarrow$ NO + C2O+ which is reported in both KIDA and UMIST and taken from the experimental work by Raksit et al. (1984). However, this channel is endothermic by more than 138 kJ/mol, as found by our computations (which are also validated by thermochemical data111111https://atct.anl.gov/) considering either of the two linear isomers, CCO+ and COC+. In this case, therefore, the problem is not due to the possible isomer structure but to the experiment itself. Indeed, another experiment by McEwan et al. (1983) did not find the NO + C2O+ channel but two other channels, also found by Raksit et al. (1984), i.e. CN + CO2+ and CO + NCO2, both of them exothermic.

3.2.4 Missing destruction reactions

In cold molecular clouds, with few exceptions, neutral species are mainly destroyed by the most abundant cations, namely H${}_{3}^{+}$, HCO+, H3O+, He+ and H+. Therefore, every neutral species should have at least one destruction channel for each of the above cations. Table 5 lists all the neutral species having no destruction channels with the aforementioned cations in the cleaned reaction network. Please note that some reaction channels were removed during the cleaning process. We mark them with a “$\bullet$” symbol. The Tab. 5 list is meant to identify the reactions that need experiments and/or QM computations. In addition, our database of absolute electronic energies, ZPE and $m_{s}$ (Tinacci et al., 2021) can be used to test the possible reaction channels of formation and/or destruction routes.

As an illustrative example, it is worth mentioning the destruction reactions involving SiO. The reaction of SiO with H3O+ $\rightarrow$ H2O + HSiO+ has an enthalpy of +152.1 kJ/mol. A similar endothermicity is found for the reaction SiO + HCO+ $\rightarrow$ CO + HSiO+. However, if the product of both reactions is not HSiO+, as reported in the KIDA database, but its constitutional isomer SiOH+, as reported in UMIST, the reactions become respectively exothermic of -112.2 and -208.7 kJ/mol. However, since in KIDA HSiO+ is produced from H2SiO+, which is in turn produced from H3SiO+, and HSiO+ is involved in a dozen reactions (see Fig. 3), a simple editing may lead to wrong results. A more detailed study on the reliability of the SiOH+ isomer as a product is postponed to a future work by our group focusing on the improvement of the Silicon reaction network (clearly necessary, as highlighted in § 3.2.2).

3.2.5 Missing formation reactions

The aim of the present work is to provide a reaction network as reliable as possible. Identifying important missing reactions is, therefore, an aspect of this goal. In addition to the missing destruction reactions described in the previous subsection, here we list obvious missing formation reactions, namely those forming detected species. To this scope, we used the recent McGuire (2022) census of the detected species in the ISM and extragalactic medium. Table 6 reports the detected species which do not have a formation route in our cleaned network GRETOBAPE. Most of these molecules are also missing from commonly used databases (e.g., KIDA). Hopefully, Tab. 6 can be used to motivate future works and strengthen the interdisciplinary collaboration between chemists and astronomers.

3.3 Reduced network: GRETOBAPE-red

Following the criteria described in Section 2.4, the reduced network GRETOBAPE-red contains 204 chemical species from the 488 existent in the (cleaned) network GRETOBAPE. In addition, the removal of reactions associated with pruned species reduces the size of the network from 7240 reactions to 2810. Compared to GRETOBAPE, GRETOBAPE-red has less than half the amount of reactions and species, implying a significant decrease in computation time when used in any chemical code. This performance boost is seen in the tests we performed on the reduced network using the krome package (see next section, § 4), where the computing time of a typical gas-phase run decreases from 14 to 3 s using the complete and the reduced networks, respectively.

Due to the criteria imposed on the creation of this network, the main discrepancy we expect between the complete and the reduced networks concerns the chemistry of medium-sized (between 3-6 atoms) carbon chains. This is due to removing potentially important destruction routes, i.e. reactions between these medium-sized carbon chains to form larger ones when pruning out C-bearing species with more than 6 atoms. Nonetheless, a complete assessment of the reliability of this reduced network is discussed in the following section, following detailed modeling of a cold and warm case, respectively.

\startlongtable

4 Astrochemical implications

This section illustrates the impact of the new reaction networks on the predicted species abundances, obtained using the original reaction network versus the new cleaned GRETOBAPE and the reduced GRETOBAPE-red networks, respectively. To this end, we run two time-dependent gas-only astrochemical codes: MyNahoon, a modified version of the publicly available code Nahoon (Wakelam et al., 2012)121212Briefly, the modifications concern a more friendly usage of the code, in input and output, and not the core of the code itself., and krome, a publicly available package (Grassi et al., 2014). MyNahoon is used to evaluate the difference in the species abundances predicted with the original and the new GRETOBAPE reaction network in a typical cold molecular cloud (§ 4.1.1) and a typical warm molecular outflow shock (§ 4.1.2), respectively. krome is used to test the reliability of the reduced network GRETOBAPE-red (§ 4.2) in a large parameter space.

4.1 Predicted abundances with the original versus the cleaned reaction network GRETOBAPE

In order to assess the difference between the obtained clean reaction network GRETOBAPE and the original one, we performed two simulations: the first one simulates a typical cold molecular cloud (§ 4.1.1) and the second one a typical warm molecular outflow shock (§ 4.1.2).

4.1.1 Cold molecular cloud

We adopted the parameters of a typical molecular cloud, namely: temperature (gas and dust) equal to 10 K, H nuclei number density equal to $2\times 10^{4}$ cm*-3*, visual extinction Av equal to 20 mag, cosmic-ray (CR) ionization rate for H2 $\zeta$ equal to $3\times 10^{-17}$ s*-1*. The initial elemental abundances are as follows (from Jenkins, 2009, with the heavy metals depleted by a factor 100): He/H = $9.0\times 10^{-2}$, C/H = $2.0\times 10^{-4}$, O/H = $2.6\times 10^{-4}$, N/H = $6.2\times 10^{-5}$, S/H = $8\times 10^{-8}$, Si/H = $8\times 10^{-9}$, Fe/H = $3\times 10^{-9}$, Mg/H = $7\times 10^{-9}$. Finally, we started with the conditions of a pseudo-translucent cloud, namely with hydrogen in the molecular form and the other elements in atomic (ionic or neutral) form, and we let evolve the chemistry up to $1\times 10^{7}$ years.

Figure 5 shows the ratio of the predicted abundances obtained with the cleaned GRETOBAPE and original reaction networks as a function of time, for the species where this ratio is larger than 3 or smaller than 0.3 at any time $\geq 10^{3}$ yr of the simulation and whose abundances are larger than $10^{-13}$ at steady-state. The figure also shows the abundances of these species, to further appreciate the impact. Table 7 reports the minimum and maximum ratios in the time intervals $10^{4}$–$10^{5}$ yr, $10^{5}$–$10^{6}$ yr and at steady-state, respectively.

Overall, 26 species are significantly affected by the reaction network cleaning process, meaning that their abundances have changed by more than a factor three during the evolution using the two reaction networks. The affected species represent only 5% of the species in the clean network GRETOBAPE, which is a reassuring result.

The species impacted by the reaction network cleaning process can be divided into three groups, where we only consider those whose abundance is larger than $10^{-13}$ during the evolution:

The first group (upper panels of Fig. 5) contains species where the post- over pre- cleaning ratio is larger than $10^{3}$ or lower than $10^{-3}$ during the evolution: HCSi, SiO+, CHSi+ and HCCSi. None of these species have been detected in space. Note that five of these six species contain Si and all of them have an abundance lower than $2\times 10^{-11}$. The most abundant one is SiO+, with a predicted steady-state abundance of $8\times 10^{-12}$, $3\times 10^{3}$ larger with respect to the original reaction network. The species whose ratio is the largest at steady-state is CHSi+, when its abundances is lower than $10^{-13}$. HCCSi presents a ratio of 0.3 and an abundance of $5\times 10^{-13}$ at steady-state, while at $\sim 2\times 10^{5}$ yr the minimum ratio value is 0.008. We will discuss in § 4.1.3 the Si chemistry and why it is so impacted by the reaction network cleaning process.

• 2.

The second group (middle panels of Fig. 5) contains species where the ratio is within 10 and $10^{3}$ or $10^{-3}$ and 0.1 during the evolution: SiH, HNSi, HSiNH+, Si+, SiF+, SiH${}_{2}^{+}$ and H2CS. H2CS and Si+ are the only species detected in space belonging to this group. H2CS has a ratio ranging from $0.1$ to $0.03$ during the evolution, namely its abundance is lower in the cleaned network simulation by about a factor 10 at steady-state ($3.4\times 10^{-12}$ instead of $4\times 10^{-13}$). The cleaning process affects identically, with the same ratio profile and increasing the abundances by one order of magnitude at steady-state, the following Si-bearing species: SiH, Si+, SiF+ and SiH${}_{2}^{+}$. HNSi and HSiNH+ share the same behavior, but with a steady-state reduced abundance by about a factor 100. The predicted SO${}_{2}^{+}$ abundance is a factor $\sim 10$ larger with the new network.

• 3.

The third group (bottom panels of Fig. 5) contains species where the ratio is within 3 and 10 or 0.1 and 0.3: SiC, SiH+, HSiO+ and SiN. Three species of this group have been detected in the ISM: SiC, SiN and SiNC. At steady-state, the abundance of SiC is about a factor 10 larger with the new reaction network, while that of SiN is a factor $\sim$3 larger.

Overall, 5 out of 26 species affected by the cleaning process have been detected in the ISM: Si+, SiC, SiN, SiNC and H2CS. Of those, only thioformaldehyde (H2CS) has been detected so far in cold molecular clouds (Irvine et al., 1989; Minowa et al., 1997; Marcelino et al., 2005). We will discuss the case of H2CS in § 4.1.3. Finally, 18 out of 20 species affected by the cleaning process contain Si. We will discuss in some detail this case in § 4.1.3.

4.1.2 Warm molecular outflow shock

Young protostars are known to have spectacular ejections of matter which cause shocks when they hit the surrounding quiet environment. From the chemical point of view, the passage of the shock has two major effects: (i) it heats and compresses the gas and (ii) species frozen on the icy grain mantles are sputtered and injected into the gas phase. As a result, molecular shocks have gas warm ($\sim$80–100 K) and dense ($\geq 10^{5}$ cm*-3*) regions with much larger abundances of some species with respect to cold molecular clouds.

In order to simulate the chemical composition of the gas after the passage of a shock, we adopt a two-step procedure where we first compute the chemical abundances of the cold molecular cloud (previous subsection) and then we suddenly (i) increase the gas temperature and the H nuclei number density and (ii) inject into the gas large abundances of species that are known to be present on the grain icy mantles. Although the physics is slightly different, a similar modeling can be also considered a fair enough simulation of what happens in hot cores and hot corinos, which are regions heated by the central forming star and where the grain icy mantles sublimate (e.g. see the recent review by Ceccarelli et al., 2022). Also, we chose to assume the pre-shock abundances from a steady state cold cloud because our goal is not to reproduce a specific case, but just to approximately understand the impact of the endothermicity correction in a “warm” case. In this respect, the steady state choice has two advantages: it does not depend on the choice of the time (i.e. choice of initial abundances) and it only depends on the change in the reaction network.

In practice, our model has two steps, as follows.

Step 1: We first compute the steady-state chemical composition of a cold molecular cloud at 10 K with a H nuclei number density of $2\times 10^{4}$ cm*-3* and with the parameters described in § 4.1.1.

Step 2: We increase the gas temperature to 90 K, H nuclei number density to $8\times 10^{5}$ cm*-3* and the gaseous abundance of grain-mantle species as follows: CO2/H = $3\times 10^{-5}$, H2O/H = $2\times 10^{-4}$, NH3/H = $2\times 10^{-5}$, CH3OH/H = $6.5\times 10^{-6}$, CH3CHO/H = $3.8\times 10^{-8}$, C2H5/H = $8\times 10^{-8}$, SiO/H = $1\times 10^{-6}$, Si/H = $1\times 10^{-6}$ and OCS/H = $2\times 10^{-6}$. We then leave the chemical composition evolve for $10^{4}$ yr.

The values adopted for the physical parameters and abundances of the injected species are based on the astronomical observations of the prototypical molecular outflow shock L1157-B1, for which our group carried the same model described here (e.g. Codella et al., 2017; Podio et al., 2017; Codella et al., 2020).

The results of the modeling are reported in Tab. 8, which reports the ratio of the predicted abundances obtained with the cleaned GRETOBAPE and original reaction networks in the time intervals $1\times 10^{2}$–$1\times 10^{3}$ yr and $1\times 10^{3}$–$5\times 10^{3}$ yr for species with a ratio larger than 3 or less than 0.3 and abundances larger than $10^{-13}$.

Overall, 22 species are significantly affected by the reaction network cleaning process, namely their abundances change by more than a factor 3 during the evolution from the passage of the shock and $5\times 10^{3}$ yr. The species impacted by the reaction network cleaning process can be divided into three groups with respect to their ratio at $5\times 10^{3}$ yr (e.g. the approximate L1157-B1 age):

The first group contains species whose ratio is larger than $10^{3}$ or lower $10^{-3}$: SiS+, CHSi+, SiO+ and SiN. SiN is the only species in this group that has been detected in outflows and will be discuss in detail in a dedicated paragraph of § 4.1.3. The species with the largest variation due to the cleaning process is SiS+, with an almost constant ratio of $4.2\times 10^{7}$ and reaching an abundance of $\sim 4\times 10^{-11}$ at $5\times 10^{3}$ yr.

• 2.

The second group contains species where the ratio is within $10^{2}$ and $10^{3}$ or $10^{-3}$ and $10^{-2}$: HSiO+, CH3CHCH2, SiH${}_{2}^{+}$, Si+, SiH, HNSi, SiNC+. Two species of this group have been detected in the ISM: CH3CHCH2 and Si+. Only the latter has been detected in outflows and we discuss its case in § 4.1.3. CH3CHCH2 has an almost constant ratio of $\sim 5\times 10^{-2}$ and a maximum abundance of $\sim 5\times 10^{-13}$ at the beginning of the simulation. On the other hand Si+ has a ratio of almost $4.5\times 10^{2}$ that increases the abundance in the cleaned simulation with respect to the original.

• 3.

The third group contains species where the ratio is within 3 and $10^{2}$ or $10^{-2}$ and 0.3 during the evolution: HSS, c-HCCHSi, c-SiC2, SiS, SO${}_{2}^{+}$, HCCSi, Si, HSSH, SiH+, SiC and C3O. Five species belonging to this group have been detected in the ISM: HSS, c-SiC2, SiS, SiC and C3O. Only SiS has been detected in outflows and we discuss its case in § 4.1.3. HSS has an abundance ratio between 0.2 and 0.9 in the $1\times 10^{3}$–$5\times 10^{3}$ yr interval and a maximum abundance of 4.7$\times 10^{-13}$ at the end of the simulation. c-SiC2 at the L1157-B1 age has a ratio of 2$\times 10^{-1}$ and an abundance of 7.2$\times 10^{-13}$. SiS has an almost constant ratio of 0.15 meanwhile C3O ratio oscillates between $\sim 20$ and 2. SiC is the less affected species by the cleaning process with a factor 3 of difference.

As in the cold molecular cloud simulation (§ 4.1.1), the cleaning process mostly affected the Si reaction network. Indeed, 17 of the 22 species present in Tab. 8 are Si-bearing species. Four of them, SiN (Schilke et al., 2003), SiS (Podio et al., 2017) and Si+ (Haas et al., 1986) are detected in outflows. They will be discussed in some detail in § 4.1.3.

4.1.3 Discussion of specific cases

In this section, we discuss in some details the species that are more significantly affected by the reaction network cleaning process.

Thioformaldehyde

In cold molecular gas (§ 4.1.1), the thioformaldehyde (H2CS) abundance predicted after the cleaning process at the steady-state is about 10 times smaller than in pre-cleaning. Here we first summarize the observations and then we comment on what causes the difference in the predicted abundances.

Thioformaldehyde has been detected in cold molecular clouds, notably TMC-1 and L134, with an abundance (wrt to H nuclei) of about $\sim 10^{-9}$ (e.g. Ohishi & Kaifu, 1998). Vastel et al. (2018) and Spezzano et al. (2022) carried out a detailed study of thioformaldehyde toward the prototypical cold prestellar core L1544, reporting a column density of $\sim 6\times 10^{12}$ cm*-2*, corresponding to an average abundance of $\leq 10^{-10}$, depending on where the detected H2CS emission comes from. These authors attribute the relatively low abundance of all the S-bearing species observed in these objects to a general depletion of the gaseous S elemental abundance. Likewise, Esplugues et al. (2022) have carried out a detailed study of H2CS towards several cold starless cores in the Taurus, Perseus, and Orion regions and found H2CS abundances ranging from 0.8 to 14 $\times 10^{-11}$.

As discussed in detail in Esplugues et al. (2022), in cold gas, H2CS is mainly formed by two reactions: S + CH3 and H3CS+ + e-. While H2CS is not directly involved in the cleaning process (Tab. 2), the H3CS+ post-cleaning abundance at steady-state is lower by a factor $\sim 18$ than the pre-cleaned one because one dominant reaction in the original network forming H3CS+, CH$4$ + S+ $\rightarrow$ H + H3CS+, is endothermic by 48 kJ/mol. Please note that the steady-state post-cleaning abundance of H3CS+ is $\leq 10^{-13}$ so that this species does not appear in Fig. 5 and Tab. 7. However, the low abundance ($6\times 10^{-15}$) is counterbalanced by the large recombination rate of H3CS+, whose product is H2CS.

Si chemistry and Si+

As already anticipated, the reaction network cleaning process significantly affects the Si chemistry, with the removal of 48 reactions and 6 species that represent the $\sim 9$% and $\sim 10$% of total Si-bearing reactions and species, respectively. In addition, 18 out of 20 most affected species of the cold molecular cloud simulation (§ 4.1.1) contain Si, of which the vast majority (12) have less than four atoms. Similarly, 17 of the 22 most affected species in the the warm molecular outflow shock simulation (§ 4.1.2) are Si-bearing species.

The six removed Si-bearing species (Tab. 2) are (neutral) chains and are removed because their formation or destruction routes in the original network are endothermic (§ 3.2.2). Their removal very weakly affects two other Si-bearing species, c-HCCHSi and SC3H, respectively (Figs. 3 and 4) and, consequently, the rest of the Si-bearing species.

Therefore, the large impact of the cleaning process on the Si chemical network is due to the removal of the 48 reactions involving Si-bearing species. Figure 4 shows that those reactions connect 36 Si-bearing species out of 61, namely more than half of them. It is, therefore, not surprising that 17 Si-bearing species are among the most affected ones in the warm molecular outflow shock model. For example, SiS+, the most affected species, has three removed reactions connecting it to Si+, SiC3H5 and SiS, respectively. The first two are destruction routes of SiS+; their removal is, therefore, responsible of the augmented predicted SiS+ abundance. A similar argument applies to CHSi+, the second most affected species in the modeling: a reaction of destruction is removed, causing the increase in the predicted abundance.

Probably most important, Si+ has 10 removed reactions linking it to other Si-bearing species. Because of that, the predicted abundance of Si+ is larger more than 400 times the predicted abundance before the cleaning process. This causes a domino effect, propagating as a wave: 11 of the other 16 impacted Si-bearing species are one step away from Si+ (i.e. 1 node distance in Fig. 3), while the other 5 are two steps away. The species one step away are also those whose predicted abundances are the most impacted.

Silicon nitride

Silicon nitride, SiN, is the detected Si-bearing species most affected by the cleaning process (Tab. 8), so we will analyze its case in some detail here. SiN was first detected in the circumstellar envelope of the evolved star IRC+10216 (Turner, 1992) and, later towards the Galactic Center cloud SgrB2(M) (Schilke et al., 2003), probably in the warm shocked gas of the region. IRC+10216 and SgrB2(M) remain the only two objects where SiN has so far been detected, despite sensitive spectral surveys towards the high-mass star forming region Orion (Tercero et al., 2011) and several low-mass star forming regions, including the prototype of the warm molecular shocks L1157-B1 (Ceccarelli et al., 2017).

The predicted SiN abundance is more than 100 times lower when using the post-cleaning reaction network than the pre-cleaning one. Carefully looking at Figs. 3 and 4, one can notice that the SiN major formation route in the pre-cleaning network is the reaction involving HSiNH+, followed by those involving SiC and HNSi+ at a lesser extent. Although this reaction is not affected by the cleaning process, two reactions forming HSiNH+ are removed after it (from Si+ and HNSi, respectively). Indeed, the abundance of HSiNH+ changes by three orders of magnitude (but it is not listed in Tab. 8 because of its low abundance, $\leq 10^{-13}$). About these reactions, as discussed for the HSiO+ case (see § 3.2.4), it is probable that KIDA erroneously considered the product HSiNH+ instead of its more stable isomer SiNH${}_{2}^{+}$ (Parisel et al., 1996; Glosik et al., 1995).

As said, the new reaction network predicts a much lower SiN abundance than before, which may explain the paucity of the SiN detections. For example, Tercero et al. (2011) found a lower limit of [SiO]/[SiN]$\geq 121$. In our post-cleaning modeling, the SiO abundance at $5\times 10^{3}$ yr is $\sim 2\times 10^{-6}$, implying a predicted [SiO]/[SiN]$\sim 5\times 10^{-6}$, therefore largely compatible with the non detection of Tercero et al. (2011).

SiS

Silicon monosulfide, SiS, is another detected Si-bearing molecules affected by the cleaning process. Specifically, SiS has been so far detected towards two warm molecular shocks, Orion KL (Tercero et al., 2011) and L1157-B1 (Podio et al., 2017). The factor between the pre- and post- cleaning predicted abundances is $\sim 0.15$, again probably explaining the paucity of SiS detections. Peering into Fig. 3 and 4, in the pre- and post- cleaning networks, SiS is mainly formed from HSiS+, whose abundance ($\leq 10^{-13}$, so the species not present in Tab. 8) is affected by the cleaning process, with a decrease factor of $\sim 0.06$.

4.2 Predicted abundances with the reaction network GRETOBAPE-red versus GRETOBAPE

To explore the reliability of our reduced network, GRETOBAPE-red, we performed an exhaustive parameter-space exploration and compared the steady-state abundances of all species with those obtained using the complete network, GRETOBAPE.

To this end, we initialized different combinations of parameters that can affect the resulting species abundances in the cold and warm conditions present during the star-formation process. Namely, the varied parameters are: the CR ionization rate for H2, with values of $3\times 10^{-16}$, $3\times 10^{-17}$, and $3\times 10^{-18}$ s*-1*; the gas temperature, with values of 10, 50, 100, and 500 K; and the H nuclei number density, with values of $10^{2}$, $10^{4}$, $10^{6}$, and $10^{8}$ cm*-3*. All possible combinations of these parameters yield 48 different runs, the results of which we will explain in what follows. We used the krome chemistry package (Grassi et al., 2014) to perform the parameter space exploration, considering only gas-phase chemistry. We keep the rest of the input parameters (visual extinction, dust-to-gas mass ratio, grain radius and grain density), equal to those used to perform the run described in § 4.1.1, and left the chemistry evolve up to $1\times 10^{7}$ years.

In all runs, the reduced network GRETOBAPE-red reproduces almost perfectly the steady-state abundances of GRETOBAPE, for all species with abundances $>10^{-10}$, with differences smaller than a factor 3 in general. However, the differences progressively worsen for trace species. A few species with abundances $\leq 10^{-10}$ have differences larger than a factor 3 when comparing the results of GRETOBAPE-red with those of GRETOBAPE. The main species affected are the carbon chains C6, C3, H2CCCO, c-C3H2O, C3O, H2CCN, H2CCO, in addition to CH2NH and O2H. Larger differences might be present for other trace species with abundances $<10^{-15}$ but, as these abundances are already very low, these molecules may not be of interest in the astrochemical simulations where the reduced network GRETOBAPE-red is suitable. As expected, the effect on carbon chains can be attributed to the removal of larger species that can either be part of their destruction or formation routes in the case of dissociative recombination of large cations that could produce the species listed here.

These results overall indicate that the reduced network GRETOBAPE-red is reliable to be used when following the chemistry of species with fractional abundances larger than $\sim 10^{-10}$, while the complete network GRETOBAPE is needed when exploring the chemistry of trace species.

5 Conclusions

Gas-phase reaction networks are a crucial element of any astrochemical modeling. Present day publicly available networks, e.g. KIDA and UMIST, are made up of more than 7000 reactions, of which only a tiny fraction (10–20 %) has been studied in laboratory or theoretical works.

In this work, we present a new gas-phase reaction network, GRETOBAPE, built from the publicly available KIDA2014 (Wakelam et al., 2015) with the addition of several new reactions from more recent studies. The most important novelty is that GRETOBAPE is cleaned by the most obvious source of error, the presence of endothermic reactions not recognized as such. We also present a reduced network, GRETOBAPE-red, to be used in problems that do not require following the chemical evolution of trace species.

To this end, we performed an extended and systematic theoretical characterization of more than 500 species at CCSD(T)/aug-cc-pVTZ//M06-2X/cc-pVTZ computational level, providing for each of them: electronic state, electronic spin multiplicity, geometry, harmonic frequency, absolute electronic energy and dipole moment.

We then computed the enthalpy at 0 K of each neutral-neutral, neutral-ion and cation-anion reaction in the network. Finally, we identified the endothermic reactions not recognized as such, i.e. originally reported as barrierless or having an exponential factor $\gamma$ not taking into account the level of endothermicity.

Following these two criteria, we deleted about 5% of the studied reactions in the original pre-cleaned network leading to the removal of 11 species. The final cleaned network GRETOBAPE consists of 6911 reactions and 488 species, while the reduced network GRETOBAPE-red contains 2810 reaction and 204 species.

We also reported a list of probably missing reactions, namely: (i) neutral species lacking reactions with the most abundant interstellar cations (H${}_{3}^{+}$, HCO+, H3O+, He+ and H+), and (ii) detected species absent in the network.

Using astrochemical model simulations of a typical cold molecular cloud and a warm molecular outflow shock, respectively, we measured the impact of the new cleaned GRETOBAPE on the predicted abundances. Overall, only about 5% of the species abundances are affected by more than a factor 3 with respect to the original pre-cleaned network predictions. Thioformaldehyde (H2CS), Si+, SiN and SiS are the detected species whose predicted abundances are affected by the cleaning process. Of the 5% affected abundances, the immense majority concerns Si-bearing species. We discuss in detail the origin of this large impact on the Si chemistry and conclude that a systematic review of this chemistry is needed.

We also verified that the abundances predicted by GRETOBAPE-red of the species with abundances larger than $10^{-10}$ are not affected by the reduction of the GRETOBAPE network, on a large parameter space. Therefore, GRETOBAPE-red can be used when the goal is to reliably reproduce only the most abundant species and make the simulations less computational demanding.

In conclusion, we carried out the first ever systematic study of the exo/endothermicity of the reactions included in an astrochemical network. The process removing the endothermic reactions not recognized as such and the sink/source species led to the new network GRETOBAPE and its reduced version GRETOBAPE-red. The control of the exo/endothermicity and sink/source species is only the first and obvious step to make the astrochemical gas-phase reaction networks reliable. Inclusion of the reactions of neutrals with abundant interstellar cations is a second obvious but still not complete step. The next step would be to carefully review all the reactions, rate constants and products, which is obviously unfeasible given the too large number of reactions present in the networks. A compromise could be to start with the control/verification of sub-networks that form commonly detected species, like the work done by Vazart et al. (2020) on acetaldehyde. In addition, the lack of gas-phase reactions forming interstellar detected species does not necessarily mean that they do not exist and a lot of work here is still necessary (e.g. Balucani et al., 2015; Skouteris et al., 2018; Sleiman et al., 2018; Cernicharo et al., 2022, just to mention some recent works). Last, but not least, in addition to devote studies on the formation routes, attention should also be dedicated to the destruction routes. For example, the study of the destruction of methyl formate and dimethyl ether by He+ has shown that the rate constants included in the KIDA and UMIST networks are wrong by at least a factor 10 (Ascenzi et al., 2019).

Finally, the database created in this work, and available on-line (§ 6), can be used to rapidly test or search existing or new possible formation and destruction reaction pathways and further refined via experiments or theoretical studies to computed rate constants.

6 Online Database

All scripts and data used in this work are publicly available in the following GitHub link: https://github.com/TinacciL/GreToBaPe_Cleaning. The networks and the database of all molecular species properties are available on the ACO (AstroChemical Origins) Database.

7 Supporting Material

All the reaction networks used and produced in this work, and cited in the main body, are reported in the Supporting Material. We also provide all the .XYZ files of each studied species as well as all the relative information in the Database_Molecules.csv file, following Tinacci et al. (2021). Table 9 summarizes the list of files with the various reaction networks along with a brief description of their content.

This project has received funding within the European Union’s Horizon 2020 research and innovation programme from the European Research Council (ERC) for the project “The Dawn of Organic Chemistry” (DOC), grant agreement No 741002, and from the Marie Skłodowska-Curie for the project “Astro-Chemical Origins” (ACO), grant agreement No 811312. SP, NB, PU acknowledge the Italian Space Agency for co-funding the Life in Space Project (ASI N. 2019-3-U.O). CINES-OCCIGEN HPC is kindly acknowledged for the generous allowance of super-computing time through the A0060810797 project. We warmly thank an anonymous referee whose valuable and expert comments and suggestions helped to clarify the original manuscript. LT is grateful to Jacopo Lupi and Stefano Ferrero for insightful discussions and to the LaTeX  community for the insights on TikZ and PGFPlots packages. Finally, we wish to acknowledge the extremely useful discussions with Prof. Gretobape.