Using a neural network approach to accelerate disequilibrium chemistry
  calculations in exoplanet atmospheres
In this era of exoplanet characterisation with JWST, the need for a fast
implementation of classical forward models to understand the chemical and
physical processes in exoplanet atmospheres is more important than ever.
Notably, the time-dependent ordinary differential equations to be solved by
chemical kinetics codes are very time-consuming to compute. In this study, we
focus on the implementation of neural networks to replace mathematical
frameworks in one-dimensional chemical kinetics codes. Using the gravity
profile, temperature-pressure profiles, initial mixing ratios, and stellar flux
of a sample of hot-Jupiters atmospheres as free parameters, the neural network
is built to predict the mixing ratio outputs in steady state. The architecture
of the network is composed of individual autoencoders for each input variable
to reduce the input dimensionality, which is then used as the input training
data for an LSTM-like neural network. Results show that the autoencoders for
the mixing ratios, stellar spectra, and pressure profiles are exceedingly
successful in encoding and decoding the data. Our results show that in 90% of
the cases, the fully trained model is able to predict the evolved mixing ratios
of the species in the hot-Jupiter atmosphere simulations. The fully trained
model is ~1000 times faster than the simulations done with the forward,
chemical kinetics model while making accurate predictions.
Julius L. A. M. Hendrix, Amy J. Louca, Yamila Miguel
Using a neural network approach to accelerate disequilibrium chemistry calculations in exoplanet atmospheres

###### Abstract

In this era of exoplanet characterisation with JWST, the need for a fast implementation of classical forward models to understand the chemical and physical processes in exoplanet atmospheres is more important than ever. Notably, the time-dependent ordinary differential equations to be solved by chemical kinetics codes are very time-consuming to compute. In this study, we focus on the implementation of neural networks to replace mathematical frameworks in one-dimensional chemical kinetics codes. Using the gravity profile, temperature-pressure profiles, initial mixing ratios and stellar flux of a sample of hot-Jupiter atmospheres as free parameters, the neural network is built to predict the mixing ratio outputs in steady state. The architecture of the network is composed of individual autoencoders for each input variable to reduce the input dimensionality, which is then used as the input training data for an LSTM-like neural network. Results show that the autoencoders for the mixing ratios, stellar spectra, and pressure profiles are exceedingly successful in encoding and decoding the data. Our results show that in 90% of the cases, the fully trained model is able to predict the evolved mixing ratios of the species in the hot-Jupiter atmosphere simulations. The fully trained model is \(\sim 10^{3}\) times faster than the simulations done with the forward, chemical kinetics model while making accurate predictions.

keywords: planets and satellites: gaseous planets - planets and satellites: atmospheres - exoplanets

## 1 Introduction

There are two methods commonly used for calculating the abundance of different species in an atmosphere: thermochemical equilibrium and chemical kinetics (Bahn & Zukoski, 1960; Zeleznik & Gordon, 1968). Thermochemical equilibrium calculations treat each species independently and do not require an extensive list of reactions between different species. Consequently, this method is fast for estimating the abundance of different species in an exoplanet atmosphere and has been widely used in the community (e.g., Stock et al., 2018; Woitke et al., 2018). However, the atmospheres of exoplanets are dynamic environments. Both physical and chemical processes can alter the compositions and thermal structures of the atmosphere. In particular, atmospheric processes like photochemistry, mixing and condensation of different species can affect atmospheric abundances, deviating the concentrations observed from what would be found by chemical equilibrium calculations (Cooper & Showman, 2006; Swain et al., 2008; Moses et al., 2011; Kawashima & Min, 2021; Roudier et al., 2021; Baxter et al., 2021). For example, the recent detection of SO\({}_{2}\)(Feinstein et al., 2022; Ahrer et al., 2022; Alderson et al., 2022; Rusatmarkov et al., 2022) and the determination of this species as direct evidence of photo-chemical processes shaping the atmosphere of WASP 39b (Tsai et al., 2022), suggest that certain exoplanet atmospheres are in disequilibrium and we need chemical disequilibrium models using chemical kinetics to correctly interpret the observations.

Chemical kinetics codes consider the effects that lead to a non-equilibrium state in the atmosphere. These codes incorporate a wide range of atmospheric processes such as the radiation from the host start that can drive the dissociation of molecules -or photochemistry-, the mixing of species at different pressures due to the planet's winds, or the diffusion of species, and calculate the one-dimensional abundances of species in exoplanetary atmospheres (e.g. Moses et al., 2011; Venot et al., 2012; Miguel & Kaltenegger, 2014; Tsai et al., 2017; Hobbs et al., 2019). However, To calculate the abundance of different species using chemical kinetics, a system of coupled differential equations involving all the species must be solved, and prior knowledge of reaction rates and a reaction list is necessary to estimate the production and loss of each species. Therefore, as more species and reactions are incorporated into the chemical networks, the complexity of these simulations increases, and so does the computational cost these simulations require. The result is that chemical kinetics codes have long computational times, and can not be used by more detailed calculations (e.g. circulation models) or as a fast way of interpreting observations (by retrieval codes), which are usually subject to simplifications.

For the past few decades, the use of machine learning techniques, specifically neural networks (NN), has become more prevalent in research fields outside of computer science. Within astronomy, neural networks have been used for applications like image processing (Dattilo et al., 2019), adaptive optics (Landman et al., 2021), exoplanet
detection (Shallue and Vanderburg, 2018), exoplanetary atmospheric retrieval (Cobb et al., 2019) and chemical modelling (Holdship et al., 2021), and more traditional machine learning techniques have been used for applications like exoplanetary atmospheric retrieval (Nixon and Madhusudhan, 2020) and chemistry modelling of protoplanetary disks (Smirnov-Pinchukov et al., 2022). Trained neural networks are fast to use, so a neural network trained to accurately reproduce the outcomes of chemical kinetics codes could greatly reduce computational time. Such a neural network could simulate a large amount of atmospheric conditions in a short period of time, which is for example useful for atmospheric retrievals from observational constraints. It could also be incorporated into a multi-dimensional atmospheric simulation that connects a multitude of individual one-dimensional simulations by the implementation of atmospheric mixing and other global processes.

