# A fast neural emulator for interstellar chemistry

a Asensio Ramos, C Westendorp Plaza, D. Navarro-Almaid a, P.

Rivière-Marichalar, V. Wakelam, A. Fuente

► To cite this version:

a Asensio Ramos, C Westendorp Plaza, D. Navarro-Almaid a, P. Rivière-Marichalar, V. Wakelam, et al.. A fast neural emulator for interstellar chemistry. Monthly Notices of the Royal Astronomical Society, 2024, 531 (4), pp.4930-4943. 10.1093/mnras/stae1432. hal-04778219

HAL Id: hal-04778219

https://hal.science/hal-04778219v1

Submitted on 14 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# A fast neural emulator for interstellar chemistry

A. Asensio Ramos<sup>1,1,2*</sup> C. Westendorp Plaza,<sup>1,2</sup> D. Navarro-Almaida<sup>1,3</sup> P. Rivière-Marichalar<sup>1,4</sup>
V. Wakelam<sup>1,5</sup> and A. Fuente<sup>6</sup>

Accepted 2024 June 3. Received 2024 June 3; in original form 2024 March 12

## ABSTRACT

Astrochemical models are important tools to interpret observations of molecular and atomic species in different environments. However, these models are time-consuming, precluding a thorough exploration of the parameter space, leading to uncertainties and biased results. Using neural networks to simulate the behaviour of astrochemical models is a way to circumvent this problem, providing fast calculations that are based on real astrochemical models. In this paper, we present a fast neural emulator of the astrochemical code Nautilus based on conditional neural fields. The resulting model produces the abundance of 192 species for arbitrary times between 1 and $10^7$ yr. Uncertainties well below 0.2 dex are found for all species, while the computing time is of the order of $10^4$ smaller than Nautilus. This will open up the possibility of performing much more complex forward models to better understand the physical properties of the interstellar medium. As an example of the power of these models, we ran a feature importance analysis on the electron abundance predicted by Nautilus. We found that the electron density is coupled to the initial sulphur abundance in a low-density gas. Increasing the initial sulphur abundance from a depleted scenario to the cosmic abundance leads to an enhancement of an order of magnitude of the electron density. This enhancement can potentially influence the dynamics of the gas in star formation sites.

**Key words:** astrochemistry – methods: numerical – stars: formation – ISM: abundances – ISM: clouds – ISM: molecules.

## 1 INTRODUCTION

Interstellar matter is continuously processed describing a life cycle that connects stellar and interstellar environments. Interstellar matter forms filamentary clouds that fragment and collapse to form dense cores, where protostars are born. These protostars develop accretion discs, these latter being at the origin of planetary systems. The life cycle is closed when, at the end of their lives, stars inject a good fraction of their mass back to the interstellar medium, either through gentle winds, in the case of low-mass stars, or as supernova explosions, in the case of massive stars.

Molecules serve as excellent tracers of the physical and chemical evolution of interstellar matter (Jørgensen, Belloche & Garrod 2020). To date, more than 200 different molecules have been identified in space.<fnref n="1" /> The molecules observed in space go from the simplest hydrides, such as OH, $OH^+$, CH, and $CH^+$, present in diffuse interstellar clouds (Gerin, Neufeld & Goicoechea 2016; Gerin et al. 2019), to the complex organic molecules, some of them with prebiotic

relevance, that are detected in pre-stellar cores, young protostars, and planet-forming discs (Crockett et al. 2015; Rivière-Marichalar et al. 2020; Guélin & Cernicharo 2022; Öberg, Facchini & Anderson 2023). In environments that are heavily obscured at ultraviolet, optical, and near-infrared wavelengths, the observation of molecules becomes the only means to get access to the physical and chemical processes at work. Observations have demonstrated that the chemical nature of molecules varies drastically as objects evolve and their environment change. Some molecules are routinely used as chemical diagnostics of critical aspects such as, e.g. the gas ionization degree (see e.g. Caselli et al. 1998; Bergin et al. 1999; Fuente et al. 2019; Pineda et al. 2024), the local flux of cosmic rays (Favre et al. 2017; Cabedo et al. 2023), the gas kinetic temperature (Hacar, Bosman & van Dishoeck 2020), and the density (Navarro-Almaida et al. 2021). However, the reliability of molecular diagnostics used to determine the properties of the interstellar medium depends on our ability to understand the interstellar chemistry and the predictive power of our chemical models. Complex chemical models including thousands of reactions are routinely used in astrophysics (Le Petit et al. 2006; Taquet, Ceccarelli & Kahane 2012; Ruaud, Wakelam & Hersant 2016; Holdship et al. 2017; Laas & Caselli 2019). These codes are time-consuming and computationally expensive, which limits the total number of computed models. This limitation challenges the interpretation of the data provided by current telescopes, which

© 2024 The Author(s).

<footnote n="1"><sup>1</sup>Instituto de Astrofísica de Canarias (IAC), Avda Vía Láctea S/N, E-38200 La Laguna, Tenerife, Spain</footnote>

<footnote n="2"><sup>2</sup>Departamento de Astrofísica, Universidad de La Laguna, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="3"><sup>3</sup>Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="4"><sup>4</sup>Observatorio Astronómico Nacional (OAN,IGN), Calle Alfonso XII, 3. E-28014 Madrid, Spain</footnote>

<footnote n="5"><sup>5</sup>Laboratoire d'astrophysique de Bordeaux (LAB), Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, F-33615 Pessac, France</footnote>

<footnote n="6"><sup>6</sup>Centro de Astrobiología (CSIC/INTA), Ctra. de Torrejón a Ajalvir km 4, E-28850 Torrejón de Ardoz, Spain</footnote>

<footnote>* E-mail: aasensio@iac.es</footnote>

<footnote n="1">¹See the Cologne Database for Molecular Spectroscopy (CDMS; Müller et al. 2001) found at https://cdms.astro.uni-koeln.de/classic/molecules.</footnote>
allows to simultaneously observe tens, even hundreds, of molecular lines (Fuente et al. 2014, 2019; McGuire et al. 2020; Cernicharo et al. 2023), and map large areas in affordable amounts of time (Friesen et al. 2017; Gratier et al. 2021).

Molecules have an important influence on the heating and cooling of interstellar gas and, therefore, on the cloud structure of the interstellar medium in galaxies. Molecular clouds are the site of ongoing star formation to a large extent because molecules are efficient gas coolers, thereby diminishing thermal support relative to self-gravity. Moreover, molecules suppress the degree of ionization in these clouds and, therefore, decouple the gas from any supporting magnetic fields. Molecules, thus, play a key role in the formation process of stars and planets. Steady state models have been widely used to interpret molecular observations in the last decades. Many of the processes linked to the formation of a star, such as bipolar outflows and disc formation, have time-scales ($10^3-10^5$ yr) smaller than the time required for most species to reach a steady state in their abundances (~1 Myr). Coupling dynamics and chemistry is therefore of paramount importance to progress in our understanding of the physical and chemical evolution from the initial cold pre-stellar core to a protoplanetary disc, where planets are formed.

