Chemical Evolution of CO2 Ices under Processing by Ionizing Radiation: Characterization of Nonobserved Species and Chemical Equilibrium Phase with the Employment of PROCODA Code

Astrophysical ices are being exposed to ionizing radiation in space environments, which trigger new reactions and desorption processes. In the lab, such processing by radiation has revealed the appearance of several new species and complements the study of the chemical evolution of icy astrophysical scenarios. Here, we develop a computational methodology that helps to clarify the chemical evolution of ices investigated experimentally under photolysis/radiolysis processes until reaching chemical equilibrium (CE). Briefly, the code (named PROCODA) solves a system of coupled differential equations and describes the evolution of the molecular abundances with the irradiation time for ices under processing by radiation. Two experimental ice samples containing pure CO2 and irradiated by two ionizing agents (cosmic rays and ultraviolet photons) were considered prototype systems. Here, we considered 11 different chemical species within the ice (four observed: CO2, CO, O3, and CO3; seven nonobserved or unknown: O, O2, C, C2, C2O, C2O2, and C2O3), 100 reaction routes (e.g., direct dissociation reactions, bimolecular and termolecular reactions) and radiation-induced desorption processes. The best-fit models provide the reaction rates, several desorption parameters, as well as the characterization of the CE phase. At CE, the percentage of nonobserved species in the UV model was almost triple the one calculated in the CR model (which also includes a lot of O and C atoms). The determined values can be employed in future astrochemical models to map chemical evolution embedded species in astrophysical regions under the presence of an ionizing radiation field.