In this study, we investigate the feasibility of machine learning techniques for speeding up a one-dimensional chemical kinetics code. To this end, we perform calculations on a fiducial giant planet as an example to show how this technique can be used to bring the best of these two worlds: the detailed information of chemical kinetics calculations and the speed of Neural Networks techniques. In the next section, we explain in more detail how we obtain the dataset and the specifics of the architectures used. The results of our networks are presented in the following section (section 3), which are discussed afterwards in section 4. Finally, we summarise and conclude our findings in section 5.

## 2 Methods

### Chemical Kinetics

Chemical kinetics is the most realistic way of calculating abundances and is necessary, particularly at low temperatures (T < 2000 K) and pressures (P < 10 - 100 bars), where timescales of processes such as mixing in the atmosphere are shorter than chemical equilibrium and dominate the chemistry and abundances in the atmosphere.

We make use of the one-dimensional chemical kinetics code VULCAN (Tsai et al., 2017, 2021), to create a large dataset on the atmospheres of gaseous exoplanets. The code is validated for hot-Jupiter atmospheres from 500 K to 2500 K. VULCAN calculates a set of mass differential equations:

\[\frac{\partial n_{i}}{\partial t}=\mathcal{P}_{i}-\mathcal{L}_{i}-\frac{ \partial\Phi_{i}}{\partial z}, \tag{1}\]

where \(n_{i}\) is the number density of the species \(i\), \(t\) is the time, \(\mathcal{P}_{i}\) and \(\mathcal{L}_{i}\) are the production and loss rates of the \(i\)-th species, and \(\Phi_{i}\) its transport flux that includes the effects of dynamics caused by convection and turbulence in the atmosphere. For a more complete derivation of this equation from the general diffusion equation, we refer the reader to Hu et al. (2012). VULCAN starts from initial atmospheric abundances calculated using the chemical equilibrium chemistry code _FastChem_(Stock et al., 2018), although we note that the final disequilibrium abundances are not affected by the choice of initial values adopted Tsai et al. (2017), and further evolves these abundances by solving a set of Eulerian continuity equations that includes various physical processes (e.g. vertical mixing and photochemistry). To solve these partial differential equations, VULCAN numerically transforms them into a set of stiff ordinary differential equations (ODEs). These ODEs are solved using the _Rosenbrock_ method, which is described in detail in the appendix of Tsai et al. (2017). In this study, we make use of machine learning techniques to solve these equations and hence speed up the process.

### Building the dataset

#### Parameter Space

To construct the dataset, we vary the following parameters:

1. **Planet mass, M** [M\({}_{\rm J}\)]: within the range [0.5, 20] M\({}_{\rm J}\).
2. **Orbit radius, r** [AU]: within the range [0.01, 0.5] AU.
3. **Stellar radius, R\({}_{\star}\)** [R\({}_{\odot}\)]: within the range [1, 1.5] R\({}_{\odot}\).

Other parameters such as surface gravity, irradiation temperature, and stellar effective temperature are derived from these free parameters.

1. **Planet radius** [R\({}_{\rm Jup}\)]: This is derived from the planet mass using the relation from Chen and Kipping (2017), shown in Equation 2, where \(R\) is the planet radius and \(M\) is the planet mass: \[\frac{R}{R_{\oplus}}=17.78\left(\frac{M}{M_{\oplus}}\right)^{-0.044}.\] (2)

We note that our aim is to present the results for a simple general case, and the mass-radius relation we use is suitable for this purpose. However, we must emphasize that the relation between mass and radius for giant exoplanets is not unique and depends on various factors, such as the mass of metals, core mass, irradiation received by the planet, and their effect on the inflation of radius. All of these factors can impact the evolution path of giant planets and their final radius, leading to a dispersion in the mass and radius relation.

1. **Temperature-pressure profile**: As our aim is to demonstrate the use of neural networks for calculating non-equilibrium chemical abundances in a general case, we have utilized an analytical, non-inverted temperature-pressure profile from Heng et al. (2014). While these analytical profiles are simplistic, they are widely used in the literature to explore general cases and are suitable for our purposes. However, for calculating the chemistry of a real planet, more detailed calculations that take into account the opacities of different species and their abundances in the atmosphere should be included. The assumptions for this calculation are \(T_{int}=120\) K, \(\kappa_{L}=0.1\), \(\kappa_{S}=0.02\), \(\beta_{S}=1\) and \(\beta_{L}=1\), based on the default values included in the VULCAN code Tsai et al. (2017). The pressure profile is constructed within the range [\(10^{-2}\), \(10^{9}\)] dyne cm\({}^{-2}\). This calculation is an important step, as it determines whether the set of parameters is valid for the dataset. If any part of the temperature profile falls outside of the range [500, 2500] K, the temperature range for which VULCAN is validated, the example is rejected from the dataset.
2. **Stellar flux**: the stellar spectra used for the dataset have two sources: the Measurements of the Ultraviolet Spectra Characteristics of Low-mass Exoplanetary Systems (MUSCLES) collaboration (France et al., 2016; Youngblood et al., 2016; Loyd et al., 2016) and the PHOENIX Stellar and Planetary Atmosphere Code (Baron et al., 2010). The MUSCLES database contains observations from M- and K-dwarf exoplanet host stars in the optical, UV, and X-ray regime, and is used for stars with an effective temperature lower than 6000 K. For effective temperatures of 6000 K and above, stellar spectra are generated by the PHOENIX model. Flux values below \(10^{-14}\) erg nm\({}^{-1}\) cm\({}^{-2}\) s\({}^{-1}\) are cut-off.