Numerical simulations have become a valuable tool to study, in a realistic manner, the evolution of interstellar matter during the star formation process. This dynamical evolution has been investigated following different approaches: analytical one-dimensional models (see, e.g. Priestley, Viti & Williams 2018), two-dimensional semi-analytical models (Drozdovskaya et al. 2014, 2016), and three-dimensional radiation-magnetohydrodynamical (R-MHD) simulations (see e.g. Commerçon et al. 2012; Hincelin et al. 2016; Masson et al. 2016; Zhao et al. 2018a; Gong et al. 2020; Hennebelle et al. 2020; Gong et al. 2021; Lebreuilly et al. 2024). MHD simulations have become the state-of-the-art approach for the analysis of the gravitational collapse of pre-stellar cores into protostellar objects probing the decisive role of the magnetic field in the process of collapse, fragmentation, and disc formation. A few of these state-of-the-art simulations already incorporate a reduced chemical network that provides important parameters for MHD calculations, such as the gas ionization degree (Marchand et al. 2016; Zhao, Caselli & Li 2018b). However, full coupling of the MHD evolution of the gas with a full chemical reaction network with thousands of reactions is challenging with the current computational capabilities (Hsu et al. 2021, 2023). As a feasible approximation, an alternative method based on the usage of virtual particles has been used. These particles store the key parameters required for the computation of the chemistry (namely temperature and density) and a full state-of-the-art chemical model is then used to predict their chemical composition as a function of time (Hincelin et al. 2016; Coutens, Commerçon & Wakelam 2020; Navarro-Almaida et al. 2024). Yet, these calculations are expensive in terms of computing time since a large number of particles is needed to sample all the protostellar components and only a handful of cases have been completed, thus far. Consequently, a full comparison of these 3D chemo-MHD simulations with observations has not been done yet, so that their predictive power remains largely unexplored.

Astrochemistry has become a necessary tool for understanding the interstellar medium. Chemical models consist of a system of coupled ordinary differential equations (ODE) that represent a chemical network. This system is resolved for a set of input parameters which include the gas physical conditions (density, temperature), cosmic-ray ionization rate per H atom, incident UV field, visual extinction, grain properties, initial chemical composition, and elemental

abundances. Computing large grids of complex chemical models, varying input parameters to describe all astrophysical environments, and predicting the chemistry for each time-step is unaffordable. It is therefore interesting to develop fast chemistry emulators that allow to address some problems in a fast way, which is precisely the purpose of this work. We are not the first team to explore emulators in chemistry. The vast majority of accurate enough emulators are based on neural networks, although other methods like random forest regression have been applied (Bron et al. 2021). de Mijolla et al. (2019) used a neural network to emulate astrochemical models in equilibrium after $10^7$ yr obtained with UCLCHEM (Holdship et al. 2017). Holdship et al. (2021) trained a neural network to predict the change in abundance during a small time-step, allowing the user to evolve the chemistry for long times by iteratively applying the neural network. Grassi et al. (2022) and Sulzer & Buck (2023) used auto-encoders to project large chemical networks into a low-dimensional latent space, where the chemistry is potentially evolved faster. A similar approach is followed by Maes et al. (2024) for emulating the chemistry in dynamical environments. Branca & Pallottini (2023) proposed neural networks as *ansatzs* for the time evolution of the chemistry and trained them using the ODEs in the loss function. Along similar lines, Branca & Pallottini (2024) used neural operators to emulate the chemical network (Lu et al. 2021).

Despite the amount of work done in the field, we think our architectural decisions, training strategy, and size of the chemical network can be of impact in the community. Moreover, we have trained this emulator to explore the impact of the sulphur elemental abundance on the physics and chemistry of the interstellar medium. In this paper, we focus on the gas ionization degree, an essential parameter for the dynamics of interstellar clouds.

## 2 CHEMICAL MODEL