Over the last decade, several authors have undertaken reaction rate studies with applications in astrochemical fields. Most of these works have considered mainly gas-phase reactions (e.g., Wakelam et al. 2010) but some of them also include surface reactions (e.g., Holdship et al. 2018;Shingledecker et al. 2019;Heyl et al. 2020). Such studies have been very useful for astrochemical modeling and the interpretation of astronomical observations. Shingledecker et al. (2019) have presented a model to include radiation chemistry in a rateequation-based astrochemical code, named MONte cArlo COde (MONACO), described originally by Vasyunin et al. (2017) to describe pure O 2 and pure H 2 O ices under processing by kiloelectronvolt protons. In their model, the authors consider that the first-order rate coefficients for the radiolysis process are calculated directly from inelastic/electronic stopping cross sections and the radiochemical yield, and then they estimate the chemical abundances within the ices as a function of time. They also concluded that the radicals produced via radiolysis have a low diffusion rate and react quickly with their neighbors in the ice. Bearing the uncertainties on the reactions rates in the solid phase, the accurate quantification of these parameters is essential to the correct interpretation of ice photochemistry and its posterior use in astrochemical models.
The current manuscript introduces a new computational code, PROgram for solving COupled Differential equations in Astrochemistry (PROCODA), that calculates the reaction rates of suggested chemical systems involved in ice processing by ionizing radiation in the laboratory, directly from the experimental data (abundances or column densities obtained from absorbance IR spectra). In some cases, an interpolation of the data, which should not add or miss any structure on it, might be employed to improve the behavior of the optimization algorithm. In addition to the direct radiation-induced reactions, bimolecular and termolecular collision reactions, the radiationinduced desorption reaction rates, as well as the desorption yield were also calculated. The code also provides abundances for both observed and nonobserved (unknown) species in the IR spectra during ice chemical evolution and characterizes the chemical equilibrium (CE) phase at larger radiation fluence.
Here, we consider that such collision reactions (and their reaction rates) are, in some sense, affected by the type of the incoming radiation, since the radiation allows for the following physicochemical processes within the ices: (i) local heating mainly induced by secondary electrons (stronger in the case of CRs), (ii) molecular excitation, (iii) ionization and dissociation processes, and (iv) ice compaction (mainly in the case of CRs). Additionally, the large energy delivered by ions within the ion track highly enhances the local ice temperature, allowing for changes in the reaction rates since they also depend on the temperature (see also Andrade et al. 2008;Pilling et al. 2010b;Mainitz et al. 2016;Flowers et al. 2019). Two experimental ice samples containing pure CO 2 and that were irradiated by two ionizing agents (CRs analog and UV photons) were considered in this manuscript as prototype systems.
As reported by Öberg et al. (2011), CO 2 is the second most abundant frozen molecule with respect to H 2 O ice as observed toward low-mass protostars (29%), and quiescent molecular clouds (38%), and it has an abundance similar to CO ice toward high-mass protostars (13%). The formation of CO 2 ice is still under debate, and its detection toward molecular clouds and star-forming regions indicates that it can be formed without the presence of ionizing sources. Experimental studies by Grim & d'Hendecourt (1986) concluded that one of the direct formation routes of solid CO 2 via CO + O has a significant activation barrier (290-1000 K; d 'Hendecourt et al. 1985;Roser et al. 2001), and thus would not be efficient on cold grains. Successful experiments performed by Ioppolo et al. (2011) have shown that CO 2 is efficiently formed via CO + OH → CO 2 + H. Moreover, this reaction is a factor of 10 more efficient than the CO + O reaction (Ioppolo et al. 2013). Adding the CO + OH reaction to chemical models of molecular clouds enhances the abundance of CO 2 (Garrod & Pauly 2011). Another reaction pathway starting with formaldehyde (H 2 CO) has been proposed by Minissale et al. (2015), i.e., H 2 CO + O → CO 2 + H 2 . Alternatively, irradiation experiments involving solid H 2 O and CO indicate that CO 2 is readily formed (d 'Hendecourt et al. 1985;Gerakines et al. 1996;Jimenez-Escobar et al. 2016), which can explain its formation in stellar nurseries. Infrared Space Observatory observations of CO 2 at 4.25 and 15.5 μm analyzed by Gerakines et al. (1999) indicate that CO 2 ice may reside in two chemical environments, namely, polar (H 2 O-rich) and nonpolar (CO 2 -rich). However, good profile fits are also found when nonpolar components are replaced by thermally annealed ices, thus indicating that the chemical evolution of CO 2 is mostly due to thermal annealing compared to energetic processing of the ice mantle. Pontoppidan et al. (2008) analyzed the CO 2 band at 15.5 μm observed with Spitzer/ InfraRed Spectrograph as part of the program "Core to Disks" (c2d) Legacy program (see also Evans et al. 2009). A similar conclusion was found compared with Gerakines et al. (1999), that two-thirds of the CO 2 ice is found in H 2 O-rich environments, whereas one-third resides in CO-rich ice mantles. Signatures of ice thermal processing were also found, although via the distillation process, namely, that CO desorbs from the mixture CO:CO 2 after warm-up, leaving behind pure CO 2 , as observed by the split absorption band at 15.5 μm.
Finally, experimental studies addressing the CO 2 ice UV photolysis indicate that CO is the main product, and CO 3 and O 3 are also formed (Gerakines et al. 1996). When solid CO 2 is mixed with CH 4 , the abundances of CO, CO 3 , and O 3 decrease by an order of magnitude, which shows the importance of the chemical environment on the formation of photoproducts (Öberg et al. 2009). Although these experiments improve the qualitative knowledge of the CO 2 photodissociation process, astrochemical models are required to completely quantify the CO 2 photochemistry. This knowledge is achieved by calculating the photodissociation branching ratios and rates. In this context, the rates of possible back-reactions and other CO 2 reformation routes have to be quantified.
Section 2 describes the computational methodology in this work and some details of the employed experimental data (taken from literature). Results are shown in Section 3 with an emphasis on the estimate of reaction rates and molecular abundances at equilibrium chemistry in the two sets of data considered (one from experiments employing CR bombardment and the other from experiments employing UV irradiation). Astrophysical implications were also provided. The conclusions and final remarks are listed in Section 4.

Methodology
In this manuscript, we developed a computer code called PROCODA to solve a system of coupled differential equations (SCDE) to describe the chemical evolution of typical astrophysical ices during the processing by ionizing radiation as a function of time. The solution provides the reaction rates (ks) of each reaction considered in the chemical network. Here, we consider a set of coupled differential equations containing 100 chemical reactions (plus a reaction set to take into account the molecular desorption from ice to gas phase induced by radiation) to map the chemical evolution of selected molecules in ices containing initially only the CO 2 molecules under the presence of ionizing radiation (CRs and UV photons).

The Coupled Chemical Reactions
The proposed reaction network involves 11 molecular species (the parent species CO 2 and the new produced species CO, O 3 , CO 3 , O, O 2 , C, C 2 , C 2 O, C 2 O 2 , and C 2 O 3 ). Two experimental data sets were considered: (i) CO 2 ice at 13 K bombarded by 52 Ni ions (taken from Pilling et al. 2010a), and (ii) CO 2 ice at 8 K irradiated by UV photons of ∼ 10 eV (taken from Martín-Doménech et al. 2015). It is worth noting that, in both cases, the experimental data shows the abundances of only four observable species (CO 2 , CO, O 3 , and CO 3 ). With the current methodology, we are able to map and characterize, for the first time, not only the abundances of the observable species but also the abundances of other proposed seven species nonobservable (unknown species) in the experimental IR spectra, the newly produced species O, O 2 , C, C 2 , C 2 O, C 2 O 2 , and C 2 O 3 .
It is worth mentioning that since the sample temperature is very low, in the absence of incoming radiation, the chemical abundances would virtually not change with time. Therefore the chemical reactions modeled and described in this manuscript are in some sense induced by radiation. The incoming radiation basically allows: (i) local heating (mainly in the case of CRs); (ii) excitations; (iii) ionization and dissociation processes; and (iv) ice compaction (mainly in the case of CRs). All of these processes affect the reaction rates of embedded species. Since the studied ices are kept at constant low temperature and are exposed to constant radiation flux, the chemical abundances within the ice changed with the passing of time until CE was reached at longer radiation fluence. The CE phase occurs when virtually the molecular abundances do not change anymore with time even during additional irradiation (considering the same incoming radiation flux and ice temperature). Additional details on CE are provided below. Figure 1 shows the selected reactions considered in the present chemical network, involving potential products of the CO 2 irradiation by CRs and UV photons (the left panel shows the main reactions involving CO 2 , and the right panel shows the main reactions involving CO). The full reaction set is listed in Table 1. The segmentations in each panel indicate specific domains of collisions between the reactants (direct reactions, collisions with twins, with daughters, and with granddaughters). Most of the reactions consider bimolecular collisions, but some termolecular collisions were also considered.
The current model describes the evolution of the chemical abundance of 11 molecular species and the number of species desorbed from ice to gas phase due to the incoming X-rays in frozen CO 2 -rich ice. In total, 100 reactions were considered in the ice phase (half in the forward direction and the other half in the reverse direction to sustain the CE hypothesis at large fluences). Additionally to these chemical reaction sets, reactions with the molecular desorption processes were also considered.
The chemical reactions selected for this work consider reactions only between neutral species, radicals, or atoms, some of them already suggested in the literature (mostly at gas phase; see Table 2). Therefore, the current reactions consist of a subset of the uncountable number of the reactions that might occur within the ice bulk due to the possibility of different reactants as well neighbor species (which also increases with irradiation time, until CE is reached). Reaction with ions and electrons was not considered in this work. A good review on reaction networks for interstellar chemical modeling can be seen at Wakelam et al. (2010).
The typical equation in the chemical coupled system we solve by employing PROCODA has the following parameters: where dN i /dt is the column density changes of a given molecule i along the time t, and k values are the reaction rates for different processes. The first brackets indicate consumption (or destruction) processes, and the second brackets indicate the production terms. N i indicates the column density of a given i species and the values N a , N b , and N c indicate the column density of species a, b, and c, respectively, which might participate in the reaction to produce or consume the respective species i. In Equation (1), the values k d1 , k d2 , and k d3 indicate reaction rates for direct dissociation (or during the impact of CR/UV), bimolecular collision, and termolecular collision, respectively. The values k p1 , k p2 , and k p3 indicate reaction rates for the direct production of i from a given species a, or from bimolecular collision and termolecular collisions also involving species b and c, depending on the case. The term DES i (t) in Equation (1) is the differential column density desorption, i.e., the number of molecules (or atomic species) that desorbs from ice to gas phase per centimeter squared and per second due to incoming radiation that is  Notes. a Input parameters: p1 (SF's weight for the CO2), p2 (SF's weight for the CO), p3 (SF´s weight for the O3), p4 (SF's weight for the CO3), p5 (SF's weight for the mass conservation), p6 (SF's weight for the desorption similarity), p7 (SF's weight for the slope similarity criterion; hypothesis for the chemical equilibrium (CE) stage), pDES (estimated input for the desorption yield). b Additional input parameters for the CR Model are: N_CO2_initial (molecules cm −2 ) = 2.14 × 10 14 , phi (ions cm −2 s) = 2 × 10 9 , d_ice (cm) = 1 × 10 −5 , s_ice (cm 2 ) = 0.8. c Additional input parameters for the UV Model are: N_CO2_initial (molecules cm −2 ) = 2.00 × 10 17 , phi (photons cm −2 s) = 2 × 10 14 , d_ice (cm) = 8.1 × 10 −6 , s_ice (cm 2 ) = 0.8. d Outputs parameters: SF = Score function; sumCHI2 = summed chi-squared function for the observed species; TotDES = total summed desorption yield; MSS (%) = Mass similarity criterion (similarity between initial mass and final mass) in percentage; SSC (%) = Slope similarity criterion in percentage between summed observed mass and summed unknown mass (considering the last 200 s); hypothesis for the CE stage; Average reaction rates: For direct (k d ), bimolecular (k A+B ) and termolecular (k A+B+C ) reactions. Other output parameters are: individual ks, individual molecular abundance at CE, individual desorption column density and individual intrinsic molecular desorption in units of s −1 .
obtained by the equation: where k des,i is the intrinsic desorption rate (time independent), in units of hertz, and Ω i (t) is the dimensionless surface coverage of the species i as a function of time, which ranges from 0 to 1 (e.g., Ω CO2 (0) = 1 and Ω daughters (0) = 0; see also Pilling et al. 2010b). The desorption itself was considered a first-order type reaction. By summing all individual column density desorption, we calculate the total column density desorption. From such a quantity and by also considering the total incoming radiation fluency and the sample area, we calculate the total desorption yield, in units of molecules per ion or molecules per photon (depending on the case).

The PROCODA Code
The code was written in Python language by using packages Numpy, Scipy, and Matplotlib, 3 and was structured in three blocks: the first one was designed to handle input data, the second block was designed to find the best solution for the coupling differential equations set, and the last block was designed to make final calculations, including the molecular abundances at the CE (at large radiation fluences), the individual desorption column densities, the average values for the rates, and generating the outputs and plots.

First Block: Reading Input Data and Parameters
The first block was designed to read the experimental data set, to create an interpolated set for the experimental data, and to read the input parameters (e.g., the type of species considered, the initial column density of the ice, sample size and thickness, and initial guess for the reaction rates and molecular desorption). The input experimental data set is a typical ASCII file with the experimental column density and exposure time in each line (no experimental error bars were considered, but they are usually very small with respect to the trends). To allow for better accuracy in the minimization function (to improve the behavior of the optimization algorithm), we solve the chemical coupled system considering a larger and interpolated data set (generated by employing a linear or/cubic-spline interpolation from the lowresolution experimental data sets). The sizes of the experimental data (for each observed species) considered in the CR model and UV model were 12 and 11, respectively. In the current methodology, for each observed species in the CR model, 710 nonequally spaced interpolated data points were considered from 0 to 5000 s (∼1.4 hr). For each observed species in the UV model, 1010 nonequally spaced interpolated data points were considered from 0 to 14382 s (∼4 hr). No structure to the modeled data was added or missed due to the interpolation processes. See Appendix A for more details of input data interpolation.

Second Block: Solving the Coupling Chemical Reaction Set (the Exploratory and Refinement Stages)
To find the best solution for the chemical equation set, we develop a methodology in two stages (the exploratory and refinement stages) to allow the process to navigate through different local minima and to find the global minimum during the minimization process of a predefined score function (SF), which compares the model with the experimental data (details of the SF are given below). In the minimization process, we employ the L-BFGS-B 4 algorithm using the modules SCIPY.INTE-GRATE.ODEINT and SCIPY.OPTIMIZE.MINIMIZE. 5 The method L-BFGS-B is an extended version of the L-BFGS algorithm for bound constrained minimization (typically for variables with constraints between 1 and a number much less than 1). It is an optimization algorithm in the family of quasi-Newton methods that approximates the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm using a limited amount of computer memory (Byrd et al. 1995;Zhu et al. 1997;Nocedal & Wright 2006;Morales & Nocedal 2011). The L-BFGS-B algorithm is a very efficient algorithm for solving large-scale problems. The method works by identifying fixed and free variables at every step (using a simple gradient method), and then using the L-BFGS method on the free variables only to get higher accuracy, and then repeating the process. However, the BFGS method is one of the most famous quasi-Newton algorithms for optimization (pointed out by Dai 2002). Depending on the complexity of the reaction set and the initial input conditions, the convergence to a globally optimal solution might be difficult. Other examples of coupled chemical systems solved by employing differential equations with applications in astrochemistry are shown by Kim et al. (2010), Bennett et al. (2006), Bennett et al. (2007), Jamieson et al. (2006), and Frenklach et al. (1992).
In the exploratory stage, we start with the initial guess for the reaction rates (typically of the order of 1 × 10 −6 ) and run, recursively, the minimization process of the SF considering as the input of a given run the best solution of the previous run. In this stage, we employ a multiplicative factor (which considers a random value within different ranges depending on the step) for the k values at the beginning of every loop. This multiplicative factor, at a specific moment during the run, varies from 0.001 to 1000, which allows the system to skip a given local minimum in the minimization processes and permits the search of a new minimum in the SF.
In the refinement stage, we start with the best solution obtained in the exploratory stage and proceed with the solution refinement protocol by considering several cycles of solving the differential equation solutions and minimization of SF, in which every step, a lower random multiplicative factor was involved to fine-tune the best solution. In this part, the initial solution for the system considers the previous best solution multiplied by a random factor (from 50% to 5%) that, after several iterations, yields a smaller value of the SF, as well as, in the summed chi-square value (also calculated separately from the SF). In the exploratory stage, we navigated through different local minima to find the global minimum in the considered SF, and in the refinement stage, we assumed that we were already in a global minimum and just finetuned the values of the reaction rates providing the best solution for the coupled chemical system employing the current.
A diagram with the schematic flow chart of PROCODA is given in Figure 2 (top panel). The exploratory and refinement stages are shown in the black and red colors, respectively. In this figure, SCDE is the system of coupled differential equations (the "heart of the code"), SF is the score function calculation and minimization processes (the "brain of the code"), and L_stop is the total number of loops inside the code. Some diagnostic plots are presented at the bottom of this figure. The main panel at the right presents the evolution of the SF and the summed chi-square for the observed species in a typical code run, sumCHI2 (as we will discuss further, this parameter can be obtained from Equation (3), considering p1 = p2 = p3 = p4 = 1 and p5 = p6 = p7 = 0). From this panel, we observe that, as the program runs, the system passes through several local minima, reaching a global minimum at a large loop number. In particular, this figure shows a calculation in which the best fit was obtained at the loop equal to 42 (some info about outputs is shown in the figure header). Usually, the SF is larger than the summed chi-square for the observed values. The four panels at the right present the evolution of the summed desorption yield and average values for different types of reaction rates after each calculation loop in a typical run. Such quantities have small variations as the SF reaches a global minimum during the minimization processes. The uncertainty on the calculations was estimated to be below 20% (done by comparing different models with close input parameters).
The SF employed in the minimization processes to find the best solution for the reactions rates is given by where the MSC f is the total mass similarity criterion (SC; =mass conservation criterion), calculated considering the similarity between the initial column mass of the system and the total column mass of the model at final modeling time (largest time = 5000 s), MSCo f and MSCo m are the column mass SC between the experimentally observed column mass and the observed column mass in the model at the final modeling time and at the middle of the modeling time (2500 s), respectively. The column mass of a given molecular species (in units of daltons per centimeter squared) is easily obtained by multiplying the column density, in units of molecules per centimeter squared, by its respective molecular mass, in units of daltons per molecule. The parameter DSC in Equation (3) is the desorption similarity criterion, calculated considering the similarity of the input expected desorption yield (based on or taken from the literature) and the total molecular desorption yield computed by the model (summed individual molecular desorption divided by total incoming ions). Additionally, the parameter SSC is the slope similarity criterion, which is calculated between the similarity of two vertical segments on column mass plots considering the last 200 s of modeling (see also Figures 3(b), and 4(b)). Each vertical segment is obtained by the difference between the summed observed column mass and summed unknown column mass. The large similarity in this criterion implies that the column masses of both observed and nonobserved reach a sloped plateau at a large irradiation time. In other words, that system has reached the CE stage. The fact that the reached plateau is not truly horizontal is due to the desorption of molecules to gas phase (larger in the case of CRs), which decreases the column mass of frozen species. All of these criteria were built to range from 0 to 1 (1 = maximum similarity). The first four terms in Equation (3) include the chisquare function between the experimental concentration data of the observed molecules (CO 2 , CO, O 3 , and CO 3 ) and their modeled concentration (here, oi data represents the experimental column density data of observed species i and oi model is its respective modeled column density; the letter o in this term is an abbreviation for the observed species).
The parameters p1 to p7 employed in Equation (3) (given in the input) are the weight (dimensionless) of each term described in the equation to allow for the search of the best solution during the computational minimization processes (considering the L_BFGS-B algorithm). These parameters are fine-tuned manually to allow for the lowest value of the SF during its minimization process, as well as a very small value for the summed chi-squared function for observed molecular species at the end of the calculation. The p1 to p4 are the employed weights for the chi-squared function for specific Figure 4. Evolution of column density (panel (a)) and selected column mass (panel (b)) as a function of time obtained for the model of pure CO 2 ices at 8 K processed by UV photons (10 eV; considering data from Martín-Doménech et al. 2015). Lines are the best-fit model. The two vertical segments around time 14,000 s were considered to calculate slope SC and to guarantee the CE stage. See details in the text. observed molecular species (important for the less abundant species), large p5 are employed to force mass conservation of the system, large p6 also guarantees that the model desorption yield is close to the estimated value based on or taken from the experiments, and large p7 guarantees the CE stage at final fluences. In Appendix B, we present a comparison between models obtained with different input parameters to discuss the employed technique and the importance of using the considered SF instead of only an ordinary summed chi-square function in the minimization process to find the best solution of the reactions set. It is worth noting that the computational code was not designed to fit directly the low-resolution experimental data set (taken directly from the published papers).
Summarizing, the SF considered both the experimental data and some input parameters (the weights p1 to p7 and, similarly, the criteria for the mass MSCs, for DSC, and for slope SSC) to allow for the minimization process to converge to a global minimum during the calculations. This provides the best values for the reaction rates that describe the evolution of the experimental system with time. It is worth mentioning that besides the calculation of SF in every step, we also calculate the summed chi-square function for the observed molecular species (sumCHI2), and both values are listed in the outputs of the code. Finally, the employ of the minimization processes directly and exclusively on the sumCHI2 does not provide the best solution for the chemical systems with the current computational methodology. More details on this issue are given in Appendix B.
The similarity criterion (SC) between a given two numbers (n1 and n2), employed in the calculation of several SCs written in Equation (3), is given by the following equation: and ranges from zero to one (1 = maximum similarity). In the case of two numbers (n1 and n2) with different signals, we add an equal number to both to allow the negative ones to become equal to 1 and then calculate the SC value. In the case that one number is zero, we add one to both numbers and then calculate the SC value. In the case that both numbers are zero, the SC value is set to 1. The criterion for model convergence and to minimize degeneracy in the solutions involves the mass conservation, the similarity between calculated desorption values and the input desorption values (taken by literature), and the presence of CE at large irradiation time (or larger radiation fluences). The CE stage has been characterized experimentally by different authors in the literature (e.g., Vasconcelos et al. 2017c;Pilling et al. 2019, and references therein) and briefly assures that the molecular abundances of the system do not change significantly with time even when more radiation is added in the systems considering a constant ice temperature.
Typically, a code run has 60 main loops (customized if needed), 15 in the exploratory stage plus 45 loops in the refinement stage. Each main loop has its own internal loops in the routine of minimization and in the routine for solving the system of chemical equations (see Equation (1)). To obtain better accuracy in the solution, the models are run two times (the second run considers as input the best solution of the first run). The calculations were performed in a 64 bit Windows 7 based computer, with Intel i7-2600 CPU at 2.40 GHz with eight cores and 16 Gb of RAM, running python 3 (Anaconda distribution) via Jupiter Notebook at Chrome webpage browser. Each calculation run takes around 900 minutes.

Advantages and Limitations of the Code
The main advantage of the current methodology is to map, in addition to the observed species, the nonobserved (unknown) species in the IR spectra obtained in the experiments employing ionizing radiation at astrophysical ice analogs. In addition to the calculation of direct radiation-induced reactions, the code also calculates bimolecular and termolecular collision reactions, radiation-induced desorption rates, as well as the total molecular desorption yield. The code also provides abundances for both observed and no-observed (unknown) species in the IR spectra during ice chemical evolution and characterizes the CE phase at larger radiation fluence.
It is worth noting that the current code, in addition to its use in the solid-state astrochemistry field as shown in this study, can also be employed (with some minor changes) to quantify the reaction rates and the abundances in the CE stage in other types of chemical reactors (e.g., considering gas phase or liquids or temperature processing).
One of the limitations of the code is related to the small number of chemical species and reaction routes considered within the chemical network. The employ of the Quadrupole Mass Spectrometer in the experimental setup to map the desorption of nonobserved species in the IR might also help place a constraint on the molecular species that could exist within the ice. This should be investigated in future studies. The second limitation is related to the fact that the current code does not take into consideration the activation energies of the selected reactions, which might induce some artificial favoritism in the reaction pathways. The next version of the code will try to handle this activation energy issue. However, it is worth noting that the energy input in the current models (mainly in the CR model) is high enough to allow reactions with very large activation energies. Finally, the third main limitation is related to the degeneracy of the solutions and/or convergence to a global minimum. As pointed out by Guitton (2004) and also by Dai (2002), the difficulty in such high-dimensional chemical network calculations is to evaluate whether convergence to a global optimal solution is achieved, or whether multiple reasonable solutions are present.
To minimize these degeneracy and convergence issues, the current methodology employs a recursive solution of the chemical system (which involves the exploratory and refinement stage with random multiplicative factors applied to the input values during a typical code run) and employing the SF with CE criterion (as well as other SCs) as described previously. Additionally, we run the code twice (monitoring the values of the SF and sumCHI2) to allow for the search of the global minimum (or at least the best possible solution for the system). It is worth noting that the global minimum might not be necessarily a single point in parameter space, a large portion of parameter space can have an approximately equal value for the SF; this should be investigated in future research.
To evaluate the uncertainty of the current methodology, we calculate the relative standard error (RSE) of reaction rates and other output parameters considering 30 different calculations of the CR Model with the same initial input (see Appendix C). We can consider, as a good estimate (within a safety margin), that the current methodology to calculate the solid-state reaction rates for the irradiated ices has an error below 20%.

Evolution of Molecular Abundances and Reaction Rates
The best-fit models employing the PROCODA code to describe the chemical evolution of pure CO 2 ices under the presence of ionizing radiation are shown in Figures 3 and 4, for the employed CR (Pilling et al. 2010a) and UV (Martín-Doménech et al. 2015) experimental data, respectively. In both figures, panel (a) indicates the evolution of column density as a function of time, and panel (b) shows the evolution of selected column mass (e.g., summed column mass of observed species, summed column mass of nonobserved species (unknown species), summed column mass of desorbed species, and total summed column mass). The experimental data are represented by filled black symbols (only the four observable species CO 2 , CO, O 3 , and CO 3 ) and the lines are the models. The two vertical segments presented in panel (b) (at larger times) quantify selected differences between the summed column mass of observed and nonobserved species and are employed in the calculation of slope SC to guarantee CE of the system at larger radiation fluences. In all figures, important output parameters for the best-fit models are displayed in the header.
The models have a very good match with the experimental data and also reproduce the CE behavior (a slighted sloped plateau for the column densities observed at large fluences), for both observed and nonobserved species. In both studied systems, the main changes in the molecular abundances occur at the beginning of the irradiation phase (roughly in the first 1000 s of irradiation). With the employed methodology, we quantify not only the abundances of observable species (in the IR spectra irradiated ices) but also the abundances of other nonobservable species (also named as unknown species). Such characterization helps us to better understand the experimental data.
A direct comparison between both experimental data sets allows us to observe three main important differences: (i) the molecular desorption from ice to gas phase by CRs (also known as sputtering) is much larger than the one observed in ice irradiated by UV, as expected; (ii) the final abundance of the daughter species CO is higher in experiments employing UV (which might be related to the previous item or/and due to small desorption temperature of CO); and (iii) the amounts of other observed daughter species O 3 and CO 3 are also larger in the case of the UV experiment. Table 1 lists the input and output parameters of the best-fit models described in this manuscript. Once more, we reinforce that the best solution is the one that has, simultaneously, mass conservation, the presence of CE stage at final fluences, a good similarity between the calculated desorption yield with literature values, and a low value for the summed chi-square function for the observed molecular species (sumCHI2). It is worth noting that calculations considering only the sumCHI2 in the SF (considering simultaneously p1 = p2 = p3 = p4 = 1 and p5 = p6 = p7 = 0 in the SF), despite good chi-square values, do not produce CE scenarios or/and mass conservation. In the Appendix B, we present eight additional models (bad solution models) for a comparison purpose with the best-fit CR model. These models were calculated considering different weights in the scores function or in the CE criterion or different desorption yield depending on the case. Table 2 presents the complete chemical reactions network employed in the current model to map the chemical evolution of pure CO 2 ices during processing by ionizing radiation, with the best-fit values for the reaction rates for each reaction. In this table, we present the modeled parameters employing data from CR analog experiments (13 K ice bombarded by 52 MeV Ni ions; Pilling et al. 2010a) and data for the UV irradiation experiment (8 K ice irradiated by ∼10 eV photons; Martín-Doménech et al. 2015). The reaction order considered for each chemical reaction is also given (considered in this work equal to the number of reactant molecules). Values for the intrinsic desorption rates are listed at the bottom of the table, and they will be discussed later on. The last line presents the total molecular desorption rate and desorption yield obtained in the models. For some reactions, the model provides only upper limits of reaction rates (values with the symbol "" in the table), which indicate that values below this threshold virtually do not affect the solution of the system. Table 2 also lists references for some reactions, mostly from gas phase, taken from different literature databases. As discussed before, in this work, only reactions involving neutral species are considered.
From the best-fit models, we can draw some average values for the reaction rates of direct dissociation (the first-order reactions, triggered by the incoming radiation, not considering desorption reactions), and bimolecular and termolecular reactions for the two studied ices. The average radiolysisinduced (dissociation) rate and its respective average radiolysis-induced dissociation cross section, considering the experimental radiation flux in each case, were In the case of the CR model, the average radiolysis-induced dissociation cross section is about one order of magnitude lower than the estimated value from the experimental data, with the employ of a semiempirical methodology such as an associative exponential equation to best fit the data (Pilling et al. 2010a). In the case of the UV experiment, the value determined by Martín-Doménech et al. (2015) for the destruction of CO 2 was 9.5 × 10 −18 cm 2 , about less than four times higher than the average value obtained with the current methodology. For comparison purposes, we observe that the values determined for photodissociation cross section at the gas phase are much lower (e.g., 0.6 × 10 −19 cm 2 ) than the ice phase, as reported by van Dishoeck et al. (2006). As also discussed by Pilling et al. (2011), the presence of neighbors (chemical environment influence) increases the photodissociation cross section due to the opening of other reaction pathways during the excitation/ionization processes.
It is worth noting that the considered radiation-induced reactions in the models are the ones triggered by the direct interaction between the radiation (CR or UV) with the target species (see Table 1). In this group of reactions triggered by direct radiation, we should also consider the related processes such as the ones promoted by secondary electrons within the ice, and electromagnetic-induced excitations (e.g., secondary UV or X-ray photons in the case of CR irradiation; Prasad & Tarafdar 1983) to the frozen molecular species. Therefore, such first-order reactions represent an average set of reactions triggered by radiation, which include: (i) the direct processing of a given molecular species by the incoming radiation, (ii) electron-induced reactions, and (iii) reaction triggered by excited states, all also considering the different molecular neighbors of the given species (the effect of the average dielectric filed). For a discussion on the effects of molecular neighbors in some properties of embedded species, such as its reactivity and IR spectra, see also Pilling et al. (2021) and Alves et al. (2021).
The average bimolecular and termolecular reaction rates in the processing of pure CO 2 ices by ionizing radiation obtained in this work were The best-fit models also show that, typically, the average bimolecular and termolecular rates in the case of the CR experiment are roughly double of those obtained in the case of the UV experiment. Such differences might be related to the presence of a large number of secondary electrons within the ices, as well as the large energy delivered by ions within the ion track, which enhance the local temperature of ice, inducing changes in the reaction rates (see also Pilling et al. 2010b;Andrade et al. 2008;Mainitz et al. 2016). These average values might be employed in future astrochemical models to map chemical evolution embedded species in astrophysical environments under the presence of an ionizing radiation field. Future investigation will calculate the dependence of temperature and initial composition of ices in these average rates. The current work, in combination with the temperature dependence investigations, will help to generate a big picture of the reaction network within ices during processing by radiation in the astrophysical analog scenario.
Despite the reduced number of reactions, with the current methodology, we can model the chemical evolution of the CO 2 ice system under processing by radiation with a considerable accuracy comparing the experimental data and provide, for the first time, an estimate for bimolecular and termolecular reactions in astrophysical ice analogs under processing by radiation. Future work may have an extension to the reaction set to include several other species as well as electrons and ionic reactants.

Branching Ratio of Selected Reactions within CO 2 -rich Ice Processed by Cosmic-Ray Analog
From Table 2 we calculate the branching ratio (in percentage) of selected direct radiation-induced reactions and selected bimolecular reactions involving CO 2 or CO. The values are presented in the reaction scheme below for the experimental data employing CR and UV. The branching ratio of a given reaction was obtained by dividing the value of the rate constant of the given reaction pathway by the sum of the rate constant of the other reaction, which has the same reactants. Analyzing these branching ratios, we observe that some reaction pathways do not have a clear preference to occur in the presence of a specific radiation source, and others do. For example, the models show that the reaction CO 2 + CR yields preferentially CO + O, but when UV is employed, there is no clear preferentiality among the studied reaction pathways. In both cases, the reaction CO 3 + CR/UV yields preferentially CO + O+O instead of other studied channels, and the reaction O 3 + CR/UV yields mainly O 2 + O rather than three O atoms.

Direct radiation-induced reactions
The reactions considering the dissociation of daughter species C 2 O 2 and C 2 O 3 have different main reaction routes depending on the employed radiation. The reaction C 2 O 2 + CR yields preferentially 2 CO in contrast with the reaction C 2 O 2 + UV, which yields preferentially C + CO + O. The reaction C 2 O 2 + CR yields preferentially CO + CO + O in contrast with the reaction C 2 O 2 + UV, which yields preferentially C 2 O+O 2 .
Since the environment affects the reactivity of reactants, as discussed for example by Pilling et al. (2021), Alves et al. (2021), and references therein, it is expected that the reaction branching ratio within ice bulk or even at ice surfaces would be very different from that at gas phase (not only for direct reactions but also for bimolecular and termolecular reactions). This is observed, for example, comparing the photodissociation of CO 2 by UV in the gas phase, which yields about 16% of the reaction channel containing O + CO (Bittner et al. 2020) with the same reaction route in the ices (this work), for which the calculated branch ratio was roughly double.
The precise explanation about why specific incoming ionizing radiation induces such preference (or not) in some reaction routes is beyond the scope of this manuscript; however, we can suggest that such phenomenon might be related to the presence of energy input necessary to overpass a given activation barrier for some reactions. This energy input might be linked to the local temperature enhancement induced by the incoming radiation, which is different when we compare the different types of ionizing radiation (very important when CRs are considered) and the number of secondary electrons (and X-rays) produced/released within the ice when specific radiation is employed.

Bimolecular reactions involving collisions with CO 2
As observed for some radiation-induced reactions, depending on the type of incoming radiation source (which also reflects the amount of energy excess available within the ice), some bimolecular reaction channels have a preference to occur. Selected bimolecular reactions involving collisions with CO 2 with calculated branching ratio are presented in the reaction scheme below, for the experimental data employing CRs and UV radiation on pure CO 2 ices.
In the calculation of the branching ratio for the bimolecular reactions presented below, we do not take into consideration the estimated uncertainty of the employed methodology (see Appendix C) but since the average uncertainty for the bimolecular reactions rates was 12.1%, we consider, as a good estimate, that the average uncertainty for the bimolecular branching ratio might have roughly a similar value. However, it is worth noting that despite most of the uncertainties being below 10%, some rates presented much larger values such as the k7 (∼45%), k11 (∼90%), k31 (∼40%), k35 (∼95%), k37 (∼60%), and k58 (∼40%); therefore, the branching ratio value considering these reactions should be considered with caution. From the selected bimolecular reactions, we observe that, in both experiments, the reaction between CO 2 + O yields preferentially CO + O 2 , and the reaction between CO 2 + O 2 yields mainly the CO + O 3 species, among the considered different reaction routes. However, in the case of the reactions between CO 2 + O, CO 2 + O 3 , and CO 2 + CO 2 , depending on the type of incoming radiation, we have a preferential reaction route. For example, when CR is employed, the preferential routes for these reactions were CO 2 + O → CO + O 2 , CO 2 + O 3 → CO + 2 O 2 , and CO 2 + CO 2 → CO + CO + O 2 . For these reactants, when UV is present in the system, the preferential reaction pathways were CO 2 + O → CO 3 , CO 2 + O 3 → CO 3 + O 2 , and CO 2 + CO 2 → CO 3 + CO.
It is worth noting that the preference of some reaction routes may not depend only on the type of ionizing incoming radiation but also on the local ice temperature and excess of energy available, initial molecular abundances, as well as ice temperature and porosity (these last three parameters influence the dielectric constant of ices; see additional discussion at Pilling et al. 2021).

Bimolecular reactions involving collisions with CO
The reaction scheme below presents selected bimolecular reactions involving collisions with CO with the calculated branching ratio for the experimental data employing CR and UV. From the considered reaction set, we observe that, depending on the type of energy input, the bimolecular reactions involving collisions with CO have different preferential reaction routes. For example, when CR is employed, the reactions between CO + C, CO + O, CO + O 2 , and CO + O 3 yield preferentially C 2 + O, C+O 2 , CO 2 + O, and CO 3 + O, respectively. For these reactions, when UV is employed, the products are C 2 O, CO 2 , CO 3 , and CO 2 + O 2 , respectively. The best-fit UV model shows that some bimolecular reactions, such as the reaction between CO + CO and CO + CO 2 , do not have any preference among the studied channels. Future works may focus on this issue, for example, with the addition of other reaction routes and molecular species into the chemical system.

Molecular Desorption from Ice to the Gas Phase Induced by Radiation
The best-fit model for the experimental data employing CRs gives a desorption yield (this case is also named the sputtering yield), which is in good agreement (about four times higher) with the value determined experimentally by Pilling et al. (2010a) of 2.2 × 10 4 molecules ion −1 employing an associative exponential equation to describe the experimental data. For the data employing UV photons (Martín-Doménech et al. 2015), the modeled desorption yield also presents good similarity (about only 2.5 times higher) with the photon-induced desorption values. For the sake of comparison, the authors of the UV experiment have determined that for CO (the most abundant daughter species), the desorption yield is 1.2 × 10 −2 molecules photon −1 , and for minor products such as O 3 or CO 3 , the value is on the order of 10 −3 -10 −4 molecules photon −1 , respectively. The reader should bear in mind that, in the current methodology, we adopt the experimental desorption values (without considering its uncertainty) in the DSC inside the SF. This also helps to minimize the degeneracy of the solutions during the navigation processes to find the best solution for the coupled chemical reaction set. Figure 5 presents the molecular desorption plots of best-fit models (top panels: CR model; bottom panels: UV models). The left panels present the desorption column density, in units of molecules per centimeter squared, for each modeled chemical species as a function of time. The right panels present the intrinsic molecular desorption rate, in units of hertz, for each molecular species. The desorption column density (DES(t)) and the intrinsic desorption rate (k des_i ) were described by Equation (2). As discussed previously, k des_i and the other reaction rates considered (listed in Table 2) are time independent in this manuscript. In these panels, we observe that, as the individual surface coverage and molecular abundances changes with time, the desorption column density also changes. However, curiously, we notice that depending on the time, the desorption might be ruled by one molecular species or another. For example, in the case of the CR experiment (panel (a1)), the desorption column density at the beginning is mainly ruled by the CO 2 and at the final is mainly due to the O atoms. However, in the case of the UV experiment (panel (b1)), we also observed CO 2 dominating the desorption column density at the beginning of the experiment but the desorption is mainly ruled by CO at the end. Interesting competitions between desorption are observed for the other species, and a discussion of this phenomena will be the subject of future study.
The intrinsic molecular desorption rate, presented in the right panels of Figure 5, shows several similarities between the two modeled experimental data sets. For example, CO 3 and C 2 O were among the species with the highest values in the two models, as CO 2 , O 3 , and C 2 O 3 were among the species with the lowest values. The larger difference occurs for O 2 , which, in the CR model, has a value of 0.01 s −1 and about 0.03 s −1 for the UV model. In general, the intrinsic desorption rate of the observed species, with the exception of CO, was larger in the UV model in comparison with the CR model. The calculated average intrinsic desorption rates were about one order of magnitude larger in the CR model than in the UV model (1.6 × 10 −3 s −1 and 1.6 × 10 −4 s −1 for the CR and UV models, respectively).
The current results on the molecular desorption induced by radiation may help us to understand the type and the quantity of freeze-out species from astrophysical ices induced by radiation in colder regions of space. These newly produced gas-phase molecules will join the gas-phase reaction network also enhancing the chemical complexity in the gas phase. Future manuscripts, employing the current methodology, will address ices with different temperatures and different compositions to provide a clear view of reaction routes (and physicochemical processes) within the ice bulks induced by incoming radiation in such colder astrophysical scenarios. Additionally, the next manuscript will also focus on the relationship between the intrinsic desorption and intermolecular binding energies (as well as the activation energies) to help to clarify these features.

Molecular Abundances at Chemical Equilibrium
As discussed by Pilling et al. (2019) and Rachid et al. (2020), the CE phase or stage occurs when the variation of summed chemical abundances is negligible (<1% of the total abundance). In this regime, the processing by incoming radiation induces equal reaction rates for the destruction and production of a given species, which keeps molecular abundances nearly constant even during continuously and extended further irradiation (see also Almeida et al. 2017;Bonfim et al. 2017;Rachid et al. 2017;Vasconcelos et al. 2017aVasconcelos et al. , 2017bVasconcelos et al. , 2017cCarvalho & Pilling 2020a, 2020bFreitas & Pilling 2020). Precisely speaking, due to continuous desorption-induced processes by the incoming radiation, the concentration of molecules in the ice (observed and nonobserved species) does not have an exactly constant value as the fluence increases; in fact, the column densities of all frozen species decrease slowly (with similar slope) after CE is reached. Additionally, the concept of CE can be also employed to estimate the type and number of specific molecules that go to the gas phase during desorption processes (usually considered to be the same proportion that exists in the ice at the CE stage). See additional discussions at Pilling et al. (2019) and Carvalho & Pilling (2020a, 2020b. Figure 6, left panels (panels (a1) and (b1) for CR and UV models, respectively), presents the molecular abundances in percentage at the CE phase (also called the CE branching ratio -EBR%) obtained from the best-fit models employing the current PROCODA code in the two experimental data sets (CO 2 ices at 13 K irradiated by Ni ions from Pilling et al. 2010a, and at 8 K irradiated by UV photons from Martín-Doménech et al. 2015). Observed and nonobserved (also called unknown) species are displayed in the bar plots by the blue and red colors, respectively. It is worth noting that the identification of nonobservable species was only possible with the current methodology, which solves the coupled system of chemical reactions. This achievement helps to clarify the experimental data of frozen samples in the IR. The percentage of nonobserved species at equilibrium chemistry in the CR experiment was almost one-third the value of the UV model. Considering the observed species, all of them presented larger molecular abundances in percentage at the UV model; however, the bigger difference was noticed for CO 3 (∼70 times higher), followed by O 3 , CO, and CO 2 (only ∼1.2 times higher). For the nonobserved species, we observe that the molecular abundances of O and C atoms in the CR model are much larger than those observed in the UV model, which indicates that the larger energy deposited by CR induces larger atomization in the sample. Curiously, the molecular abundances, in percentage, of the species O 2 , C 2 O, and C 2 O 3 were larger in the UV model than in CR model. Virtually no C 2 O 2 and C 2 O 3 are present at the equilibrium phase for the CR model. The molecular abundances at CE for the two studied ices are also listed in Table 3.
Since the employed experimental papers do not handle the nonobserved species, a direct comparison between the current molecular abundances at equilibrium (Figure 6) with experimental data is not possible. However, it is worth noting that in the case of the UV data from Martín-Doménech et al. (2015), the authors provide some upper limits for the nonobservable O 2 species (assuming also that all remaining O atoms are locked in O 2 molecules). Considering their upper limit abundance for O 2 , the estimated O 2 /O 3 ratio in their final UV fluence was roughly 7 (see their Figure 8). In the current manuscript, the calculated O 2 /O 3 ratio and (O+O 2 )/O 3 ratio for the UV model at CE were 2.6 and 3.2, respectively (see also Table 3).
The percentage molecular desorption at CE obtained from the best-fit models is shown in the right panels of Figure 6 (panels (a2) and (b2) for the CR and UV models, respectively). As discussed previously, the total desorption yields obtained in these models (see also Table 2) are close to the estimated values provided in the experimental papers. Some interesting results can be withdrawn from these figures; for example, the largest desorbed species in the CR model was the O atoms, while in the UV model, it was the CO molecule. In both models, the second larger desorbed species was the parent species CO 2 . In the CR model, other than O atoms, the C atoms also presented considerable desorption in percentage (the thirdlargest value). The percentage desorption of O 2 in the UV model was roughly 30 times larger than that observed in the CR model, which indicates that the irradiation of pure CO 2 ices by UV can also release (in addition to the larger amount of CO to the gas phase) a considerable amount of oxygen molecules to the gas phase.

Astrophysical Implications
To explain the wealth of complex organic molecules in hot cores and hot corinos observed with sensitive and resolved submillimeter observations (see review by Jørgensen et al. 2020), different astrochemical models employing three-phase chemical processes have been employed. Among the input parameters, these models adopt reaction rates in both the gas and solid phases, which are calculated from laboratory experiments. Nevertheless, the reaction rates in the solid phase represent a huge uncertainty source in the models because often they are guessed or assumed to be the same as gas-phase reactions, as discussed by Cuppen et al. (2017). Because of these issues, a dubious reputation has been given to surface chemistry, namely, "the last refuge of the scoundrel" (Charnley et al. 1992). However, substantial efforts have been made to improve the accuracy of the parameters used in chemical models (e.g., Cazaux et al. 2016;Minissale et al. 2016;Lamberts et al. 2017;Suhasaria et al. 2018;Qasim et al. 2020).
This paper contributes to reducing the inaccuracies of chemical reactions formed via photodissociation or radiolysis of the frozen CO 2 in astrochemical models. Carbon dioxide holds a large reservoir of depleted carbon in icy reservoirs located in star-forming regions. Quantifying its chemical pathways in the solid phase triggered by photons or CRs is essential to better understand the redistribution of carbon in the ISM. In particular, this work shows that IR inactive O 2 is formed from the carbon dioxide processing by photons (mainly) or CRs. O 2 was detected by Rosetta in the coma of comet 67P/Churyumov-Gerasimenko with abundances relative to water of around 3.8 ± 0.85 % (Bieler et al. 2015). In the 67P comet, the observations suggest that the O 2 signal might be With the current methodology, we estimate for the first time average values for the reaction rates in bimolecular and termolecular reactions within astrophysical ice analogs in the presence (and triggered) of ionizing radiation. Additionally, with the PROCODA code, we can have a better picture of the CE phase that happens at astrophysical ices when irradiated, at a constant temperature, for a long time in space. In addition to the characterization of molecular abundances at equilibrium, we estimate the desorption rates of frozen species to the gas phase in the presence of incoming radiation. These values might be employed in future astrochemical models to map chemical evolution embedded species in astrophysical environments under the presence of an ionizing radiation field.
The upcoming launching of the James Webb Space Telescope will improve the sensitivity and resolution by factors of 50 and 20, respectively, of ice bands in the protostellar spectra, compared with previous observations performed by Spitzer/IRS (Wright et al. 2015). To shrink the gap between experiments and models, this work provides an affordable methodology to derive physicochemical parameters associated with the chemical reactions of several molecules, ions, and radicals. In particular, we expect that the employment of the current code with ices irradiated with the same initial composition, irradiated by the same projectile but at different temperatures, may also help to determine other thermodynamical properties within the ices such as the activation barrier of reactions, entropy, enthalpy, and Gibbs energy. Future works will map the chemical changes of such ices for different initial compositions and ice porosity, which permits the determination of the dependence of k of a given reaction on the molecular environments (which could be written in terms of the grain compaction, types of neighbor species, or the average electric field of the medium, i.e., the dielectric constant). An application of PROCODA at the experimental data set containing ice heating only will allow us also to find dependencies of the reaction rate on the temperature.

Conclusions
We present a computational methodology (PROCODA code) to clarify the chemical evolution of ices investigated experimentally under photolysis/radiolysis processes until reaching the chemical equilibrium (CE) phase after receiving large radiation fluences. Here, we described the chemical evolution of two experimental studies involving the processing of pure CO 2 ices at low temperature, one employing a CR analog in 13 K ice (52 MeV Ni ions; Pilling et al. 2010a) and another employing UV in 8 K ice (10 eV; Martín-Doménech et al. 2015). The code was developed using the Python language and briefly consists of a mathematical process to solve a set of coupled differential equations that describe the chemical evolution as a function of time. For the current data sets, there were 11 different chemical species (four observed: CO 2 , CO, O 3 , and CO 3 ; seven nonobserved or unknown: O, O 2 , C, C 2 , C 2 O, C 2 O 2 , and C 2 O 3 ), 100 reaction routes (e.g., direct dissociation reactions, bimolecular and termolecular reactions, and desorption processes). The main conclusions were as follows: (i) With the employed methodology, we quantified not only the abundances of observable species (in the IR spectra irradiated ices) but also, and for the first time, the abundances of other nonobservable (unknown) species. Such characterization helps us understand the experimental data.
(ii) We characterized and quantified the radiation-induced desorption process, of both observed and nonobserved species, which includes the evolution of the individual desorption column density as a function of time and intrinsic desorption rates. The desorption rates (and desorption yield) from the best-fit models were 1.4 × 10 14 molecules s −1 (8.7 × 10 14 molecules ion −1 ) for the ice bombarded by CR, and 4.8 × 10 12 molecules s −1 (3.0 × 10 −2 molecules photon −1 ) for the ice irradiated with UV. Such values are in good agreement with the values listed by the authors in the employed experimental papers. The largest desorbed species (in percentage) at the CE phase were the O atoms and the CO molecule for the CR and UV models, respectively. In both models, the CO 2 was the second largest desorbed species. The percentage desorption of O 2 in the UV model was roughly 30 times larger than that observed in the CR model, which also indicates that the irradiation of pure CO 2 ices by UV can also release, other than the larger amount of CO to the gas phase, a considerable about of oxygen molecules to the gas phase. The calculated average intrinsic desorption rates were 1.6 × 10 −3 s −1 and 1.6 × 10 −4 s −1 for the CR and UV models, respectively.
(iii) The formation of O 2 after energetic processing of CO 2 ice suggests that oxygen molecules may not likely start out as O 2 , in the coma of 67P/C-G comet. It is conceivable that photolysis or radiolysis of molecules as carbon dioxide leads to its formation. Additionally, this work supports the conclusion that dicarbon trioxide (C 2 O 3 ) can be formed in the ice (mainly during UV irradiation) and break up in the comet coma upon sublimation, as an alternative route for the detected O 2 .
(iv) The best-fit models provide numerical values for the reaction rates for the considered chemical reactions within the ices, as well as the branching ratio for selected reaction routes. The average reaction rates for the direct dissociation (radiation-induced reactions) were =´k 3.6 10 d 3 and 5.0 × 10 −4 s −1 for the CR model and UV model, respectively. The average calculated dissociation cross sections were s =´-1.8 10 d 12 cm −2 and 2.5 × 10 −18 cm −2 for the CR model and UV model, respectively. These values are less than one order of magnitude lower than the values estimated from the experiments. For the bimolecular reactions within the irradiated ices, the average values for the reaction rates were =+ k 1.1 10 A B 25 and 5.3 × 10 −26 cm 3 molecule −1 s −1 for the CR model and UV model, respectively. In the case of the termolecular reactions, the average values for the reaction rates were + + k A B C = 1.1 × 10 −48 and 7.2 × 10 −49 cm 6 molecule −2 s −1 for the CR and UV models, respectively. These average values can be employed in future astrochemical models to map chemical evolution embedded species in astrophysical environments under the presence of an ionizing radiation field.
(v) The molecular abundances at the CE phase were characterized, considering both observed and nonobserved species. At this phase, the percentage of nonobserved species in the CR model was almost one-third the value for the UV model. All observed species in the UV model presented larger molecular abundances, in percentage, than in the CR model. For the nonobserved species, the best-fit model shows that the percentage molecular abundances of O and C atoms in the CR model are much larger than those observed in the UV model, which indicates that the larger energy deposited by CR induces larger atomization in the sample. Curiously, the percentage molecular abundances of O 2 , C 2 O, and C 2 O 3 were larger in the UV model than in the CR model.
The current methodology helps to better characterize the chemical evolution of ices during its processing by radiation and also can place important constraints on the astrochemical models. We hope this methodology may also help illuminate the molecular complexity of astrophysical ices that will be observed with high resolution by the forthcoming James Webb Space Telescope.

Data availability
The data underlying this article will be shared on reasonable request to the corresponding author.

Appendix A Interpolation of Input Data
The current work, we applied an interpolation (generated by employing a linear or/cubic-spline interpolation) to the input experimental data, which improves the behavior of optimization algorithm. The interpolation processed did not add or miss any structure on the modeled data.
In the current methodology, for each observed species in the CR model, 710 nonequally spaced interpolated data points were considered from 0 to 5000 s (∼1.4 hr). For each observed species in the UV model, 1010 nonequally spaced interpolated data points were considered from 0 to 14,382 s (∼4 hr). Figure A1 presents the original experimental data (black points) overlaid with its interpolated points (red points) employed in each model Figure A1. Comparison between experimental data (black points) and interpolated experimental data (red points) employed in the models. Panel (a) presents the data from pure CO 2 ice at 13 K irradiated by 52 MeV Ni ions (from Pilling et al. 2010a), and panel (b) presents data from pure CO 2 ice at 8 K irradiated by UV photons (from Martín-Doménech et al. 2015). Models Considering Different Input Parameters for the CR Data (Bad-fit Models) for Comparison with the Best-fit Model Figure B1 presents the evolution of column density and selected column mass obtained from eight models with different values as weights in the SF (bad-fit models), for comparison with the best-fit model presented in Figure 3 for CR data (pure CO 2 ice at 13 K irradiated by 52 MeV ions from Pilling et at. 2010a). Panels (a1)-(a8) indicate the evolution of column density, as a function of time, and panels (b1)-(b8) show the evolution of selected column mass (e.g., summed column mass of observed species, summed column mass of nonobserved species (unknown species), summed column mass of desorbed species, and total summed column mass). The experimental data are represented by filled black symbols (only the four observable species CO 2 , CO, O 3 , and CO 3 ), and the lines are the models. The two vertical segments presented in panels (b1)-(b8) (at larger times) quantify selected differences between summed column mass of observed and nonobserved species and are employed in the calculation of the slope similarity criterion to guarantee CE of the system at larger radiation fluences. In all panels, some important output parameters for the best-fit models are displayed in the headers.
We often observe in these models' large values for the summed chi-square of the observed species, as well as low values for the CE criterion (also called "slope similarity" in the panel headers). The parameters employed (weights in the SF) and the selected outputs of these models are presented in Table B1. This Appendix helps to clarify the importance of the SF (and its weights) in constraining and decreasing the degeneracy of solutions of the chemical system. Figure B1. Comparison between models considering different input parameters in the SF applied to the CR data. Panels (a1)-(a8) present the column density as a function of time, and panels (b1)-(b8) present the evolution of selected column mass as a function of time. Lines are the models, and symbols are the experimental data (pure CO 2 ice at 13 K irradiated by 52 MeV Ni ions from Pilling et al. 2010a). Model parameters are described in Table B1.    Notes. a Input parameters: p1 (SF's weight for the CO2), p2 (SF's weight for the CO), p3 (SF's weight for the O3), p4 (SF's weight for the CO3), p5 (SF's weight for the mass conservation), p6 (SF's weight for the desorption similarity), p7 (SF's weight for the slope similarity criterion; hypothesis for the CE stage), and pDES (estimated input for the desorption yield). b Additional input parameters for the CR Model are: N_CO2_initial (molecules cm −2 ) = 2.14 × 10 14 , phi (ions cm −2 s) = 2 × 10 9 , d_ice (cm) = 1 × 10 −5 , and s_ice (cm 2 ) = 0.8. c Additional input parameters for the UV Model are: N_CO2_initial (molecules cm −2 ) = 2.00 × 10 17 , phi (photons cm −2 s) = 2 × 10 14 , d_ice (cm) = 8.1 × 10 −6 , and s_ice (cm 2 ) = 0.8. d Outputs parameters: SF = Score function; sumCHI2 = summed chi-squared function for the observed species; TotDES = total summed desorption; MSS (%) = Mass similarity criterion (similarity between initial mass and final mass) in percentage; SSC (%) = Slope similarity criterion in percentage between summed observed mass and summed unknown mass (considering the last 200s); hypothesis for the CE stage; Average reaction rates: For direct (k d ), bimolecular (k A+B ) and termolecular (k A+B+C ) reactions. Other output parameters are: individual ks, individual molecular abundance at CE, individual desorption column density, and individual intrinsic molecular desorption in units of hertz. In this table we mark some very incorrect output values (much smaller or larger than values obtained in the best-fit model PROCODA V3b-2 CR) in bold and red color (not applied for SF).

Appendix C Relative Standard Error of Selected Output Parameters Considering 30 Different Calculations of the CR Model
To evaluate the uncertainty of the current methodology, we calculate the relative standard error (RSE), in percentage, of the reaction rates and other output parameters considering 30 different calculations of the CR Model (with the same initial input). The RSE is a measure that shows how large the standard error is relative to the estimated parameter value (here, it is the average value of the given parameter). In other words, it provides the level of fluctuation of selected output parameters with respect to its average value computed from the different models. It is calculated by dividing the standard error of an estimated value by the estimated value itself, and then multiplied by 100 and expressed as a percent. The RSE, in percentage, for a given parameter in terms of its average value is given by: 100% 5 x x Figure C1. RSEs for the obtained reaction rates considering 30 different calculations of the CR model with the same initial conditions. Panels (a), (b), (c), and (d) show the RSEs for rates in the direct reactions, in the bimolecular reactions, in the termolecular reaction, and in the intrinsic molecular desorption from ice to gas phase induced by radiation, respectively.
where s s x , , x x are the standard deviation, the standard error, and the average value of a given parameter x in the data set, respectively, and n is the number of values used for average calculation (in this case, n = 30). Figure C1 shows the RSE for the reaction rates considering 30 different calculations of the CR model. Panels (a), (b), (c), and (d) show the RSE for the rates of the direct reaction, bimolecular, reactions, termolecular reactions, and molecular desorption from ice phase to gas phase triggered by the incoming radiation, respectively. The horizontal dashed line indicates the average RSE for each shown subset. We observe for this model that the average standard error was about 6.4% for the direct reactions, with most values of RSE below 3%. For this reaction rate set, the largest errors found were for k20 (∼40%) and k91 (∼25%). The average RSE for the bimolecular reactions and termolecular reactions were 12.1%, and 16.9%, respectively. Despite most RSE values being below 10%, some rates presented much larger RSEs, such as k7 (∼45%), k11 (∼90%), k31 (∼40%), k35 (∼95%), k37 (∼60%), and k58 (∼40%). Finally, the average standard error for the intrinsic molecular desorption was 1.8% (the lowest relative error of the rate set). Considering all of these RSEs, we can consider, as a good estimate (within the safety margin), that the current methodology to calculate the solid-state reaction rates for the irradiated ices has an error below 20%. We expect similar behavior for the UV model. Figure C2 shows the RSEs for other selected outputs considering 30 different calculations of the CR model with the same initial conditions. Panel (a) shows the RSE for the molecular abundances at CE (both observed and nonobserved species) and, despite individual fluctuations, the average error was around 3.9%. Despite most of the modeled species having RSEs for the abundances at equilibrium below 2%, a larger error was observed in the calculation of the nonobserved species C 2 O 2 (∼20%). Therefore, we can consider once more, as a good estimate (within the safety margin), that the molecular abundances at CE obtained by the models also have RSEs below 20%. Panel (b) of Figure C2 shows the RSE for the SF value (0.4%), the chi-squared value from the observed species (3.2%), the mass similarity criterion (<0.1%), the summed desorption yield (0.5%), and the slope similarity criterion for the CE (<0.1%). The average error for these parameter sets was 0.8%. We expect similar behavior for the UV model.
Finally, it is worth noting that this appendix is an effort to quantify the convergence of the BFGS optimization algorithm for a fixed set of inputs (it is not based on the errors of the input experimental data or the details of the interpolation between).