The remaining parameters of the VULCAN configuration files are kept constant throughout the dataset. Eddy- and molecular diffusion are both taken into account as fixed parameters. For the eddy diffusion constant, we make use of a constant of \(K_{zz}=10^{10}\) cm\({}^{2}\)/s. The molecular diffusion constant is taken for a hydrogen-dominated gas as described in Banks and Kockarts (1973). As standard chemical
network, we make use of VULCAN's reduced default N-C-H-O network that includes photochemistry. We assume experimental abundances for the hot-Jupiters, and we make use of 150 pressure levels for the height layers. The output of VULCAN is saved every 10 steps. In total, 13291 valid configurations are generated within the parameter space.

#### 2.2.2 Formatting

In order to limit the computation times when training the network, the input-output pairs do not contain all of the information supplied by VULCAN.

For the inputs, a selection of six properties is made. These properties are extracted from VULCAN before time integration starts, so they can be interpreted as the initial conditions of the simulation. The six properties are:

1. **Initial mixing ratios**: the initial mixing ratios of the species in the simulation. These mixing ratios are calculated by VULCAN using FastChem (Stock et al., 2018). The shape of the array containing the mixing ratios is (69, 150), as the mixing ratios are defined for 69 species for 150 height layers each.
2. **Temperature profile**: the temperature profile as calculated by the analytical expression from Heng et al. (2014). The temperature is defined for every height layer, so it has a shape of (150).
3. **Pressure profile**: the pressure profile that is calculated as part of the temperature profile calculation. It is of the same shape, (150).
4. **Gravitational profile**: the gravitational acceleration per height layer. It has the shape (150,).
5. **Stellar flux component**: one of the two components that make up the stellar spectrum contains the flux values. This is generated from either the MUSCLES database and/or the PHOENIX model and interpolated to a shape of (2500,).
6. **Stellar wavelength component**: the second component of the stellar spectrum contains the wavelengths corresponding to the flux values. It has the same shape, (2500,).

For the outputs, we make use of the time-dependent mixing ratios.

Because not every simulation takes the same amount of time to converge to a solution, the number of saved abundances differs per VULCAN simulations. To include the information contained in the evolution of the abundances through time, 10 sets of abundances, including the steady-state abundances, are saved in each output. This set of abundances is evenly spaced through time, so the simulation time between abundances will vary for different VULCAN simulation runs. Before the abundances are saved, they are converted to mixing ratios. The shape of the outputs is (10, 69, 150).

#### 2.2.3 Data Standardisation

The inputs and outputs of the various components differ by several orders in magnitude. To ensure that the neural network trained on the data set is not biased towards higher-valued parameters, the data has to be standardised. First, the distributions of the properties are standardised according to Equation 3:

\[p_{s}=\frac{\log_{10}(p)-\mu}{\sigma}, \tag{3}\]

with

\[\mu=\frac{1}{n}\sum_{i=0}^{n}\log_{10}(p_{i}), \tag{4}\]

and

\[\sigma=\sqrt{\frac{1}{n}\sum_{i=0}^{n}\left(\log_{10}(p_{i})-\mu\right)^{2}}, \tag{5}\]

where \(p\) is the property to be scaled, \(n\) is the size of the dataset and \(p_{s}\) is the standardised property. After standardisation, the properties are normalised in the range [0, 1]:

\[p_{s,n}=\frac{p_{s}-\min(p_{s})}{\max(p_{s})-\min(p_{s})}, \tag{6}\]

where \(p_{s,n}\) is the final normalised property.

Once the input properties are normalised, the output mixing ratios are normalised with the same scaling parameters as were used for the input mixing ratios. When the trained neural network is presented with an input for which to predict the mixing ratios, it only has information about the scaling parameters of the inputs. To be able to unnormalise the outputs, they need to be scaled with the same scaling parameters.

### Model Architecture

#### 2.3.1 Autoencoder Structure

The input of each configuration within the dataset consists of roughly 15800 values. To speed up the training process and complexity of the neural network we make use of an _autoencoder_ (AE) for reducing the dimensionality of the examples in the dataset. In previous studies this approach has been shown an effective way to reduce dimensionality within chemical kinetics (e.g. Grassi et al., 2022). An autoencoder consists of two collaborating neural networks: an _encoder_ and a _decoder_. The task of the encoder is to reduce the dimensionality of the input data by extracting characterising features from the example and encoding them in a lower dimensionality representation called the _latent representation_. The task of the decoder is to take the latent representation and use it to reconstruct the original input data with as little loss of information as possible. The encoder and decoder are trained simultaneously, and no restraints are placed on the way the autoencoder uses its _latent space_, apart from the size of the latent representations.

As is discussed in Section 2.2.2, the inputs of the model consist of six properties. Because these properties do not share the same shape, we cannot encode and decode them using a single autoencoder. Instead, we construct six unique autoencoders, one for each property of the model inputs. Figure 1 shows an overview of the process of encoding the initial conditions. The decoding process is not shown but is symmetrical to the encoding process. Each encoder conceals a specific property into a corresponding latent representation. To get the latent presentation of the entire input example, \(l_{i}\), the property latent representations \(\{l_{MR},l_{F},l_{W},l_{T},l_{P},l_{G}\}\) are concatenated. When decoding the latent representation of the input, the latent vector \(l_{i}\) is split back into the different property latent representations, and given to that property's decoder. Every encoder-decoder pair is trained separately.

The hyperparameters of each autoencoder are optimized by trial and error. A summary of each set of hyperparameters is shown in table 1.
### Mixing Ratios Autoencoder