We performed the abundance calculations using Nautilus 1.1 (Ruaud et al. 2016), a three-phase model, in which gas, grain surface and grain mantle phases, and their interactions, are considered. It solves the kinetic equations for the gas-phase and the solid species at the surface of interstellar dust grains. Species from the gas-phase can physisorb at the surface of dust grains with an energy that depends on the binding energy of surface species (Wakelam et al. 2017). They can be released to the gas phase by thermal desorption (Hasegawa, Herbst & Leung 1992), grain heating induced by cosmic rays (following Hasegawa & Herbst 1993, impact of UV photons (photodesorption, Ruaud et al. 2016), exothermicity of surface reactions (chemical desorption, Minissale et al. 2016), and sputtering by cosmic rays (following Wakelam et al. 2021). Detailed descriptions of the non-thermal desorption processes included in Nautilus can be found in Wakelam et al. (2021). The ice surface phase is composed of the first few monolayers of species on top of the grains (four in our case), while the rest of the molecules below these surface layers represents the bulk of the ice. The refractory part of the grains (below the bulk) is chemically inactive. The species on the surface can desorb into the gas phase, contrary to the species in the bulk. Only sputtering by cosmic rays can directly desorb species from the bulk ice. The version used in this paper allows to calculate the dust temperature using the expression of Hocuk et al. (2017) and the cosmic-ray ionization rate per H$_2$ molecule ($\zeta$) as a function of $A_V$ following the expression of Neufeld & Wolfire (2017) as described in Clément et al. (2023). However, we have disabled these options to create our data base since we were interested in determining the variation of the chemical abundances with these parameters.
## 3 NEURAL EMULATOR

Learning emulators for computationally heavy numerical procedures via neural networks is a well-established technique in the field of machine learning. The idea exploits the well-known fact that neural networks are universal function approximators (Cybenko 1988) to train them to learn the mapping between the input (physical parameters) and the output (the result of the numerical procedure). In theory, they can learn the mapping between input and output spaces with arbitrary precision provided that the network is large enough and the training set is sufficiently rich. In practice, this is often limited and only a limited precision is reachable.

In our model, the solution to the kinetic equations gives the abundance of a set of 192 species as a function of time $t$. The abundance is also dependent on the external physical conditions, that are characterized with the following parameters: the temperature ($T$), the initial gas density ($n_H$), the cosmic ray ionization rate ($\zeta_{H_2}$), the sulphur abundance ($[S]$) and the scaling factor for the surrounding UV flux ($\chi$). As a consequence, our aim is to end up with a neural network $A(t|M, \theta)$ that approximates the abundance for molecule M at time $t$, given the physical conditions summarized in the vector $\theta = (T, n_H, \zeta_{H_2}, [S], \chi)$.

### 3.1 Architecture

A crucial ingredient to end up with a successful emulator is to find a suitable architecture. We propose to emulate the interstellar chemistry using a neural field (NF), also known as implicit neural representations (INR) or coordinate-based representations (CBR) (see Mildenhall et al. 2020; Jarolim et al. 2024; Jarolim et al. 2022; Asensio Ramos 2023). These are fully connected (FC) neural networks that are used to map coordinates on the space (or space-time) to coordinate-dependent field quantities. NFs have desirable properties for serving as emulators over other parametrizations. They are very efficient in terms of the number of free parameters. Most important, they produce continuous (and differentiable$^2$) fields, which can then be evaluated at arbitrary inputs. Finally, they have a strong implicit bias, favouring specific signals.

In our case, we use an NF in which the input is the time at which we want to predict the abundance of a specific species and the output is the abundance of that species. Obviously, the NF needs to be conditioned on the physical parameters of the cloud and on the specific species. We show later how we do this conditioning technically. A typical NF gives the abundance at time $t$ of a certain species M and for a set of physical conditions $\theta$ by the following function composition:

$$A(t|M, \theta) = \phi_N^{M,\theta} \circ \dots \circ \phi_1^{M,\theta} \circ \phi_0^{M,\theta}(t), \quad (1)$$

where we explicitly state that each layer of the NF is conditioned on the vector $\theta$, that contains information about the physical conditions and the species, and the specific species to be predicted, M. Note that $(f \circ g)(\mathbf{x}) = f(g(\mathbf{x}))$ refers to the function composition operation. Typically, each layer of the NF is given by:

$$\phi_i^{M,\theta}(\mathbf{x}) = \sigma(\mathbf{W}_i \mathbf{x} + \mathbf{b}_i), \quad (2)$$

where $\mathbf{W}_i$ and $\mathbf{b}_i$ are weight matrices and bias terms, respectively, while $\sigma$ is a nonlinear activation function. In our case, we choose the rectified linear unit (ReLU; Nair & Hinton 2010) as the activation function.

The conditioning is done via two FC neural networks, one that deals with the physical parameters and the other one with the species. The input of the first neural network is the vector $\theta$, all of them properly normalized as described in Section 3.2. The input of the second neural network is a representation of the molecule in a 1-hot encoding$^3$ (i.e. a vector of dimension equal to the number of species, with a 1 in the component corresponding to the species and 0 in the rest). The outputs of each one of these neural networks are two scaling vectors, $\gamma_{\text{phys}}^{\theta}$ and $\gamma_{\text{sp}}^{M}$, and two bias vectors, $\beta_{\text{phys}}^{\theta}$ and $\beta_{\text{sp}}^{M}$. We follow the ideas behind feature-wise linear modulation (FiLM; Perez et al. 2017), so that each layer of the NF is modified using the previous vectors according to:

$$\phi_i^{M,\theta}(\mathbf{x}) = \sigma [\gamma_{\text{phys}}^{M,\theta}(\mathbf{W}_i\mathbf{x} + \mathbf{b}_i) + \beta_{\text{sp}}^{M,\theta}], \quad (3)$$

where $\gamma_{\text{phys}}^{M,\theta} = \gamma_{\text{phys}}^{\theta} \oplus \gamma_{\text{sp}}^{M}$ and $\beta_{\text{phys}}^{M,\theta} = \beta_{\text{phys}}^{\theta} \oplus \beta_{\text{sp}}^{M}$ are obtained via concatenation ($\oplus$ is the vector concatenation operator). One of the advantages of the FiLM approach to conditioning is that it is very easy to implement and it does not require a significant modification of the NF architecture. Additionally, all layers of the NF receive conditioning information, which facilitates the adaptation of the NF to the complex dependence of the output on the physical conditions.

NFs, with the previous definitions, are known to suffer from the so-called spectral bias (Rahaman et al. 2019; Wang, Wang & Perdikaris 2021), which prevents them from learning high-frequency functions. We follow Tancik et al. (2020) and solve this problem by first passing the time coordinate through a Fourier feature mapping that increases the dimensionality from 1 to $2m$ (with $m$ the number of Fourier features):

$$A(t|M, \theta) = \phi_N^{M,\theta} \circ \dots \circ \phi_1^{M,\theta} \circ \phi_0^{M,\theta}(\gamma(t)), \quad (4)$$

where

$$\gamma(t) = [t, \cos(2\pi \mathbf{B}_t t), \sin(2\pi \mathbf{B}_t t)]^{\dagger}. \quad (5)$$

We also include the time coordinate itself, to allow for outputs with very low-frequency time dependencies. The matrix elements of the matrix $\mathbf{B} \in \mathbb{R}^{m \times 1}$ are randomly sampled from a normal distribution with zero mean and standard deviation $\sigma_F$. After a non-exhaustive parameter search, we choose $\sigma_F = 2$ and $m = 512$, which produces a Fourier feature mapping of dimensionality 1025.

Although we have found success with adding Fourier features and using ReLU activation functions, other solutions like the SIRENs architecture (Sinusoidal Representation Networks; Sitzmann et al. 2020), which uses periodic activation functions, can also be found in the literature and could be explored for this problem in the future.

The graphical representation of the final model is depicted in Fig. 1. The two conditioning neural networks contain three hidden layers with 512 neurons. They produce the two conditioning vectors $\gamma_i$ and $\beta_i$ with 256 elements each. After concatenating them, they are used in the NF emulator to condition every hidden layer of the NF with 512 neurons. The final output is a single number, the abundance of the molecule for the specified physical conditions. The total number of trainable parameter is 4.83M, with 2.62M associated with the NF, and 1.05M and 1.15M associated with the conditioning networks for the physical conditions and the species, respectively.

Thanks to the interpolation capabilities of the NF, we can predict the abundance of any molecule at any arbitrary time between 1 and $10^7$ yr. Additionally, given that the NF is subdifferentiable, one can seamlessly compute derivatives of the abundance of any arbitrary

<footnote>$^2$Or, more precisely, subdifferentiable in case of using activation functions that are not differentiable everywhere.</footnote>

<footnote>$^3$Another option is to directly learn vector embeddings for all species, instead of them being computed from 1-hot vectors with the aid of neural network.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Architecture used in the emulator. The 'Emulator' block encapsulates the NF that predicts the abundance of a certain species at a given time. The two conditioning blocks produce vectors that are concatenated together and used in the neural field. Yellow blocks are input layers. Orange blocks are fully connected layers, with the number of neurons indicated for clarity. Light purple colour refers to the Fourier encoding layer explained in equation (5). The conditioning vectors enter into the fully connected layers following the FiLM paradigm.</figcaption></figure>

species with respect to both the physical conditions and time. This can be of great potential interest for the analysis of the chemistry and also for using the chemistry emulator in other downstream tasks.

## 3.2 Training set

The training set is built by solving the interstellar chemistry using Nautilus for a set of physical conditions and a total of 192 species. We add atomic and molecular species, as well as ions and ices. The chemistry network is solved for 64 times equispaced in log space between 1 and $10^7$ yr. A total of 150k models are run for building the training set. Roughly 16 per cent of these models do not pass a hard timeout threshold of 5 min and are then discarded<fnref n="4" /> This leaves us with 126k models, which is still a sufficiently large amount of models for training. The total number of training examples is, therefore ~126k×192, which amounts to 24.2M. The space covering of the training models is depicted in Fig. 2, where we witness a region of a reduced number of models for low densities and high temperatures, where Nautilus has problems converging the kinetic equations.

The physical conditions are selected randomly in the intervals defined in Table 1. The temperature is sampled uniformly without transformation, while the rest of the parameters are sampled uniformly in log space given that they are positive quantities and span several orders of magnitude. A subset of 10 per cent of the models is left out of the training set and used as a validation set. A final set of 5k models is obtained for testing purposes and it is used to compute the final performance of the emulator. We point out that other methods to

sample the physical conditions, such as Latin hypercube sampling, could have been used, but we have found that random sampling is sufficient for this problem. Additionally, it allows us to seamlessly extend the training set outside the currently considered regions.

The input physical conditions (in linear or log scales) are normalized to the interval $[-1, 1]$ when used as input to the conditioning neural networks. Likewise, $t$ is also normalized to the same interval after being transformed to log scale. The output abundances are similarly treated in log scale. However, since a large variability is found in the abundances, we normalize each time bin and molecule separately. This normalization is stored and applied in reverse when the model is used in production.

## 3.3 Training

The training proceeds by minimizing the mean-squared error between the output of the NF and the abundance for all considered molecules. The initial learning rate is set to $3 \times 10^{-4}$ and it is slowly reduced to $3 \times 10^{-5}$ using a cosine annealing schedule (Loshchilov & Hutter 2017). The optimization is carried out using the Adam optimizer (Kingma & Ba 2014) with a batch size of 1024 for 100 epochs. We utilize an NVIDIA GPU 4090 with 24 GB of memory. Given that the training takes around 24 h, we did not do an exhaustive hyperparameter search. The evolution of the training and validation losses during training are displayed in Fig. 3, showing no evidence of overtraining. We found the saturation of the validation loss but with a tiny difference with respect to the training loss. Although after epoch ~40 the decrease in the validation loss is relatively slow, we preferred to continue training for 100 epochs. The reason is that the diagnostics shown in the following sections carried out during intermediate epochs showed that the training was still improving.

In order to adapt the complexity of the NF to the specific problem, we opt for starting the optimization taking only into account the low-frequency components in the Fourier feature mapping, while gradually turning on higher frequency components during training. To this end, we extend the approach used by Lin et al. (2021) as follows. We reorder the rows of the matrix **B** in terms of their $\ell_2$ norm, with low frequencies appearing first. We then weight the Fourier features as follows:

$$ \gamma(t) = [t, w(\alpha) \cos(2\pi B_t), w(\alpha) \sin(2\pi B_t)]^\dagger, \quad (6) $$

where the $k$-th element ($k = 0, 1, \dots, m-1$) of the weighting function is given by

$$ w_k(\alpha) = \begin{cases} 0 & \text{if } \alpha < k \\ \frac{1-\cos((\alpha-k)\pi)}{2} & \text{if } 0 \le \alpha - k < 1 \\ 1 & \text{if } \alpha - k \ge 1, \end{cases} \quad (7) $$

and $\alpha \in [0, m]$ is slowly increased during training from 0 to $m$. We choose to linearly increase $\alpha$ from 0 to $m$ during the first 25 000 batches, and set $\alpha = m$ (no weighting) for the rest of the training. This procedure allows the NF to start learning low frequencies first and then gradually incorporate high frequencies. We have found that this procedure gives more consistent and faster results during training than simply starting with all frequencies active in the Fourier encoding from the beginning.

# 4 VALIDATION

Once the emulator is trained, we evaluate the performance using the 5k Nautilus models computed specifically for this purpose. Fig. 4 shows the time evolution of five randomly selected models. We only display a subset of species of certain relevance, although the output

<footnote n="4">⁴This hard threshold is applied to allow computing the data base in a reasonable time. Even with this threshold, the calculation requires ~4 h in 512 cores.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Two-dimensional histogram showing the covering of all 126k models of the training set shown as combinations of temperatures and the rest of parameters.</figcaption></figure>

<table><caption>**Table 1.** Physical conditions considered in the emulator, together with the range of values and the transformation used for sampling them.</caption><thead><tr><th>Parameter</th><th>Minimum</th><th>Maximum</th><th>Sampling</th></tr></thead><tbody><tr><td>t [yr]</td><td>1</td><td>10<sup>7</sup></td><td>log</td></tr><tr><td>T [K]</td><td>10</td><td>80</td><td>linear</td></tr><tr><td>n<sub>H</sub> [cm<sup>-3</sup>]</td><td>10<sup>4</sup></td><td>10<sup>7</sup></td><td>log</td></tr><tr><td>ζ<sub>H</sub> [s<sup>-1</sup>]</td><td>10<sup>-17</sup></td><td>10<sup>-15</sup></td><td>log</td></tr><tr><td>[S] [cm<sup>-3</sup>]</td><td>7.5 × 10<sup>-8</sup></td><td>1.5 × 10<sup>-5</sup></td><td>log</td></tr><tr><td>χ</td><td>0.1</td><td>10<sup>4</sup></td><td>log</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Evolution of the training and validation losses during training. An exponential average with a window size equal to the number of batches in an epoch has been applied.</figcaption></figure>

is obtained for all 192 species. The time evolution of the abundance computed with the neural emulator is displayed with solid lines, while the time evolution computed with Nautilus is shown with dashed lines. A simple visual comparison demonstrates that the emulator is able to correctly reproduce the time evolution of these species.

The neural emulator provides a significant gain in computing time. A typical Nautilus model requires of the order of 120–240 s to produce the output for all considered molecules at 64 time-steps equispaced in logarithm between 1 and $10^7$ yr. This number fluctuates depending on the physical conditions, which change the stiffness of the kinetic equations and impacts the total computing time. Our emulator is able to compute 4.8 models/s in an Intel Xeon CPU E5-2660@2.60 GHz, a rather standard desktop CPU from 2017. We clarify that a ‘model’ evaluation for the emulator requires using the emulator for all 192 species at 64 time-steps. This number goes up to

87 models s⁻¹ when using a GPU if models are computed one by one. This does not saturate the computing power of the GPU, which are specially tailored to compute many models in parallel. This saturation happens, in our specific conditions, when ~32 models are evaluated simultaneously. At this rate, we are able to compute 146 models s⁻¹. Therefore, our neural emulator is between $10^4$ and $3 \times 10^4$ times faster than Nautilus. Additionally, the computing time is not sensitive to the specific physical conditions, carrying out exactly the same number of operations for the whole parametric volume defined in Table 1.

An arguably better representation of the performance of the emulator is given in Fig. 5, together with Figs A1–A3 in Appendix A. We display the distribution of differences in logarithm between the abundance predicted by the emulator and the original one computed with Nautilus. For simplicity, we characterize the distributions showing the percentiles. The black solid line shows the median difference, which is interestingly close to zero for all species. This means that deviations above and beyond are equally likely. This is confirmed with the shaded areas, which display the range between percentiles 1 and 99 in blue, 5 and 95 in orange, and 16 and 84 in green. The figure shows that the difference lies below 0.2 dex in the majority of species of relevance for all times with 98 per cent probability. Only for a few species (like H₂S₃, HC₇N,...) this difference goes up to 0.4 dex for the latest phases of the chemical evolution. Note that we also show species in the ices (characterized with the letter J at the beginning when the species is in the surface, and K for those in the ice bulk).

The difference is consistently very small for the initial phases of the evolution for all species. The reason is that the initial abundance for all species is fixed for all models so it is easy to be learned by the emulator. The specific final abundance depends on the time evolution of the chemistry, which turns out to be more difficult to learn. The abundance at very large times converges to the equilibrium value, which should only depend on the physical conditions. Despite this fact, we find that the uncertainty in the emulator is often larger in the final stages of the evolution, although still consistently below 0.4 dex. This is because the equilibrium value is often reached at very large times so the uncertainty in the emulator is larger.

## 5 IMPORTANCE OF THE IONIZATION DEGREE AND INITIAL SULPHUR ABUNDANCE

Neural networks are often described as black boxes in the sense that, even though they can approximate the results well, the structure of the neural network itself generally does not provide any insight into the model being approximated. It is a reality (although the machine
<figure><img src="image_5.png" /><figcaption>Figure 4. Five samples per species of those considered of more interest. We show the time evolution of the abundance computed with the neural emulator in solid line, while the original time evolution of the test set is shown in dashed lines. These have been selected randomly with no cherry-picking.</figcaption></figure>

learning community is actively working towards solving this issue) that it is hard to understand the learned relationship between the input parameters and the outcome of the model. This information is sometimes of great interest in astrochemistry, as questions like what physical conditions favour or hinder the production of a given chemical species are frequently faced.

A potential avenue to understand the relationship between inputs and outputs in astrochemical models was addressed in Heyl, But- terworth & Viti (2023) using Shapley values (Shapley 1953). The Shapley value is an approach from game theory that provides a way to fairly distribute the outcome of a cooperative game among the players. In the context of astrochemical emulation, the chemical emulator is the cooperative game, the outcome (or prediction) is the abundance of a given chemical species, and the players, commonly referred as features, are the input parameters of the emulator. Shapley values therefore allow us to quantify the importance of the features for a given prediction. In the following, we give a brief description of the method applied to the chemical emulator. For a full description of the Shapley value approach, we refer to Shapley (1953).

Let us define $N = \{T, n_H, \zeta_{H_2}, [S], \chi\}$ as the full set of input parameters of the neural emulator and $\theta^j$ as the vector of input

parameters of training example j. Under the Shapley value framework (Lundberg & Lee 2017; Heyl et al. 2023), the explanation of the prediction $f(\theta^j)$ is assumed to be given by the following linear model:

$$f(\theta^j) \sim \langle f \rangle + \sum_{i=1}^{|N|} \phi_i^j, \qquad (8)$$

where we made the notational simplification $f(\theta^j) = A(t|M, \theta^j)$ once t and M are fixed. Put in words, the simplified explanation model is built by adding $\langle f \rangle$, the average value of the prediction, together with the Shapley values $\phi_i^j$ associated with all the ith input parameter at the jth training example. This linear model allows one to extract the importance of each input parameter for a given prediction by simply comparing the Shapley values.

Shapley values can be obtained from the predictions of our model. Let $S \subseteq N$ be a subset of N (e.g. $S = \{T, n_H\}$). The vector $\theta_S^j$ is obtained by keeping the elements from the subset S as they are and setting the rest of parameters to a random value extracted from the data set. With these definitions, the Shapley value $\phi_i^j$ associated with the input parameter i and training example j is the weighted average of the marginal contribution of the feature i over all possible coalitions
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Difference in logarithm between the emulator and the original run for selected species. The solid line shows the median difference, while the shaded area shows the quantiles 1–99 in blue (approximately 3σ), 5–95 (approximately 2σ) in orange, and 16–84 (approximately 1σ) in green.</figcaption></figure>

of features excluding the feature *i*:

$$ \phi_i^j = \sum_{S \in N \setminus \{i\}} \frac{|S|! (|N| - |S| - 1)!}{|N|!} \left[ f(\theta_{S \cup \{i\}}^j) - f(\theta_S^j) \right]. \quad (9) $$

As indicated, the summation is carried out for all subsets that do not include the input parameter of interest. Note that $|N|$ and $|S|$ are the number of elements in sets *N* and *S*, respectively.

Once the Shapley values are computed, the global importance of feature *i* in the prediction, denoted as $I_i$, is given by the mean absolute value of the Shapley values for that feature over the entire data set:

$$ I_i = \frac{1}{d} \sum_{j=1}^{d} |\phi_i^j|, \quad (10) $$

where *d* is the number of examples in the data set. Obviously, the relative importance of feature *i* can be obtained as:

$$ \bar{I}_i = \frac{I_i}{\sum_j I_j}. \quad (11) $$

Sulphur atoms, with an ionization potential of 10.36 eV, are easily ionized in translucent gas, becoming an important donor of electrons in these areas inside molecular clouds (see, e.g. Goicoechea et al. 2006; Fuente et al. 2019; Bulut et al. 2021). We expect that in denser

and shielded gas, where sulphur is mostly found in sulphur-bearing molecules, the initial sulphur abundance becomes less important to set the ionization degree of the gas. The ionization degree is a key quantity that influences the coupling between the gas and the interstellar magnetic field. Consequently, it strongly impacts the dynamics of the gas and, potentially, the star formation process. With the advent of numerical simulations of core collapse, this quantity is present in diffusive phenomena like ambipolar diffusion, Ohmic dissipation, and Hall effect (Masson et al. 2012). These diffusive phenomena were proposed as a solution to the magnetic braking catastrophe, in which the implementation of ideal magnetohydrodynamics (MHD) prevents discs from being formed (see, e.g. Allen, Li & Shu 2003; Galli et al. 2006; Hennebelle & Fromang 2008; Commerçon et al. 2010). It is therefore interesting to examine the impact of the initial sulphur abundance on the number of electrons for physical conditions corresponding to different layers of a molecular cloud and discuss the possible implications in the star formation process.

To investigate the impact of the model features on the number of electrons, we used the Python package SHAP (Lundberg & Lee 2017). We first created a wrapper around the emulator so that the input features are scaled to the range [0,1] and the output is the decimal logarithm of the chemical abundances at a given time instant. Unless specified, the values of the output are taken at 1 Myr.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Plots summarizing the SHAP values and feature importance for the 2048 predictions of the number of electrons. Left: feature importance as defined in equation (10) for each feature of the emulator. Right: a so-called beeswarm plot where dots show the SHAP values of each prediction for all features, coloured depending on their feature value.</figcaption></figure>

We performed a random uniform sampling of the parameter space with 2048 predictions. The SHAP values of the five features for the 2048 predictions and 192 chemical species were calculated with the Explainer object of the SHAP library. The importance of the features (see equation 10) setting the number of electrons is shown in Fig. 6. According to the ranking of features (left-hand side of Fig. 6), the density is the feature that makes the greatest impact on the number of electrons, with a contribution of 51.2 per cent (equation (11)), followed by the cosmic ray ionization rate (21.6 per cent), the FUV field strength (9.6 per cent), and the initial sulphur abundance and gas temperature, both with a contribution of 8.8 per cent. In the right-hand side of Fig. 6 we show the beeswarm plot of the ionization degree, that is, the SHAP values for all the features that correspond to each one of the 2048 randomly sampled predictions, colour-coded by the feature value. Our results indicate that SHAP values for density are anticorrelated with the feature (density) value, that is, high gas densities have a negative impact on the number of electrons. This is in agreement with our previous discussion: as density increases, sulphur is no longer in its ionized form but found in sulphur bearing molecules, thus reducing its contribution to the ionization state of the gas. For the rest of the features, this behaviour is inverted, as their SHAP values are correlated with their value. The cosmic-ray ionization rate and the FUV field strength have a positive correlation with their corresponding SHAP values. As expected, stronger ionizing agents have a positive impact in the ionization of the gas. Likewise, higher temperatures are associated with a positive impact on the number of electrons of the gas, although this correlation is weaker.

Correlations are quantified in Fig. 7. A strong correlation, $\rho = 0.99$, is seen between the SHAP and feature values of the cosmic-ray ionization rate, while a strong anticorrelation, $\rho = -0.97$, is present between the SHAP and feature values of the gas density. Correlations on the rest of features are weaker, meaning that predictions with a similar feature value may have different SHAP values for that feature. This indicates that the SHAP value for a feature may not be only dependent on that feature but on several interacting features. We investigate feature interaction in Fig. 7, showing the pairs of SHAP and feature values for a given feature, colour-coded by the rest of features. Gas density $n_H$ shows a clear interaction with the initial sulphur abundance [S] and the FUV field strength $\chi$ (see plots in the second column, fourth and fifth row, respectively, in Fig. 7). In

the case of the initial sulphur abundance, for any given value of this feature, lower densities yield higher (in absolute value) SHAP values for [S] than higher gas densities. Since SHAP values tell us how each feature contributes to the deviation of a given prediction from the average outcome (equation (8)), this indicates that the number of electrons is susceptible to change more when the initial sulphur abundance is modified in low-density gas than when it is modified in denser gas, where the number of electrons is more robust against changes in [S]. A similar behaviour is also observed in the interaction between density $n_H$ and the FUV field strength $\chi$.

We can provide a quantitative estimation of how different densities alter the variations in the ionization degree when the initial sulphur abundance is modified. We first considered the electron density $n_e$ in two low-density scenarios, with low and high initial sulphur abundances respectively. The family of vectors of input parameters with low density and low initial sulphur abundance, both normalized to the range [0,1], following the standard ordering of features given by the set $N = \{T, n_H, \zeta_{H_2}, [S], \chi\}$, is given by

$$ \theta^a = [T, 0, \zeta_{H_2}, 0, \chi] \quad (12) $$

with temperature $T$, cosmic-ray ionization rate $\zeta_{H_2}$, and FUV field strength $\chi$ set as free parameters. Similarly, the family of vectors of input parameters with low density and high initial sulphur abundance is given by

$$ \theta^a = [T, 0, \zeta_{H_2}, 1, \chi]. \quad (13) $$

We estimated the electron density in each case by averaging over the free parameters using a set of 30 000 uniformly sampled vectors of input parameters following equations (12) and (13). The resulting electron densities are:

$$ \langle \log n_e(\theta^a) \rangle = -6.53, $$

$$ \langle \log n_e(\theta^b) \rangle = -5.50, $$

yielding the following relative difference in the electron density for the low density case

$$ \frac{\langle n_e(\theta^b) \rangle}{\langle n_e(\theta^a) \rangle} = 10^{-5.50+6.53} \sim 10. \quad (14) $$

That is, in the low-density scenario ($n_H \sim 10^4 \text{ cm}^{-3}$), increasing the initial sulphur abundance from $[S] = 8 \times 10^{-8}$ to $[S] = 1.5 \times 10^{-5}$
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Correlations between SHAP and feature values of the 2048 predictions of the number of electrons for gas temperature (first row), gas density (second row), cosmic-ray ionization rate (third row), initial sulphur abundance (fourth row), and FUV field strength (fifth row). In each row, the 2048 predictions are colour-coded by the rest of the features. The histogram on the x-axis shows the uniform sampling in the range of feature values.</figcaption></figure>

leads to an enhancement of the electron density by an order or magnitude. We can perform a similar analysis in a high-density scenario ($n_H = 10^7$ cm$^{-3}$, normalized density of 1). Given the following two families of input vectors:

$$ \theta^c = [T, 1, \zeta_{H_2}, 0, \chi], $$

and

$$ \theta^d = [T, 1, \zeta_{H_2}, 1, \chi], $$

with temperature $T$, cosmic-ray ionization rate $\zeta_{H_2}$, and FUV field strength $\chi$ set as free parameters and averaged over 30 000 random samples, their respective electron densities are:

$$ \langle \log n_e(\theta^c) \rangle = -8.60, $$

$$ \langle \log n_e(\theta^d) \rangle = -8.57. $$

This results in a relative difference of

$$ \frac{\langle n_e(\theta^d) \rangle}{\langle n_e(\theta^c) \rangle} = 10^{-8.57+8.60} \sim 1.07, $$

that is, in the high-density scenario, the density of electrons is only increased by a 7 per cent when the initial sulphur abundance is enhanced from $[S] = 7.5 \times 10^{-8}$ to $[S] = 1.5 \times 10^{-5}$.

It is now worth discussing whether the link between sulphur abundance, density, and the ionization degree has relevant implications in the star formation process. As mentioned above, non-ideal MHD simulations are the state-of-the-art tool to model core collapse and star formation. They include diffusive phenomena like ambipolar diffusion, Ohmic dissipation, and Hall effect that regulate the dissipation of the magnetic flux (Masson et al. 2012). These resistivities depend on the ionization degree of gas and dust (see, e.g. Marchand et al. 2016). At densities $n_H < 10^7$ cm$^{-3}$, electrons are the dominant charge carriers while at higher densities metallic ions and grains overcome electrons becoming the main charged species (Marchand et al. 2016). This density is precisely the largest considered in our parameter space and, accordingly, the possible effects of the initial sulphur abundance on the number of electrons and the dynamics of star formation should start early on in the isothermal core collapse. In this stage, in the low density scenario $n_H \sim 10^4$ cm$^{-3}$, the initial sulphur abundance can be more effective increasing
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Magnetic resistivities and their relevance at different stages of the star formation process as a function of density. Light blue: negative Hall resistivity, dark blue: positive Hall resistivity, red: ambipolar diffusion resistivity, and green: Ohm resistivity. Taken from Marchand et al. (2016).</figcaption></figure>

or decreasing the number of electrons and, consequently, modifying the magnetic resistivities. This is particularly relevant for the role of ambipolar diffusion in the star formation process (see Fig. 8): a high initial sulphur abundance in a low-density medium leads to a higher ionization of the gas that becomes better collisionally coupled with the neutral species, reducing the ambipolar diffusion resistivity. Likewise, low initial sulphur abundances could lead to a lower electron abundance that enhances ambipolar resistivity. As density increases, the initial sulphur abundance loses efficiency at setting the amount of electrons, and eventually, metallic ions and grains become the main charge carriers. We conclude that, at the early times during the isothermal collapse that precedes the First Hydrostatic Core phase, the initial sulphur abundance can modify the ambipolar resistivity effectively in low density environments. This ability is greatly reduced as collapse proceeds and the amount of electrons decreases. While the ambipolar diffusion does not seem to have an impact on the mass and radius of the first Larson core when compared to the ideal MHD case (Masson et al. 2016), it defines the ambipolar diffusion time-scale, the time-scale in which filaments fragment and cores form (Ciolek & Mouschovias 1993; Tassis & Mouschovias 2004; Kudoh & Basu 2014; Burge et al. 2016). In the simulations carried out by Burge et al. (2016), the collapse rate is inversely proportional to the ionization degree. According to their results, a difference of an order of magnitude in the ionization degree, which is the difference we predict between the high and low initial sulphur abundance scenarios for the low density case (see equation (14)), leads to different collapse dynamics: a magnetically regulated collapse, where gas is strongly coupled to the magnetic field and the gas collapses at velocities much lower than the sound speed, and a gravity-dominated one where gas is weakly coupled to the magnetic field and reaches higher velocities. Consequently, the initial sulphur abundance and its impact on the ionization state of the gas in low density environments could potentially play a role in the dynamics of core collapse.

## 6 CONCLUSIONS

Astrochemistry is a fundamental tool for the interpretation of observations in the interstellar medium. The solution to the system of nonlinear coupled ordinary differential equations describing the time-dependent chemistry in large chemical networks is feasible

nowadays but requires lots of computing power. Having access to a fast emulator will open the possibility of improving our understanding of the physical conditions in the interstellar medium.

We have developed such an emulator making use of NFs conditioned on the physical conditions and trained with computations carried out with the Nautilus code. We have demonstrated, through validation studies, that the emulator does an excellent job in approximating the time evolution of 192 species. The logarithm of the abundance of almost all species can be obtained with uncertainties well below 0.2 dex for all times below $10^7$ yr. The resulting model is able to provide the time evolution for all species ~$10^4$ times faster than Nautilus. In addition, the computing time is not sensitive to the specific physical conditions, carrying out exactly the same number of operations for any combination of the input parameters inside the ranges defined in Table 1. This allows us to perform thousands of simulations in a reasonable computing time to fit spectroscopic observations or explore the influence of different input parameters in chemical predictions. In this paper, we have used Shapley values to investigate the influence of the sulphur elemental abundance on the gas ionization degree, key parameter in the dynamical evolution of interstellar clouds. We conclude that the initial sulphur abundance can modify the ambipolar diffusion resistivity in low density environments ($n_H \sim 10^4$ cm$^{-3}$ to a few $10^5$ cm$^{-3}$), hence determining the time-scale for fragmentation and collapse. This ability is reduced as collapse proceeds and density achieves values close to $n_H \sim 10^7$ cm$^{-3}$. In this high-density regime the influence of sulphur abundance on the dynamics is negligible.

Additionally, the model is fully differentiable. This opens the possibility of including the chemistry in more elaborate forward models, to better understand the observations. One example is the Bayesian inference of physical properties in the interstellar medium via the interpretation of observed spectral lines. The speed of our model make it feasible to systematically determine the physical and chemical conditions over large regions of the sky to characterize star-forming regions of ours and external galaxies.

Emulators, obviously, trade speed for flexibility. Since these models require a training set, the chemical network, the number of species, the range of physical parameters and all remaining hyperparameters have to be fixed in advance. Updating the chemical network (either changing the reaction rates or adding/removing reactions), including more species or exploring different regions of the space of parameters requires a retraining of the neural model. Anyway, the training time is negligible when compared with the amount of time saved in any downstream task carried out with the emulator.

## ACKNOWLEDGEMENTS

AAR thanks Carlos Díaz Baso for helpful discussions on the training of NFs. AAR acknowledges support from the Agencia Estatal de Investigación del Ministerio de Ciencia, Innovación y Universidades (MCIU/AEI) and the European Regional Development Fund (ERDF) through project PID2022-136563NB-I00. DNA acknowledges funding support from Fundación Ramón Areces through their international post-doc grant program. AF and PRM thank Spanish MICIN for funding through project PID2019-106235GB-I00 and grant PID2022-137980NB-I00 by the Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/ 10.13039/501100011033 and by ‘ERDF A way of making Europe’. AF thanks the European Research Council (ERC) for funding under the Advanced Grant project SUL4LIFE, grant agreement no. 101096293. This research made use of the IAC HTCondor facility (http://research.cs.wisc.edu/htcondor/), partly fi-
nanced by the Ministry of Economy and Competitiveness with FEDER funds, code IACA13-3E-2493. This research has made use of NASA's Astrophysics Data System Bibliographic Services. We acknowledge the community effort devoted to the development of the following open-source packages that were used in this work: numpy ( numpy.org, Harris et al. 2020), matplotlib ( matplotlib.org, Hunter 2007), PyTorch ( pytorch.org, Paszke et al. 2019), zarr ( github.com/zarr-developers/zarr-python) and SHAP ( Lundberg & Lee 2017, github.com/shap/shap).

DATA AVAILABILITY

Access to the data and codes used in this work is described in the repository https://github.com/aasensio/neuralchemistry.

REFERENCES

Allen A., Li Z.-Y., Shu F. H., 2003, ApJ, 599, 363

Asensio Ramos A., 2023, Sol. Phys., 298, 135

Bergin E. A., Plume R., Williams J. P., Myers P. C., 1999, ApJ, 512, 724

Branca L., Pallottini A., 2023, MNRAS, 518, 5718

Branca L., Pallottini A., 2024, A&A, 684, 203

Bron E. et al., 2021, A&A, 645, A28

Bulut N. et al., 2021, A&A, 646, A5

Burge C. A., Van Loo S., Falle S. A. E. G., Hartquist T. W., 2016, A&A, 596, A28

Cabedo V., Maury A., Girart J. M., Padovani M., Hennebelle P., Houde M., Zhang Q., 2023, A&A, 669, A90

Caselli P., Walmsley C. M., Terzieva R., Herbst E., 1998, ApJ, 499, 234

Cernicharo J. et al., 2023, in European Conference on Laboratory Astrophysics ECLA2020. The Interplay of Dust. The Interplay of Dust, Ice, and Gas in Space, New York, p. 35

Ciolek G. E., Mouschovias T. C., 1993, ApJ, 418, 774

Clément A. et al., 2023, A&A, 675, A165

Commerçon B., Hennebelle P., Audit E., Chabrier G., Teyssier R., 2010, A&A, 510, L3

Commerçon B., Launhardt R., Dullemond C., Henning T., 2012, A&A, 545, A98

Coutens A., Commerçon B., Wakelam V., 2020, A&A, 643, A108

Crockett N. R., Bergin E. A., Neill J. L., Favre C., Blake G. A., Herbst E., Anderson D. E., Hassel G. E., 2015, ApJ, 806, 239

Cybenko G., 1988, Technical Report, Approximation by Superpositions of a Sigmoidal Function. inst-tufts-cs

de Mijolla D., Viti S., Holdship J., Manolopoulou I., Yates J., 2019, A&A, 630, A117

Drozdovskaya M. N., Walsh C., van Dishoeck E. F., Furuya K., Marboeuf U., Thiabaud A., Harsono D., Visser R., 2016, MNRAS, 462, 977

Drozdovskaya M. N., Walsh C., Visser R., Harsono D., van Dishoeck E. F., 2014, MNRAS, 445, 913

Favre C. et al., 2017, A&A, 608, A82

Friesen R. K. et al., 2017, ApJ, 843, 63

Fuente A. et al., 2014, A&A, 568, A65

Fuente A. et al., 2019, A&A, 624, A105

Galli D., Lizano S., Shu F. H., Allen A., 2006, ApJ, 647, 374

Gerin M., Liszt H., Neufeld D., Godard B., Sonnentrucker P., Pety J., Roueff E., 2019, A&A, 622, A26

Gerin M., Neufeld D. A., Goicoechea J. R., 2016, ARA&A, 54, 181

Goicoechea J. R., Pety J., Gerin M., Teyssier D., Roueff E., Hily-Blant P., Baek S., 2006, A&A, 456, 565

Gong M., Ivlev A. V., Akimkin V., Caselli P., 2021, ApJ, 917, 82

Gong M., Ivlev A. V., Zhao B., Caselli P., 2020, ApJ, 891, 172

Grassi T., Nauman F., Ramsey J. P., Bovino S., Picogna G., Ercolano B., 2022, A&A, 668, A139

Gratier P. et al., 2021, A&A, 645, A27

Guélin M., Cernicharo J., 2022, Front. Astron. Space Sci., 9, 787567

Hacar A., Bosman A. D., van Dishoeck E. F., 2020, A&A, 635, A4

Harris C. R. et al., 2020, Nature, 585, 357

Hasegawa T. I., Herbst E., 1993, MNRAS, 261, 83

Hasegawa T. I., Herbst E., Leung C. M., 1992, ApJS, 82, 167

Hennebelle P., Commerçon B., Lee Y.-N., Charnoz S., 2020, A&A, 635, A67

Hennebelle P., Fromang S., 2008, A&A, 477, 9

Heyl J., Butterworth J., Viti S., 2023, MNRAS, 526, 404

Hincelin U., Commerçon B., Wakelam V., Hersant F., Guilloteau S., Herbst E., 2016, ApJ, 822, 12

Hocuk S., Szűcs L., Caselli P., Cazaux S., Spaans M., Esplugues G. B., 2017, A&A, 604, A58

Holdship J., Viti S., Haworth T. J., Ilee J. D., 2021, A&A, 653, A76

Holdship J., Viti S., Jiménez-Serra I., Makrymallis A., Priestley F., 2017, AJ, 154, 38

Hsu C.-J., Tan J. C., Goodson M. D., Caselli P., Körtgen B., Cheng Y., 2021, MNRAS, 502, 1104

Hsu C.-J., Tan J. C., Holdship J., Duo, Xu, Viti S., Wu B., Gaches B., 2023, preprint (arXiv:2308.11803)

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Jarolim R. et al., 2022, in AGU Fall Meeting Abstracts, Chicago. p. SH45D-2368

Jarolim R. et al., 2024, ApJL, 961, L31

Jørgensen J. K., Belloche A., Garrod R. T., 2020, ARA&A, 58, 727

Kingma D. P., Ba J., 2014, preprint (arXiv:1412.6980)

Kudoh T., Basu S., 2014, ApJ, 794, 127

Laas J. C., Caselli P., 2019, A&A, 624, A108

Le Petit F., Nehmé C., Le Bourlot J., Roueff E., 2006, ApJS, 164, 506

Lebreuilly U. et al., 2024, A&A, 682, A30

Lin C.-H., Ma W.-C., Torralba A., Lucey S., 2021, inarxiv:2104.06405

Loshchilov I., Hutter F., 2017, SGDR: Stochastic Gradient Descent with Warm Restarts, preprint (arXiv:1608.03983)

Lu L., Jin P., Pang G., Zhang Z., Karniadakis G. E., 2021, Nat. Mach. Intell., 3, 218

Lundberg S. M., Lee S.-I., 2017, in von Luxburg U., Guyon I., Bengio S., Wallach H., Fergus R., eds, Proceedings of the 31st International Conference on Neural Information Processing Systems. NIPS'17. Curran Associates Inc., Red Hook, NY, USA, p. 4768

Maes S., Ceuster F. D., de Sande M. V., Decin L., 2024, MACE: A Machine learning Approach to Chemistry Emulation. preprint (arXiv:2405.03274)

Marchand P., Masson J., Chabrier G., Hennebelle P., Commerçon B., Vaytet N., 2016, A&A, 592, A18

Masson J., Chabrier G., Hennebelle P., Vaytet N., Commerçon B., 2016, A&A, 587, A32

Masson J., Teyssier R., Mulet-Marquis C., Hennebelle P., Chabrier G., 2012, ApJS, 201, 24

McGuire B. A. et al., 2020, ApJ, 900, L10

Mildenhall B., Srinivasan P. P., Tancik M., Barron J. T., Ramamoorthi R., Ng R., 2020, in European Conference on Computer Vision, Glasgow. p. 405

Minissale M., Dulieu F., Cazaux S., Hocuk S., 2016, A&A, 585, A24

Müller H. S. P., Thorwirth S., Roth D. A., Winnewisser G., 2001, A&A, 370, L49

Nair V., Hinton G. E., 2010, in Proceedings of the 27th International Conference on Machine Learning (ICML-10), 2010 June 21–24, Haifa, Israel, p. 807

Navarro-Almaida D. et al., 2021, A&A, 653, A15

Navarro-Almaida D. et al., 2024, A&A, 685, A112

Neufeld D. A., Wolfire M. G., 2017, ApJ, 845, 163

Öberg K. I., Facchini S., Anderson D. E., 2023, ARA&A, 61, 287

Paszke A. et al., 2019, in Wallach H., Larochelle H., Beygelzimer A., d'é Buc F., Fox E., Garnett R. eds, Advances in Neural Information Processing Systems 32. Curran Associates, Inc., p. 8024, http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf
Perez E., Strub F., de Vries H., Dumoulin V., Courville A., 2017, preprint (arXiv:1709.07871)

Pineda J. E. et al., 2024, preprint (arXiv:2402.16202)

Priestley F. D., Viti S., Williams D. A., 2018, AJ, 156, 51

Rahaman N., Baratin A., Arpit D., Draxler F., Lin M., Hamprecht F., Bengio Y., Courville A., 2019, in Chaudhuri K., Salakhutdinov R.eds, Proceedings of the 36th International Conference on Machine Learning, vol. 97. PMLR, p. 5301, https://proceedings.mlr.press/v97/rahaman19a.html

Rivière-Marichalar P. et al., 2020, A&A, 642, A32

Ruaud M., Wakelam V., Hersant F., 2016, MNRAS, 459, 3756

Shapley L. S., 1953, 17. A Value for n-Person Games. Princeton University Press, Princeton, p. 307, https://doi.org/10.1515/9781400881970-018

Sitzmann V., Martel J., Bergman A., Lindell D., Wetzstein G., 2020, Adv. Neural Inf. Process. Syst., 33, 7462

Sulzer I., Buck T., 2023, preprint (arXiv:2312.06015)

Tancik M. et al., 2020, Adv. Neural Inf. Process. Syst., 33, 7537

Taquet V., Ceccarelli C., Kahane C., 2012, A&A, 538, A42

Tassis K., Mouschovias T. C., 2004, ApJ, 616, 283

Wakelam V., Dartois E., Chabot M., Spezzano S., Navarro-Almaida D., Loison J. C., Fuente A., 2021, A&A, 652, A63

Wakelam V., Loison J. C., Mereau R., Ruaud M., 2017, Mol. Astrophys., 6, 22

Wang S., Wang H., Perdikaris P., 2021, Comput. Methods Appl. Mech. Eng., 384, 113938

Zhao B., Caselli P., Li Z.-Y., 2018b, MNRAS, 478, 2723

Zhao B., Caselli P., Li Z.-Y., Krasnopolsky R., 2018a, MNRAS, 473, 4868

## APPENDIX A: QUANTITATIVE RESULTS FOR ALL MOLECULES

Figs A1–A3 display the distribution of differences in logarithm for all species in the validation set, which complement those shown in Fig. 5.

This paper has been typeset from a TEX/LATEX file prepared by the author.

<figure><img src="image_10.png" /><figcaption>**Figure A1.** Same as Fig. 5 for additional species.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure A2. Continuation of Fig. A1.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure A3.** Continuation of Fig. A2.</figcaption></figure>