The mixing ratios are the largest contributor to the size of the model inputs. Compressing each of the species' mixing ratios efficiently reduces the size of the input latent representation, \(l_{i}\), by a substantial amount. Because of the limited size of the training dataset, a compromise has to be made to successfully train this autoencoder. Rather than concurrently encoding all 69 species for each example, each species is encoded individually. As a result, the training dataset expands by a factor of 69, while disregarding any potential correlations in species abundances during the encoding procedure.

Figure 2 shows the application of such an autoencoder: for a given input, each of the 69 species' mixing ratios is encoded into corresponding latent vectors \(\{l_{1},l_{2},l_{3},...,l_{69}\}\). The concatenation of these 69 latent vectors then makes up the latent representation of the mixing ratios \(l_{MR}\).

All encoders and decoders are multilayer perceptron (MLP) neural networks. For the mixing ratio autoencoder (MRAE), the encoder and decoder both consist of 7 fully connected layers, followed by hyperbolic tangent activation functions. The encoder input layer has a size of 150 and the output layer has a size of 30. The hidden layers have a size of 256. Adversely, the decoder has an input layer size of 30 and an output layer size of 150. The compression factor of the MRAE is therefore \(150/30=5\).

To train the MRAE, the dataset is split into a train dataset (70%), a validation dataset (20%), and a test dataset (10%). To increase the size of the dataset, the MRAE is trained on a shuffled set1 of the mixing ratios of both the inputs and the output of the chemical kinetics simulations. The performance of the autoencoder is measured using the loss function in equation 7:

Footnote 1: Using a random sampler function within the PyTorch package.

\[\mathcal{L}=\frac{1}{N}\sum_{i=1}^{N}\left(\frac{p_{i}-a_{i}}{a_{i}}\right)^{ 2}. \tag{7}\]

Figure 1: The different properties and their corresponding encoders that together encode the input data. The decoding process is symmetrical to this encoding process, where every property has a corresponding property decoder.

Figure 2: A more detailed sketch of the architecture of the mixing ratio autoencoder. MR \(i\) denotes the mixing ratio of a certain species, \(i\), for all height layers, and \(l_{i}\) denotes the encoded mixing ratios for species \(i\).
where \(\mathcal{L}\) is the loss, \(N\) is the number of elements in the actual/predicted vector, and \(p_{i}\) and \(a_{i}\) are \(i\)-th elements of the predicted and actual vectors, respectively. The MRAE is optimised using the Adam optimiser (Kingma and Ba, 2014), with a learning rate of \(10^{-5}\). A batch size of 32 is used, and the model is trained for 200 epochs. These hyperparameters can also be found in table 1.

### Atmospheric profile Autoencoders

The temperature, pressure, and gravity profiles all have the same shape of (150,), so their autoencoders can use the same architecture. Moreover, the atmospheric profiles share their shape with the mixing ratios of individual species (i.e. the height layers in the atmosphere). Therefore, a very similar neural network structure to that of the MRAE is used for the atmospheric profile autoencoders. The encoder input layer shape, the decoder output layer shape, and the hidden layer shapes are taken directly from the MRAE for all atmospheric profile autoencoders. An important parameter to tune for each atmospheric profile autoencoder separately is the size of the latent representations.

The pressure profile is set logarithmically for all examples in the dataset. By taking the logarithm of the pressures the spacing becomes linear, \(\log(P_{i})-\log(P_{i+1})\). Theoretically, we then only need two values to fully describe the pressure profile: the pressure at the first and last height layers. To encode these values, no autoencoder is needed. One could take this one step further, and provide input parameters like mass and radius, from which the pressure and gravity profile are dependent, directly to the core model as inputs. While this more specialised approach is suitable for these input parameters, it is not generalisable to other input parameters. To keep the model architecture more general and adaptable to different input parameters, an autoencoder is used nonetheless. The size of the pressure profile autoencoder (PAE) latent representations is set to 2. This corresponds to a compression factor of \(150/2=75\).

The temperature and gravity profiles are not linear. For both the temperature autoencoder (TAE) and gravity autoencoder (GAE), a latent representation size of 30 is used. This corresponds to a compression factor of \(150/30=5\), the same as for the MRAE.

All profile autoencoders are evaluated using the loss function previously defined in equation 7 and are optimised using the Adam optimiser. The TAE and GAE use a learning rate of \(10^{-5}\), and the PAE uses a learning rate of \(10^{-6}\). All profile autoencoders are trained with a batch size of 4, for 100 epochs (see also table 1).

### Stellar Spectrum Autoencoders

After the mixing ratios, the stellar spectrum components contribute predominantly to the input data size. The stellar spectrum is comprised of a flux and a wavelength component. These components share the same shape, so one NN structure can be used for both autoencoders. The structure of the encoder and decoder is similar to that of the MRAE: a 7-layer, fully connected MLP, with hyperbolic tangent activation functions after each layer. The encoder input layer and decoder output layer have a size of 2500, and the hidden layers have a size of 1024.

Similarly to the PAE, the wavelength bins are spaced logarithmically. Again only two values are needed to fully describe the wavelength range. The latent representation size for the wavelength autoencoder (WAE) is, therefore, also 2. The compression factor for this network is \(2500/2=1250\). The flux autoencoder (FAE) has a latent representation size of 256, which gives it a compression factor of \(2500/256\approx 10\).

Both autoencoders are evaluated using the loss function from equation 7. They are optimised using the Adam optimiser, the WAE with a learning rate of \(10^{-7}\), and the FAE with a learning rate of \(10^{-5}\). They are both trained for 200 epochs with batches of 4 examples (see also table 1).

#### 2.3.2 Core Network

As mentioned before, the outputs are also large in dimensionality. Because the outputs contain mixing ratios for all species, for 10-time steps (Section 2.2), they can be encoded using the MRAE. Figure 3 shows how the autoencoder would encode both the inputs and the last time step of the outputs to their latent representations \(l_{i}\) and \(l_{o}\), respectively. Note that even though the autoencoder is shown twice in this figure, the same autoencoder is used to encode both the inputs and the outputs. In the middle of the figure, connecting the two latent spaces, a second neural network called the _core network_ is located.

The function of the core network is to learn a mapping between the latent representations of the inputs and the evolved outputs. The design of the core network takes advantage of some of the characteristics of VULCAN. From Section 2.1 we know that VULCAN solves ODEs for specific atmospheric configurations over a simulated period of time. To impart this sense of time in the core neural network, a _Long-Short Term Memory_ (LSTM) is used as the base of the design. The LSTM was chosen for its proven performance in numerous applications, from stellar variability (e.g. Jamal and Bloom, 2020) to Core-collapse supernovae search (e.g. Iess et al., 2023) and solar radio spectrum classification (e.g. Xu et al., 2019), as well as the ease of implementation. The LSTM has known shortcomings like the vanishing gradient problem and long training times when dealing with long sequences. However, with the short sequence length used with our model (i.e. 10 timesteps), these shortcomings are not considered problematic for this proof of concept.

The input of the core network is not sequential in nature. With some changes, we can use the LSTM in a 'one-to-many' configuration. In this configuration, the initial output of the LSTM \(h_{0}\) is given to an MLP. This MLP produces a vector with the same shape as the initial input \(x_{0}\), which can be interpreted as the 'evolved' input \(x_{1}\). This evolved input is fed back into the LSTM to produce \(h_{1}\), from which the MLP produces \(x_{2}\), and so forth. This can be repeated for an arbitrary number of steps.

Figure 3: An overview of the model architecture. The core neural network maps between the latent representations of the VULCAN inputs \(l_{i}\) and outputs \(l_{o}\).
The design of the core network is visualised in Figure 4. We interpret the latent representation of the inputs \(l_{i}\) as the initial value \(x_{0}\). The LSTM and MLP configuration produces 9 intermediary 'evolved' latent representations \(\{x_{1},...,x_{9}\}\) before arriving at the final evolved latent representation \(x_{10}\). We interpret this latent representation as the prediction of the latent representation of the evolved output \(I_{o}\).

#### Training

When the core model predicts a sequence of 10 latent representations, it is essentially traversing the latent space. We can guide the network to learn to traverse the latent space similarly to how VULCAN simulations evolve by using the sequence of outputs saved in the dataset (Section 2.2). We do this in two ways: first, we construct a loss function that not only depends on the accuracy of the prediction of the latent representation of the final output \(I_{o}\), but also on the accuracy of the intermittent latent representation predictions:

\[\mathcal{L}=\sum_{t=1}^{10}\left(\frac{1}{N}\sum_{i=1}^{N}\left(p_{t,i}-a_{t, t}\right)^{2}\right), \tag{8}\]

where \(\mathcal{L}\) is the loss, \(N\) is the number of elements in the actual/predicted vector, \(p_{t,i}\) is the \(i\)-th element of the latent representation prediction vector at time step \(t\), and \(a_{t,i}\) is the \(i\)-th element of the latent representation vector of the output at time step \(t\). With this notation, \(a_{10}=I_{o}\). By training a network with this loss function, we force the core network to evolve the latent mixing ratios similarly to how VULCAN evolves mixing ratios. It should be noted that the latent representation of the inputs \(l_{i}\) is larger than the latent representation of the outputs \(I_{o}\), as it contains more properties than just mixing ratios. The predicted latent representations \(x_{t}\) are therefore also larger than \(I_{o}\). To account for this, we only look at the elements corresponding to the encoded mixing ratios in \(l_{i}\) when comparing the predicted latent representations \(x_{t}\) and the output mixing ratios \(a_{t}/I_{o}\).

To further incentivise the core network to adhere to VULCAN's evolution patterns, we can intercept the predicted latent representations \(x_{t}\) before they get fed back into the LSTM, and replace them with the latent representation of the actual output of the corresponding time step. This way, the core network is always learning from latent representations that follow VULCAN's evolution, even if the network is predicting poorly. This is only done during the training of the network when the true VULCAN outputs are known. During validation and testing, the predicted latent representations \(x_{t}\) are not altered.

The core model LSTM has a hidden- and cell size of 4096. The MLP has only two layers: an input layer of size 4096, and an output layer of the same size as the latent representation of the inputs \(l_{i}\), followed by a hyperbolic tangent function. It is optimised with the Adam optimiser, with a learning rate of \(10^{-4}\) and a batch size of 8. It is trained for 100 epochs.

#### 2.3.3 Deployment

When the trained model is deployed on the validation- and test dataset, the first step is encoding the inputs into their latent representation \(l_{i}\) using the encoder part of the autoencoder (top left

\begin{table}
\begin{tabular}{l c c c c c}
**model** & **hidden size** & **latent size** & **optimiser** & **learning rate** & **batch size** & **epochs** \\ \hline \hline MRAE & 256 & 30 & Adam & \(10^{-5}\) & 32 & 200 \\ PAE & 256 & 2 & Adam & \(10^{-6}\) & 4 & 100 \\ TAE & 256 & 30 & Adam & \(10^{-5}\) & 4 & 100 \\ GAE & 256 & 30 & Adam & \(10^{-5}\) & 4 & 100 \\ FAE & 1024 & 256 & Adam & \(10^{-5}\) & 4 & 200 \\ WAE & 1024 & 2 & Adam & \(10^{-7}\) & 4 & 200 \\ \end{tabular}
\end{table}
Table 1: Hyperparameters for the property autoencoders. Each autoencoder has an encoder and decoder neural network. These are MLPs, consisting of 7 fully connected layer with hyperbolic tangent activation functions.

Figure 4: The design of the core network. It consists of a one-to-many LSTM + MLP configuration that is run for 10 steps.
section in figure 3). The core network then predicts the latent representation of the evolved VULCAN output \(I_{o}\) by traversing the latent space in 10 steps (centre section in figure 3). The prediction of the latent representation of the VULCAN output is then decoded by the decoder part of the autoencoder to obtain the predicted mixing ratios (bottom right section in figure 3).

## 3 Results

### Autoencoders

#### 3.1.1 Mixing Ratio Autoencoder

The top row of figure 5 shows the reconstructed mixing ratio values against their actual value (left plot), for all examples from the test dataset. For the entire range of mixing ratios, the majority of the reconstructions lie within an order of magnitude of the diagonal line that marks perfect reconstructions, with an R-squared value of \(R^{2}=0.9997\).

The right plot of the top row of figure 5 shows the reconstruction error of the mixing ratios in logarithmic space. This scale is chosen because the autoencoders are trained in log space (see section 2.2). The solid line shows the median and the dashed lines show the 5th and 95th percentiles. From the right figure, we can see that 90% of the reconstructions have an error between -0.39 and 0.40 orders of magnitude.

#### 3.1.2 Flux Autoencoder

The middle row of figure 5 (left) shows the reconstructed flux values against their actual value. All reconstructed flux values are within 0.5 order magnitude of the graph diagonal. At fluxes with values around \(10^{6}\) erg m\({}^{-1}\) cm\({}^{-2}\) s\({}^{-1}\), the FAE is slightly underpredicting the actual flux values.

From the reconstruction error plot (right) we can see that 90% of the reconstructions have an error between -0.024 and 0.031 orders of magnitude. In this figure, we see a distinct underprediction of a small number of examples, which are the high flux values we see being underpredicted.

#### 3.1.3 Wavelength Autoencoder

The bottom row of figure 5 (left) shows the reconstructed wavelength values against their actual value. All reconstructed wavelength values are close to the actual values, deviating less than \(\sim 10\) nm from the graph diagonal. We can see that for wavelengths with values around 650 nm, the FAE has a tendency to underpredict.

From the reconstruction error of the wavelength values plot (right), we can see that 90% of the reconstructions have an error between -0.001 and 0.001 orders of magnitude. The slight underprediction of higher wavelength values is also visible in this figure.

#### 3.1.4 Pressure Profile Autoencoder

The top row of figure 6 (left) shows the reconstructed pressure values against their actual value. All reconstructed pressure values are well within 0.1 order magnitude of the graph diagonal.

From the reconstruction error plot of the pressure values (right) we can see that 90% of the reconstructions have an error between -0.0003 and 0.0003 orders of magnitude. This figure also shows a very minor underprediction of some pressure values, which correspond to pressure values of \(\sim 10^{3}\) bar.

#### 3.1.5 Temperature Profile Autoencoder

The middle row of figure 6 (left) shows the reconstructed temperature values against their actual value, for samples from the test dataset. It is immediately obvious that this autoencoder cannot accurately reconstruct the temperature profiles. A fraction of temperatures in the range \(\sim 750\) K \(<\) T \(\lesssim 1400\) K are reconstructed close to the graph diagonal, but the FAE is largely overpredicting temperatures below \(\sim 750\) K and underpredicting temperatures above \(\sim 750\) K.

The histogram of reconstruction errors of the temperature values (right) shows under- and over-prediction. 90% of the reconstructions have an error between 15.8 K and 754.97 K. Most predictions outside this range are underpredictions.

#### 3.1.6 Gravity Profile Autoencoder

The GAE shows similar behaviour to the FAE. The bottom row of figure 6 (left) shows the reconstructed gravity values against their actual value, for samples from the test dataset. Gravity values below \(\sim 4000\) cm s\({}^{-2}\) are underpredicted by the autoencoder, while values above \(\sim 4000\) cm s\({}^{-2}\) are underpredicted. Only values around \(\sim 4000\) cm s\({}^{-2}\) are predicted accurately by the GAE. In the reconstruction errors plot of the GAE of figure 6 (right) we can see that gravity values are consistently being over- and underpredicted. 90% of the reconstructions have an error between 197.46 and 4939.6 cm s\({}^{-2}\). Most predictions outside this range are over-predictions.

### Core Network

Because the FAE and GAE do not accurately reconstruct the temperature and gravity profiles, these profiles were not encoded for the final model. Instead, they were put directly in the latent representations of the inputs. This way, no information contained in these profiles is lost.

The left plot in figure 7 shows the mixing ratios predicted by the trained neural network model against the actual mixing ratios for the test dataset. The histogram shows that most of the model predictions lie within \(\sim 1\) order magnitude of the diagonal of the graph. A notable exception is the predictions for the few mixing ratios with values lower than \(\sim 10^{-44}\), for which the model overpredicts. These species can be neglected since they are not abundant enough to play a big role in the chemistry or to show features in the observable spectra.

Figure 7 (right) shows the mixing ratio prediction error of the neural network model, in log-space. The solid line shows the median and the dashed lines show the 5th and 95th percentiles. From the figure, we can see that 90% of the model predictions have an error between -0.66 and 0.65 orders magnitude. Outside of this range, the model does not show a clear tendency to either over- or underpredict.

Figure 8 shows selected examples (best, typical, and worst cases) of predictions by the neural network compared with the output of the VULCAN model, for a selection of seven species. The best case (top panel) shows a prediction that is almost indistinguishable from the actual mixing ratios. The examples in the typical case (middle panel) and worst case (lower panel) show larger prediction errors. In the typical case, CO\({}_{2}\), CO\({}_{3}\), and HCN have the largest prediction errors in the lower atmosphere, though still negligible. The worst case shows the largest prediction errors, with H having prediction errors of up to almost 1 order of magnitude in the upper atmosphere. Notable is that this case has very strong photochemistry in the upper atmosphere as it's positioned nearby its host star.

To compare the computational efficiency of VULCAN and the
Figure 5: The reconstructed against the actual input values (left column), and the reconstruction error in log space (right column) for the mixing ratios (top row), stellar flux (middle row), and the wavelengths (bottom row). The diagonal dashed line in the reconstructed vs. actual mixing ratios plot shows the performance of a perfectly reconstructing model. Here the colour represents the number of examples within each bin. In the reconstruction error figure, the solid line shows the median value, and the dashed lines show the 5th and 95th percentiles. The R\({}^{2}\) values of each reconstruction plot are shown in the left column.
Figure 6: The reconstructed against the actual input values (left column), and the reconstruction error in log space (right column) for the pressure profile (top row), temperature profile (middle row), and the gravity profile (bottom row). Note that the reconstruction errors plot of the temperature and gravity values are calculated in linear space.
neural network model, the computational time to calculate or predict every example in the full dataset was recorded. The results are presented in table 2. It should be noted that the VULCAN simulations were run on similar, but older hardware than the neural network model. The median computational times show a \(\sim 7.5\cdot 10^{3}\times\) decrease in computational time for the neural network model. The longest computational time required by the neural network model still shows a \(\sim 10^{3}\times\) decrease in computational time compared to the fastest VULCAN simulation.

## 4 Discussion

In this study, we successfully used autoencoders to extract most of the characterising input features and encode them into latent representations for the mixing ratios, stellar flux, wavelengths, and pressure profiles.Within these four groups, the largest prediction errors stem from the MRAE due to the high variability in input values, as opposed to the other input sources. We included initial and evolved mixing ratios of 69 species over 150 height layers. Additionally, the mixing ratio profiles among species differed significantly from one another (e.g. CH\({}_{4}\) and CO in figure 8). This made the complexity, of extracting and encoding the fundamental input features, highest for this particular autoencoder. In opposition, the variety in the flux from stellar spectra was much less. We obtained the stellar spectra from either the MUSCLES database or generated them using the PHOENIX model. The spectra from these different sources were quite distinct from each other in the EUV (0.5 - 200 nm). The PHOENIX models assume the spectra to follow blackbodies, while, in reality, M- and K stars have shown to be highly active in the EUV (Reiners and Basri (2008)), as was observed by the MUSCLES collaboration. Nonetheless, the spectra within each method seemed largely similar, which made it more straightforward for the FAE to learn how to accurately reconstruct them. For both the PAE and the WAE, the profiles were linearly spaced in logarithmic space, which made it easy for the autoencoders to learn how to encode these parameters. It is remarkable, however, that the WAE was not able to perfectly reproduce the wavelengths. A solution would be to make use of a handcrafted algorithm that encodes merely the first and last elements in the array. It is recommended to make use of such an algorithm for future use. Finally, the temperature- and gravity profile autoencoders were not successful at encoding and reconstructing their inputs. Both autoencoders produced the same solutions for each input example. The limited data set size and large variations in the temperature and gravity example cases could explain the autoencoders to be prone to errors. Future studies could focus on improving these specific autoencoders by performing root cause analysis. However, a more specialised approach to encoding the pressure and gravity profiles would be to provide hyperparameters, such as the planet mass and radius, directly to the core model. Such an approach negates the need to train autoencoders for these input parameters.

The prediction of the core network (LSTM) is within one order of magnitude for the majority (>90%) of the predictions. These errors are comparable with the discrepancies between different chemical kinetics codes (Venot et al., 2012). However, the accuracy of predictions of different examples varies. This inconsistency can arise due to some bias within the data set. Example cases similar to the best-case scenario (see 8) were more prevalent in the data set, causing the core network to produce better predictions of this type of hot-Jupiters. Additionally, by plotting the loss of each validation case against input parameters (see figure 9) it becomes apparent that some specific system parameters perform better than others. From figure 9, we see that planets with smaller orbit radii seem to have worse predictions. One explanation could be that these planets endure more irradiation from their host star, ensuring photochemistry being the dominant process in the upper atmosphere. The abrupt and severe

\begin{table}
\begin{tabular}{c c c c}
**code** & **median** & **minimum** & **maximum** \\ \hline \hline VULCAN & 5994.3 s & 1236.7 s & 102223.0 s \\ NN model & 0.77 s & 0.73 s & 0.93 s \\ \end{tabular}
\end{table}
Table 2: Median, minimum and maximum running times of VULCAN and the neural network model for all configurations in the dataset. VULCAN runs were performed on a single CPU core using an Intel(R) Xeon(R) CPU E5-4620 0 0 2.26GHz. The neural network model was run on a single core using an Intel(R) Xeon(R) W-1250 CPU @ 3.36GHz

Figure 7: The LSTM predicted mixing ratios plotted against the actual mixing ratios (left) and the LSTM mixing ratio prediction error in log space (right). The dashed diagonal line in the left plot shows the performance of a perfectly predicting model and the colour of each bin represents the number of predictions. The solid line in the right plot shows the median value, and the dashed lines show the 5th and 95th percentiles.
Figure 8: The mixing ratios per height layer for the best (top), typical (middle), and the worst (bottom) case of the validation set. The planet parameters for each case are given at the top of the plot. The solid lines show the actual mixing ratios as calculated by VULCAN, and the dashed lines show the neural network model predictions.
changes in abundances for some species due to photodissociation in the upper atmosphere could be difficult for the core network to learn with a limited dataset as provided in this study. Also noteworthy is the correlation between the planetary mass and the performance of the core network. Higher-mass planets tend to have lower losses as compared to lower-mass planets. Future work could focus on improving the prediction losses of the chemistry profiles for lower-mass planets and planets that orbit their host star close in.

We also showed that the trained model consistently over-predicts mixing ratios that have a value lower than \(10^{-44}\). This can again be explained by the lack of examples that have such lower values. Species with mixing ratios this low are small contributors to the atmospheric composition nevertheless and are not expected to affect forward models.

Finally, we want to note that the hyperparameters used in this study have been found by trial and error and have not been proven to be the most optimal values. Future studies could focus on a hyperparameter search for each individual autoencoder and core network to find the most optimal parameters.

Due to all mentioned caveats, there is room for improvement in future work. Here we detail some of the aspects that are out of the scope of this paper, but we will be looking into them in future publications. Evidently, a larger size of the data set, with which the network is trained, is expected to improve the results significantly. In order to train a neural network to be more generalised and less biased, a more diverse and extensive data set should be created. Free parameters that can be taken into account, which were not explored in this study, are variables such as the eddy diffusion coefficient 2, condensation, and composition of the atmosphere (e.g. varying metallicity and the C/O ratio).

Footnote 2: Note that vertical mixing is taken into account in every simulation, but is kept constant throughout the dataset.

Another approach to improve possibly the results is to change the model itself. The traditional autoencoders can be replaced with _variational autoencoders_, VAEs (Kingma and Welling, 2013). These types of autoencoders are based on Bayesian statistics. It is possible to regulate the latent space such that similar input examples have similar latent representations that lie close to each other within the latent space. The core network then might be able to learn how to traverse a regulated latent space and predict more accurately. The core network itself can be improved by, for example, including more time steps within the LSTM. Adding more time steps will ensure that the network predicts the solutions in a more similar way as VULCAN integrates toward the solution. A disadvantage of this is that the training time for the network will increase. Lastly, the recurrent neural network architecture could be changed to a _transformer_ design. Recently, the transformer neural network architecture, proposed by Vaswani et al. (2017), revolutionised the field of sequence transactions within machine learning. By using a so-called _attention_ mechanism, the transformer neural network outperforms recurrent neural networks in accuracy and efficiency. Because of the similarity between transformer and recurrent neural network applications, the core model may perform better when changing to this new type of architecture. However, because transformers require tokenized inputs, the autoencoders will also have to be changed to produce the expected outputs. The implementation of a transformer would therefore increase the complexity of the entire model and should be done carefully.

A different direction in approach would be, for example, to use interpolation methods in the already existing data set. The limitation of such methods is the ease of distribution. The size of the data set used in this study is 600GB, as opposed to a size of 3GB for the weights of the neural network used in this study.

## 5 Summary & Conclusions

In this study, we investigated the ability of a neural network to replace the time-dependent ordinary differential equations in the chemical kinetics code VULCAN (Tsai et al., 2017, 2021). The aim of this research was to explore the LSTM architecture for solving ordinary differential equations that include vertical mixing and photo-chemistry. We first created a data set that contains the in- and outputs of VULCAN simulations of hot-Jupiter atmospheres. We made use of the planetary mass (0.5 - 20 \(M_{J}\)), the semi-major axis (0.01 - 0.5 AU), and the stellar radius (1 - 1.5 \(R_{\odot}\)) as free parameters. Other parameters for the VULCAN configurations were derived either from analytical relations or kept constant throughout the data set. The input of the data set comprises the initial mixing ratios, the stellar spectrum, the temperature- and pressure profiles, and the gravity profiles. Note that the neural network trained in this study is limited to the chosen free parameters and can not be used for atmospheric models that include e.g. condensation. The outputs of the data set contain the mixing ratios of the species in the atmosphere, taken from 10-time steps (including the steady state) during the VULCAN simulation. This data was used to train a neural network that consists of two parts: the _autoencoder_ network and the _core_ network. The autoencoder was used to reduce the dimensionality of the input and output data from the data set by encoding them into lower dimensionality _latent representations_. The autoencoder network consisted of six smaller autoencoders, designed and trained to encode and decode the mixing ratios, flux, wavelengths, and temperature-, pressure-, and gravity profiles to and from their respective latent representations. The total input latent representation was the concatenation of these 6 smaller ones. The core network was designed to have an LSTM-based architecture and it mapped between the latent representation of the inputs to the encoded evolved output by traversing the _latent space_ in ten steps. During the training, the latent representations at these ten steps were compared to the ten sets of mixing ratios saved in the outputs of the data set to ensure that the core network is evolving the

Figure 9: The loss as a function of the semi-major axis of each validation case. The colour represents the planet mass in \(\mathbf{M_{J}}\) and the size of each scatter point represents the size of the host star which ranges between 1 \(M_{\odot}\) and 1.5 \(M_{\odot}\). The loss is calculated by making use of eq. 7
latent representation in a similar fashion as the VULCAN simulation evolves the mixing ratios. To summarise, we found that:

* the mixing ratios, flux, wavelengths, and pressure profile autoencoders were able to efficiently encode and accurately reconstruct their respective input properties
* the autoencoders were not able to encode and decode the temperature- and gravity profiles successfully. These autoencoders were, therefore, not used and instead, these profiles were put directly into the latent representation of the inputs
* the fully trained model (i.e. including the core network) was able to predict the mixing ratios of the species with the errors in the range [-0.66, 0.65] orders magnitude for 90% of the cases. Due to imbalances in the dataset, the model is biased to more accurately solve for some examples as compared to others
* the fully trained model is \(\sim 10^{3}\) times faster than the VULCAN simulations

Overall, this study has shown that machine learning is a suitable approach to accelerate chemical kinetics codes for modelling exoplanet atmospheres.

## Data Availability

All simulated data created in this study will be shared upon reasonable request to the corresponding author. The code and results are publicly available on github.com/JuliusHendrix/MRP.
