# DARWEN: Data-driven Algorithm for Reduction of Wide Exoplanetary Networks

A. Lira-Barria, J N Harvey, T. Konings, R. Baeyens, C. Henríquez, L. Decin, O. Venot, R. Veillet

► To cite this version:

A. Lira-Barria, J N Harvey, T. Konings, R. Baeyens, C. Henríquez, et al., DARWEN: Data-driven Algorithm for Reduction of Wide Exoplanetary Networks. Astronomy & Astrophysics - A&A, 2024, 692, pp.A158. 10.1051/0004-6361/202452070. [hal-05382040](https://hal.science/hal-05382040v1)

HAL Id: hal-05382040

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

Submitted on 25 Nov 2025

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>

Distributed under a Creative Commons Attribution 4.0 International License
# DARWEN: Data-driven Algorithm for Reduction of Wide Exoplanetary Networks

## An unbiased approach to accurately reducing chemical networks

A. Lira-Barria<fnref n="1" />,*, J. N. Harvey<fnref n="1" />, T. Konings<fnref n="2" />, R. Baeyens<fnref n="3" />, C. Henríquez<fnref n="4" />,
L. Decin<fnref n="2" />, O. Venot<fnref n="5" />, and R. Veillet<fnref n="5" />

Received 30 August 2024 / Accepted 18 November 2024

## ABSTRACT

*Context.* Exoplanet atmospheric modeling is advancing toward complex coupled circulation-chemistry models, from chemically diverse 1D models to 3D global circulation models (GCMs). These models are crucial for interpreting observations from facilities like JWST and ELT and understanding exoplanet atmospheres. However, maintaining chemical diversity in 1D models and especially in GCMs is computationally expensive, limiting their complexity. Optimizing the number of reactions and species in the simulated atmosphere can address this tradeoff, but there is a lack of transparent and efficient methods for this optimization in the current exoplanet literature.

*Aims.* We aim to develop a systematic approach for reducing chemical networks in exoplanetary atmospheres, balancing accuracy and computational efficiency. Our method is data-driven, meaning we do not manually add reactions or species. Instead, we test possible reduced chemical networks and select the optimal one based on metrics for accuracy and computational efficiency. Our approach can optimize a network for similar planets simultaneously, can assign weights to prioritize either accuracy or efficiency, and is applicable in the presence of photochemistry.

*Methods.* We propose an approach based on a sensitivity analysis of a typical 1D chemical kinetics model. Principal component analysis was applied to the obtained sensitivities. To achieve a fast and reliable reduction of chemical networks, we utilized a genetic algorithm (GA), a machine-learning optimization method that mimics natural selection to find solutions by evolving a population of candidate solutions.

*Results.* We present three distinct schemes tailored for different priorities: accuracy, computational efficiency, and adaptability to photochemistry. These schemes demonstrate improved performance and reduced computational costs. Our work represents the first reduction of a chemical network with photochemistry in exoplanet research.

*Conclusions.* Our GA-based method offers a versatile and efficient approach to reduce chemical networks in exoplanetary atmospheres, enhancing both accuracy and computational efficiency.

**Key words.** astrochemistry – methods: data analysis – planets and satellites: atmospheres – planets and satellites: composition – planets and satellites: gaseous planets

## 1. Introduction

The field of exoplanet atmospheric modeling is advancing toward more complex coupled circulation-chemistry models, including global circulation models (GCMs; Drummond et al. 2020; Ridgway et al. 2023; Zamyatina et al. 2023). These models are crucial for accurately interpreting observations from advanced facilities like the James Webb Space Telescope (JWST) and the Extremely Large Telescope (ELT). However, the computational expense of GCMs often limits the chemical diversity within these models, necessitating the development of reduced chemical networks (Venot et al. 2019, 2020).

Despite their simplicity and the challenge of capturing a planet's entire atmospheric chemistry in a single air column,

1D models have significantly advanced our understanding of hot exoplanets' atmospheric chemistry (Moses et al. 2011; Venot et al. 2012). These models reveal that a hot Jupiter's atmosphere has three distinct chemical regimes. In the lower atmosphere, high temperatures and pressures ensure chemical equilibrium. Ascending to the mid-atmosphere, vertical mixing outpaces chemical reactions, leading to a steady state distinct from chemical equilibrium that is invariant with altitude for the most abundant species. This transition from the lower to the middle atmosphere takes place at the quenching level, whose altitude varies depending on the species. Farther up, the highest part of the atmosphere is dominated by photochemistry due to stellar radiation. This framework, although a simplification, provides crucial insights into the dynamic processes governing atmospheric compositions across different height layers (Madhusudhan et al. 2016).

<footnote n="1">¹ Department of Chemistry, KU Leuven, Celestijnenlaan 200F, 3001 Leuven, Belgium</footnote>

<footnote n="2">² Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium</footnote>

<footnote n="3">³ Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH, Amsterdam, The Netherlands</footnote>

<footnote n="4">⁴ Millennium Institute Foundational Research on Data, Vicuña Mackenna 4860, Macul, Santiago, Chile</footnote>

<footnote n="5">⁵ Université Paris Cité and Univ Paris Est Creteil, CNRS, LISA, 75013 Paris, France</footnote>

<footnote>* Corresponding author; arturo.lira@kuleuven.be</footnote>
Observations of exoplanet atmospheres, using both ground-
based telescopes and space observatories, have detected key
molecules such as H₂O, CO₂, CO, and CH₄ (e.g., Kreidberg
et al. 2014; Xue et al. 2024; Alderson et al. 2023; Esparza-Borges
et al. 2023; Finnerty et al. 2024; Guilluy et al. 2019; Bell et al.
2023), alongside less frequently observed molecules like OH,
C₂H₂, HCN, and NH₃ (Nugroho et al. 2021; Landman et al.
2021; Brogi et al. 2023; Giacobbe et al. 2021; Cabot et al. 2019;
Basilicata et al. 2024). All these molecules are crucial for a
comprehensive understanding of exoplanet atmospheres. These
molecules are particularly significant as they are involved in piv-
otal thermal and chemical processes and are essential for precise
spectral analysis. Moreover, measuring the abundances of these
species plays a crucial role in estimating the carbon-to-oxygen
(C/O) ratio, an established marker for understanding planetary
formation and evolutionary history (Madhusudhan 2012). Recent
studies have also begun to explore nitrogen and sulfur ratios,
providing additional insights into planetary formation processes
(Kama et al. 2019; Turrini et al. 2021; Ohno & Fortney 2023;
Crossfield 2023). Together, these molecules serve as key indica-
tors not only for compositional analysis but also for the complex
interplay of chemical dynamics within exoplanet atmospheres.
While models ideally aim to be as complete as possible, priori-
tizing the accurate representation of this list of key species can
often suffice for many practical analytical purposes.

Water and carbon dioxide are common in hydrogen-
dominated, high-temperature environments (e.g., Rustamkulov
et al. 2023), while methane, ammonia, hydrogen cyanide, and
acetylene can indicate disequilibrium processes and the pres-
ence of hydrocarbons (Moses et al. 2011; Rimmer & Rugheimer
2019). The list of key molecules may expand in the future.
For instance, the recent recognition of SO₂ (Alderson et al.
2023; Dyrek et al. 2024; Tsai et al. 2023) underscores the
importance of considering disequilibrium chemistry, particu-
larly photochemistry, in atmospheric models (Tsai et al. 2024;
Baeyens et al. 2024). Observations show that key species
like SO₂ can drastically deviate from thermochemical equi-
librium and significantly impact observables, highlighting the
inadequacy of assuming chemical equilibrium when explain-
ing new data. Additionally, the number of studied exoplanet
atmospheres is rapidly increasing, enabling the simultaneous
detection of several atomic species (Prinoth et al. 2023; Mounzer
et al. 2022). This expansion has prompted further surveys
(Hoeijmakers et al. 2019; Lira-Barria et al. 2022) and the first
population studies (Welbanks et al. 2019; Changeat et al. 2022;
Edwards et al. 2023). As observational capabilities advance, the
need for comprehensive kinetic networks to accurately model
these atmospheric compositions in the presence of disequi-
librium becomes critical. The growing complexity of models
requires efficient and accurate methods for reducing the com-
putational expense of these more diverse chemical models,
ensuring they remain effective tools for interpreting the complex
dynamics and chemistry of exoplanet atmospheres.

To ensure accurate chemical modeling for exoplanets, Venot
et al. (2012, 2015, 2019, 2020) developed comprehensive chem-
ical networks applicable to various contexts, from hot Jupiters
to the Solar System's gas giants, and best suited for hot,
hydrogen-dominated atmospheres (which means they are not
well suited for temperatures below 300 K). These networks
were meticulously built and validated through comparison to
experimental combustion data (Veillet et al. 2024). In particular,

Venot et al. (2020) updated their full chemical scheme (referred to here as V20), which comprises approximately 100 species and 2000 reactions. Using the methodology introduced in Venot et al. (2019), they reduced it to a more compact form (referred to here as R20), consisting of 44 species and around 600 forward and reverse reactions. This careful reduction significantly lowered computational demands compared to the full model for a set of key molecules while preserving accuracy. The reduction was achieved through a detailed kinetic analysis using ANSYS© Chemkin-Pro Reaction Workbench and incorporating classical techniques from Lebedev et al. (2013), alongside expert assessment of the impact of removing steps and species in the reduced model.

Building on this previous work, we propose here a new and
automated method for generating optimized reduced models that
satisfy various criteria such as computational efficiency. We
achieved this by using a genetic algorithm (GA) to efficiently
reduce chemical networks. The GA simultaneously explores
multiple solutions, effectively bypassing local optima. This tech-
nique, operating analogously to natural selection, progresses
only the most promising solutions based on a predefined fitness
function, thus navigating through the most beneficial regions
of the solution space without needing derivative calculations
(Goldberg 1989). While GAs have been applied in chemical net-
work reduction (Edwards et al. 1998; Elliott et al. 2005; Sikalo
et al. 2014), including some astrochemical networks (Xu et al.
2019; Grassi et al. 2012, 2013), we are not aware of any instances
when they have been used in this way for modeling exoplanetary
atmospheres. In this work, we use V20 as a reference network
to evaluate our reduction method. We note that this network is
continually being improved (Veillet et al. 2024), but we do not
take the new updates into account here.

The efficacy of this method was tested on hot giant exo-
planets with substantial atmospheres, ideal candidates for atmo-
spheric characterization due to their numerous studies and
observable features. We specifically chose HD 209458b and HD
189733b as test cases; both planets have been subjects of com-
prehensive atmospheric composition studies (e.g., Redfield et al.
2008; Astudillo-Defru & Rojo 2013; Sánchez-López et al. 2019;
Lira-Barria et al. 2022; Xue et al. 2024; Blain et al. 2024).
These studies have examined a variety of atmospheric phenom-
ena, including dynamics, cloud formation (Moses et al. 2011;
Lines et al. 2018), thermal inversions (Diamond-Lowe et al.
2014), and atmospheric escape (Etangs et al. 2012; Lampón et al.
2020), establishing them as critical benchmarks for validating
new chemical reduction techniques.

In this paper we introduce our new GA-based method; it has
been designed to be both unbiased and adaptable and is capa-
ble of balancing accuracy and computational efficiency across
various planetary conditions, including those involving photo-
chemistry. In Sect. 2 we describe the 1D model we used and
the chemical network targeted for reduction. Section 3 details
our initial reduction technique, which involves a sensitivity
analysis and a reduction based on principal component anal-
ysis (PCA). The core component, our GA-based Data-driven
Algorithm for Reduction of Wide Exoplanetary Networks (DAR-
WEN), is explained in Sect. 4. In Sect. 5 we present three output
types tailored to different priorities – accuracy, computational
efficiency, and adaptability to photochemistry – discussing their
respective advantages and limitations. Finally, we conclude with
our findings and future directions in Sect. 6.
## 2. Modeling tools

### 2.1. 1D chemical kinetics and photochemistry model

We used and modified the 1D chemical kinetics code from Agúndez et al. (2014), which takes vertical mixing, molecular diffusion, and photochemistry into account, specifically utilizing its existing option to exclude horizontal mixing. Starting from initial abundances at chemical equilibrium at each pressure level, the model integrates the system of ordinary differential equations (ODEs) as a function of time. These ODEs are solved using the DLSODES (Hindmarsh 1983) Fortran library. This process continues until a steady state is reached, capturing the chemical evolution of the atmosphere over time. This steady state differs from chemical equilibrium because, while chemical reactions and transport processes balance each other out, the movement of molecules to different parts of the atmosphere causes the chemical composition to deviate from thermochemical equilibrium. The equation describing the formation and destruction of species is given by Eq. (1):

$$f_i t = \frac{P_i}{n} - f_i L_i - \frac{1}{nr^2} (r^2 \Phi_i) r. \quad (1)$$

In this equation, $f_i = n_i/n$ represents the molar fraction, also referred to as the abundance, of each species $i$, where $n_i$ is the species' number density, and $n$ is the total number density of all particles. The $f_i$ depends both on time ($t$) and the distance from the center of the planet ($r$), or equivalently the altitude ($h = r - r_0$); hence, a specific molar fraction is referred to as $f_{ih}^t$. The terms $P_i$ and $L_i$ denote the production and loss rates of species $i$, respectively, incorporating both chemical reactions and photochemical processes. The model also accounts for the vertical transport flux $\Phi_i$ of each species, which may be either positive (indicating upward movement) or negative (indicating downward movement), driven by atmospheric dynamics including eddy and molecular diffusion. This dual approach captures both larger-scale atmospheric movements driven by eddy diffusion and smaller-scale, molecular-level movements within the atmosphere (for details, see Agúndez et al. 2014; Bauer 1973).

### 2.2. Chemical network and other data

To properly model an atmosphere, the code needs several inputs. As the physical state of the atmosphere, we used the pressure-temperature and eddy diffusion profiles from the model grid of Baeyens et al. (2021, 2022). Specifically, we chose a model with ($T_{eff} = 1200$ K, $g = 1000$ cm s⁻², K5-star host) and a model with ($T_{eff} = 1400$ K, $g = 1000$ cm s⁻², G5-star host), as these parameters are representative of the hot Jupiters HD 189733 b and HD 209458 b, respectively. For initial elemental abundances, we assumed solar metallicity. The temperature profiles were derived from a 3D GCM (Carone et al. 2021) at the substellar point, with an isothermal upper atmosphere extension to a minimum pressure of $10^{-8}$ bar. Eddy diffusion coefficient profiles were self-consistently computed from the wind speeds (Baeyens et al. 2021).

Photochemical dissociations were computed based on the stellar UV spectra of a K5- and a G5-type star. For the former, we used the spectrum of HD 85512 from the MUSCLES survey (France et al. 2016; Youngblood et al. 2016; Loyd et al. 2016), while for the latter, we used a solar spectrum (Claire et al. 2012), following the approach in Baeyens et al. (2021). The photochemical cross sections and quantum yields used are based on laboratory data and ab initio calculations (Venot et al. 2012, 2020; Hébrard et al. 2013; Dobrijevic et al. 2014).

Finally, and key for this work, the code requires a chemical network to compute the production and loss rates of abundances at each time iteration. For this study, we used the comprehensive full V20 model from Venot et al. (2020), which we aimed to reduce. This network comprises 108 species (including C, H, N, and O atoms) with 1906 forward and reverse reactions. Appendix A explains how these reactions are treated in the code.

### 2.3. Convergence of the model

Our convergence values were measured using what we call the maximum relative deviation ($\Delta_{max}$):

$$\Delta_{\max}(f_{ih}^{n} > f_{min}) = \max_{i,h} \left( \frac{|f_{ih}^{n+1} - f_{ih}^{n}|}{f_{ih}^{n}} \right). \quad (2)$$

In this 1D model, $\Delta_{max}$ quantifies discrepancies between the abundances of species $i$ across atmospheric layers $h$ at two sequential time points, $f_{ih}^n$ and $f_{ih}^{n+1}$. For this purpose, we considered only significant abundances (i.e., those satisfying $f_{ih}^n > f_{min} = 10^{-20}$) to ensure focus on relevant chemical actors. The time points in the model are exponentially spaced after an initial integration time of one day, progressing to 10 days, 100 days, and so on. Chemical steady-state is considered reached once the $\Delta_{max}$ of Eq. (2) drops below the threshold of $10^{-4}$ for two successive integration periods<fnref n="1" />.

Although a stricter convergence to steady-state could be ideal in principle, the threshold used here is sufficient for conducting the sensitivity analysis and subsequent model reduction. Additionally, we observe that reaching a steady state was difficult when including the lowest layers of the atmosphere, with pressures higher than 10 bar or, in cases where photochemical processes are included in the model, when modeling high altitudes with pressures below $10^{-5}$ bar. As these pressure ranges do not qualitatively impact on the model outcomes, we therefore excluded them during the optimization process. More specifically, when developing models without photochemistry, we selected a range between 10 and $10^{-3}$ bar, encompassing the quench points. With photochemistry, we extended the range to $10^{-5}$ bar to include higher altitudes where photochemical processes are significant, while still achieving convergence to our chosen noise threshold.

## 3. Initial reduction technique

Our primary method is a GA, which requires an initial starting point to approximate the target area for finding a more efficient and accurate chemical scheme. To establish this starting point, we performed a sensitivity analysis (detailed in Sect. 3.1) to assess how variations in each reaction impact steady-state abundances. Based on these data, we ranked reactions using a PCA-based reduction process (detailed in Sect. 3.2), which helped us identify the most influential reactions for the subsequent optimization.

<footnote n="1">¹ The ODE solver's precision is set with a relative tolerance of $10^{-4}$, balancing computational efficiency and accuracy to match the target threshold. The absolute tolerance is set to an extremely small value ($10^{-99}$), ensuring it has no impact on the solution.</footnote>
## 3.1. Sensitivity analysis

A straightforward approach to ranking the importance of parameters in a model, such as reactions, is to perform a sensitivity analysis, either global (varying multiple parameters simultaneously) or local (varying one parameter at a time). While global sensitivity analyses are widely used to reduce models (Saltelli et al. 2005; Dobrijevic et al. 2010; Perumal et al. 2013; Nurislamova & Gubaydullin 2017; Till et al. 2019), it is well known that local sensitivity analysis is insufficient for ranking parameters or reducing models unless the model is linear, as it neglects intrinsic uncertainties. In this study we employed a local sensitivity analysis as an initial approach due to its simplicity and low computational cost; this provided a foundation for further refinement in subsequent steps.

For this analysis, we evaluated how slight variations in each rate constant, $k_j$, affect species abundances across different atmospheric layers. This process involves iteratively altering one rate constant at a time and observing the resulting changes in abundance of every species at each layer. We quantified these changes using sensitivity matrices, $S_j$, for each chemical reaction, $j$, in the model, whose elements measure the dependence of the concentration of a given species at a given altitude on the chosen rate constant.

The sensitivity coefficients, or elements of the matrix, were calculated based on a finite differences approximation to the derivative of the abundance $f_{ih}$ of species $i$ at altitude $h$, with respect to $k_j$. These differences were computed by re-propagating the model to steady-state with a modified $k_j$, obtained by multiplying it by a factor of 1.1, which is equivalent to the initial $k_j$ plus 10% of its value. This factor is large enough to induce a measurable change but small enough to be considered a finite difference.

We note that the finite-step perturbations of the rate constants are not intended here as an approach for uncertainty quantification. In fact, there will be different degrees of uncertainty for each reaction, and the impact of this rate constant inaccuracy is a separate issue that we do not tackle here. Instead, the 10% differential is used as the basis for ranking the importance of reactions, assuming that reactions causing larger changes are likely to be more significant for network reduction. Additionally, we normalized these sensitivity coefficients with respect to the reference value, $f_{ih}^0$, to ensure that both abundant and less abundant species are treated equitably. Each element of the sensitivity matrix $S_j$ is thus given by

$$S_{(i,h)}^j = \frac{1}{f_{ih}^0} f_{ih} k_j \sim \frac{\Delta f_{ih}}{f_{ih}^0 \Delta k_j} = \frac{f_{ih}^* - f_{ih}^0}{f_{ih}^0 (k_j^* - k_j^0)} \quad (3)$$

Here, $f_{ih}^0$ is the pre-variation steady-state abundance of species $i$ at each atmospheric layer $h$, while $f_{ih}^*$ represents the post-variation abundance upon reaching a new steady state. Each rate constant for the ~1000 forward reactions present in V20 is varied, and this results in a $S_j$ matrix that depends on the 107 species<fnref n="2" /> and 100 height layers for every one of them. If a reaction has an associated reverse reaction, its rate constant is changed by the same factor to preserve thermodynamic consistency (for details, see Appendix A).

## 3.2. Principal component analysis

After obtaining sensitivities for a given planet, we applied the PCA reduction method, as described by Lebedev et al. (2013).

Special emphasis was placed on including key molecules already observed in exoplanetary atmospheres (H₂O, CO₂, CO, HCN, C₂H₂, OH, CH₄, and NH₃) as they are the species we refer to throughout the text as the “key species” for evaluating our model’s performance. The steps of this method are the following:

1. We generated a correlation matrix based on the sensitivity coefficients to evaluate the impact of each reaction. This was done by using the overall sensitivity matrix, denoted as $S$. This matrix is a combination of individual matrices, $S_j$, each calculated according to Eq. (3). For any specific reaction $j$, we extracted the elements corresponding to the selected group of key molecules into a single row. These elements are indexed by a composite column label $(i, h)$. Consequently, the calculation of $SS^T$ associates each row with the relative importance of its corresponding reaction.

2. Eigenvalues and eigenvectors were calculated from this matrix, which could then be used to quantify the contribution to the model made by each reaction. The eigenvalues were utilized to filter the contribution of eigenvectors, isolating those with significant impact on chemical dynamics. For the filtered eigenvalues, $\lambda_i$, the corresponding eigenvectors, $u_i$, were projected into the reaction space, $r_j$, with the component of reaction $j$ in the new space given by $r_j^i = (u_i \cdot r_j)/\|u_i\|$.

3. We adjusted two thresholds in the PCA-reduced scheme specifically to achieve our target number of reactions and molecules. First, we selected eigenvalues, $\lambda_i$, that cumulatively represent 95% of the overall significance, which typically involves two or three principal components. Second, we set a predefined threshold for the reaction space components, $r_j^i$, based on the magnitude of sensitivities. By fine-tuning these thresholds, we intentionally limited the scheme to around 100 forward reactions, ensuring that we never included more than 50 species.

4. A species was selected for the reduced scheme only if it was involved in at least one of the qualifying reactions. Consequently, the selection of species and reactions for the final PCA scheme was constructed such that the already mentioned key molecules are present. From the 100 forward reactions mentioned in the previous step, the scheme includes approximately 40 species.

We discovered that while our schemes obtained by sensitivity analysis and subsequent PCA-based reduction proved remarkably effective, they yielded lower accuracy than R20, by about one order of magnitude. This discrepancy, we believe, stems from our sensitivity analysis being overly tuned to identify reactions with low reaction rates (see Appendix A). This approach leaves out crucial fast reactions, which in the full model precede or follow the most sensitive reactions we obtained, thereby yielding an incomplete network. However, this network serves as a suitable starting point for our core method.

## 4. A genetic algorithm for chemical network optimization

A GA is a machine learning optimization technique inspired by natural selection and genetics in biological systems. In a GA, each candidate solution is represented as a string of genes; specifically, in our work, each gene indicates whether a chemical reaction is included (1) or excluded (0) in a candidate chemical scheme<fnref n="3" />. GAs operate by iteratively selecting the fittest

<footnote n="2">² We did not take helium into account in this sensitivity analysis.</footnote>

<footnote n="3">³ Alternatively, each string can be seen as a list of reaction IDs included in the model.</footnote>
individuals in a population according to a loss function. At each iteration (or generation), each selected individual is used to generate multiple new individuals through random variation using the operators described below. Our initial population starts with a single progenitor derived from our PCA-reduced chemical scheme (detailed in Sect. 3.2).

Genetic algorithms are particularly useful for solving complex optimization problems where the solution space is vast and not easily navigable using traditional gradient-based methods (see Goldberg 1989, for a comprehensive description). According to the building blocks hypothesis (Holland 1992), the effectiveness of GAs arises from their ability to implicitly evaluate approximately $n_{pop}^3$ highly fit gene combinations, even though only $n_{pop}$ individuals are explicitly tested in each generation. This extensive evaluation occurs without additional computational or memory resources. However, a significant risk is premature convergence, where the population lacks diversity among gene combinations, leading to repetitive testing of limited solutions. This issue, described in Edwards et al. (1998), can be mitigated by maintaining sufficient genetic variation in each generation, thus ensuring continuous exploration of the solution space.

The core of our new reduction method is DARWEN, a GA designed to handle a vast space of $2^{N_{reacs}} \sim 2^{1000}$ possible combinations arising from the number of genes when trying to reduce the network described in Sect. 2.2. Although GAs are known to struggle in extremely large search spaces like $10^{1000}$ (Elliott et al. 2004), the space we are dealing with, while large, remains within a manageable range for our GA.

In our implementation, DARWEN adds and removes forward and reverse reactions together, which is why $N_{reacs}$ is 958 instead of 1906 – some reactions lack a reverse counterpart (see Table B.1 for all parameter values). Importantly, DARWEN does not test every possible reaction combination because schemes that fail to conserve mass within our 1D model are disqualified early in the process. Additionally, by limiting the addition of reactions to those associated with the initial species and a few additional molecules, we reduced the potential solution space to fewer than $2^{350}$ possibilities (the typical pool of available reactions is around 350).

At each generation, we selected $n_{fit}$ individuals as progenitors based on their fitness values calculated using the loss function $\varphi$. Every $n_{elit}$ generations, if the overall performance declines, we reintroduced progenitors from earlier generations to preserve valuable genetic information and enhance optimization robustness. Each progenitor has a probability of undergoing crossover ($p_c$), where segments of genes are exchanged with another progenitor, introducing variability and potentially beneficial combinations. There is also a probability of self-mutation ($p_m$), where genes are randomly flipped to include or exclude reactions.

To each individual, we applied a mutation rate, $\rho_m$, and a killing rate, $\rho_+$, to add or remove reactions, focusing on those associated with the initial species and a limited number of additional molecules. This targeted approach effectively reduced the solution space. To explore a broader solution space and address reactions that are consistently avoided, we introduced “curious” reactions, $n_{cur}$. These under-explored reactions helped ensure a broader search of the solution space. More details on these processes and the parameter values used are in Appendix B.

After an individual is created and meets the constraints on the number of molecules, we computed its loss function, $\varphi$, using the model detailed in Sect. 2.1. This loss function is the core of

DARWEN and employs $\Delta_{max}$ (Eq. (2)) to evaluate discrepancies between the reduced (○) and full (●) chemical schemes. It focuses on the key molecules mentioned earlier (see Sect. 3.2), as well as additional “major” species defined at each altitude as those with an abundance above a threshold ($f_{ih}^* > f_{min} = 10^{-20}$). This means that a molecule can be considered a major species at certain altitude layers but not at others, depending on its abundance. The loss function integrates those discrepancies to balance chemical accuracy against computational cost, which is roughly proportional to the number of molecules in each scheme ($n_{mlc}$). So if we consider

$$ \Delta_{max} \{\text{key}\} = \max_{\text{planets}} \left( \Delta_{\max\{\text{key}\}}^{(\circ, \bullet)} \right) \\ \Delta_{max} \{\text{maj}\} = \max_{\text{planets}} \left( \Delta_{\max\{f_{ih}^* > f_{min}\}}^{(\circ, \bullet)} \right), \qquad (4) $$

the loss function, $\varphi$, becomes

$$ \varphi = \Delta_{max} \{\text{key}\} + w_0 n_{mlc} + w_1 \log \left( \Delta_{max} \{\text{maj}\} \right). \qquad (5) $$

Weighting factors $w_0 = 0.1$ and $w_1 = 1$ are strategically set, meaning they are chosen based on careful consideration to balance the impact of species count and major species fitness within the loss function. These values were determined by aiming for desired results in the schemes. Specifically, if $w_0$ is too large, it can outweigh the other factors, making the model less effective. The inclusion of $w_1$ ensures that the model does not over-fit to just the key molecules. These factors remain adjustable to facilitate the development of faster models with relaxed accuracy requirements. While DARWEN effectively reduces the number of reactions for individual planets, this method aims to create a more general chemical scheme applicable to a group of planets with sufficiently similar chemistry. The algorithm manages multiple atmospheres by choosing the highest $\Delta_{max}$ among the planets being optimized for each chemical scheme at every generation. This approach ensures that when a planet with sufficiently similar characteristics is encountered, the existing reduced model can be applied directly, potentially reducing the need to rerun DARWEN.

As stated in Eq. (5), to optimize models for multiple planets simultaneously, we took the maximum value of each scheme metric when tested in the atmospheric models of different exoplanets. Our initial scheme is always based on a PCA-reduced chemical network. Initially, we attempted to combine the PCA-reduced schemes of our two test-case planets, HD 209458b and HD 189733b. However, this approach proved inefficient, especially when integrating photochemistry, as it introduced chemical inconsistencies. These issues arose from combining reactions from the different atmospheric environments of the planets. We found it more effective to use a PCA-reduced chemical scheme from just one planet. This approach avoided the aforementioned issues and provided a sufficient starting point, even when dealing with models without photochemistry.

## 5. Three types of reduced chemical schemes

We applied our method to generate new reduced versions of the V20 network for HD 209458b and HD 189733b's atmospheres. We employed three approaches: schemes specifically aimed at enhancing key molecule accuracy as in previous reduced chemical schemes (we refer to these as validation schemes); schemes
<table><caption>**Table 1.** Results for the two planets tested in this work, for each type of scheme.</caption><thead><tr><th rowspan="2">Scheme type</th><th rowspan="2">n<sub>mlc</sub></th><th rowspan="2">n<sub>reacs</sub></th><th colspan="4">HD 209458b</th><th colspan="4">HD 189733b</th></tr><tr><th>Δ<sup>(o,•)</sup><sub>max{key}</sub></th><th>Δ<sup>(o,•)</sup><sub>max{major}</sub></th><th>Δ<sup>(o,•)</sup><sub>max{all}</sub></th><th>Time [s]</th><th>Δ<sup>(o,•)</sup><sub>max{key}</sub></th><th>Δ<sup>(o,•)</sup><sub>max{major}</sub></th><th>Δ<sup>(o,•)</sup><sub>max{all}</sub></th><th>Time [s]</th></tr></thead><tbody><tr><td>V20</td><td>107</td><td>1904</td><td>—</td><td>—</td><td>—</td><td>500</td><td>—</td><td>—</td><td>—</td><td>650</td></tr><tr><td>R20</td><td>44</td><td>582</td><td>0.07</td><td>0.97</td><td>0.97</td><td>26</td><td>0.06</td><td>0.98</td><td>0.98</td><td>33</td></tr><tr><td>Validation</td><td>47</td><td>576</td><td>0.06</td><td>0.56</td><td>1.44</td><td>25</td><td>0.06</td><td>0.58</td><td>12.45</td><td>31</td></tr><tr><td>Low-cost</td><td>32</td><td>298</td><td>0.28</td><td>0.29</td><td>0.29</td><td>10</td><td>0.33</td><td>0.69</td><td>0.69</td><td>13</td></tr><tr><td>V20*</td><td>107</td><td>1956</td><td>—</td><td>—</td><td>—</td><td>900</td><td>—</td><td>—</td><td>—</td><td>950</td></tr><tr><td>Photoscheme</td><td>48</td><td>756</td><td>0.16</td><td>4.39</td><td>5.12</td><td>33</td><td>0.18</td><td>4.74</td><td>6.16</td><td>43</td></tr></tbody></table>

**Notes.** V20, R20 and V20* (full model by Venot et al. 2020 plus the corresponding photo reactions) are shown here for reference. *n*<sub>mlc</sub> and *n*<sub>reacs</sub> are the number of molecules and forward plus reverse reactions for each scheme. All schemes discrepancies are between the corresponding full model (•) and the reduced scheme in question (o). Time to reach steady state, according to our criteria (see Sect. 2.3), is measured via Terminal on a MacBook Pro 2022 with 16 GB of RAM and an M2 chip.

<figure><img src="image_2.png" /></figure>

**Fig. 1.** Example of DARWEN's "evolution" during simultaneous optimization of two planets. *Top*: minimum of the loss function, $\phi$, tracked across all schemes (or individuals) in each generation, showing the progress of optimization. *Middle*: proportion of individuals with a loss function below thresholds of 1.5, 1, and 0.75 per generation, providing insights into the population's overall performance distribution. This visual indicates how many schemes consistently achieve high fitness levels. *Bottom*: Heatmap visualizing the average value of each gene across all individuals per generation. If all schemes possess a value of 1 for a specific gene, the heatmap is orange, whereas for 0 it appears white. This helps us evaluate the diversity within the population and how closely the schemes align with the PCA-based initial model. It also allows us to identify which types of reactions are most prevalent during the optimization process.

that emphasize computational efficiency and thereby restrict the maximum number of species to 40 (low-cost schemes); and schemes that integrate photochemical processes, which is a new development in this field, and should provide optimal accuracy ("photoschemes").

To ensure the optimal evolution of DARWEN's algorithm, we experimented with several parameters, detailed in Table B.1. For the validation and low-cost schemes, we selected a pressure range from 10 to $10^{-3}$ bar, focusing on this interval for three reasons: it captures the entire chemically significant range when photochemistry is not considered, it aligns with the pressures probed by infrared observations, and it leads to efficient convergence of the Agúndez et al. (2014) model and, hence, rapid optimization. In contrast, for models involving photochemistry, we extended the range to 10 to $10^{-5}$ bar, which is more computationally demanding but essential for capturing the higher altitude regions where photochemical behavior is relevant. Additionally, we note that the inclusion of photochemistry requires either separate optimizations to account for the unique photolytic dynamics of each atmosphere or a less accurate reduced scheme. This is due to the increased computational expense of photochemical processes, which slows and complicates DARWEN's reduction process.

## 5.1. Validation scheme

We allowed DARWEN to evolve over 20 generations, transitioning from a PCA-reduced scheme to an optimal "validation" scheme, aiming for precision comparable to V20 for our specific planets. To obtain optimal schemes, we employed a strategic bottleneck approach: after approximately ten generations, we forced the best-performing scheme to become the initial scheme for a subsequent run of DARWEN (i.e., at generation ~10, we set $n_{fit} = 1$). An example of this process is shown in Fig. 1. Following this optimization, we identified around two dozen candidate schemes to validate our procedure. The final choice was made based on the accuracy of the broader set of the "major" (i.e., the most abundant) species. While R20 has 582 reactions (30% of the total number of reactions in the full model, and all the reactions available for that subset of 44 species), our validation scheme has 576 reactions (again 30% of the total number of reactions in the full model, but 85% of reactions available for those 47 species). As detailed in Table 1 the $\Delta_{max}$ for key molecules in the atmospheres of the two tested planets (0.05–0.06) is comparable to the $\Delta_{max}$ in the R20 scheme (0.07–0.06). Taking into account that an ideal case would have the full model abundances identical to the reduced model, resulting in a null $\Delta_{max}$, a lower $\Delta_{max}$ indicates a more accurate network reduction. This demonstrates DARWEN's capability to generate reduced chemical schemes, preserving accuracy with an error of less than 6% on key species.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Validation and low-cost schemes' performance on key and major molecules. The plots illustrate the performance of our method on key molecules for HD 209458b (a) and HD 189733b (b). *Left:* molar fractions of key molecules as a function of pressure in the atmosphere. The full V20 model is represented by a solid line, DARWEN's validation scheme by a dashed-dotted line, low-cost schemes by a dashed line, and the reduced scheme (R20) by a dotted line. *Right:* maximal percentage changes in the molar fractions of the key molecules compared to the full model. The full model's line is not shown because it serves as the reference model. Lines close to zero indicate that the predictions are nearly identical to those of the full model for that molecule. We only show changes bigger than 1%.</figcaption></figure>

Notably, the $\Delta_{\max}$ for major species in our model ($\Delta_{\max}\{maj\}$) was approximately 0.6, compared to around 1 in R20. Since major species generally outnumber key species, their abundances tend to show less precise agreement with those in the full model. While R20 maintained a similar level of discrepancies between all species and major species, we intentionally chose schemes that sacrificed accuracy for minor species to better preserve the accuracy of major species. This balance choice reflects our priorities, as we are generally less interested in minor species regarding their spectral features and overall impact on key species. By adjusting $w_1$ in our loss function $\phi$, we can control this balance; increasing $w_1$ elevates the relative importance of $\Delta_{\max}\{maj\}$, leading DARWEN to prefer schemes that are more precise for major species (for more details on important discrepancies on major species, see Appendix C).

Furthermore, our selected scheme includes more molecules but fewer reactions. While R20 is ~19 times faster than V20, our selected scheme was ~20 times faster than V20 (see Table 1). This efficiency, combined with improved accuracy for chemically significant species, demonstrates the effectiveness of our approach in optimizing chemical network reduction.

## 5.2. Low-cost schemes

We have confirmed that DARWEN is highly efficient at reducing schemes to a level of accuracy similar to R20. We then decided to explore how much we could reduce the number of molecules while maintaining this accuracy. Using our validation scheme as a starting point, we developed the “low-cost” scheme. As shown in Table 1, this scheme runs 2.5 times faster than R20 while maintaining comparable accuracy. The largest discrepancy with respect to V20, as shown in Fig. 2, occurs for HCN, which has an $\Delta_{\max}$ of 0.27 for HD 209458b and 0.33 for HD 189733b, though the shape of the abundance profiles as a function of altitude remains well described.

Achieving the presented results required approximately 30 generations, starting from our optimal validation scheme, to

balance accuracy for key molecules while minimizing the number of molecules and maintaining overall species accuracy. This published scheme has the lowest number of molecules while meeting our accuracy requirements, showcasing DARWEN's capability to efficiently optimize chemical network reduction, balancing precision and computational cost effectively.

## 5.3. Photoschemes

Our final and most challenging test involved incorporating photochemistry into the model reduction process, and it yielded highly promising results. Although DARWEN does not currently alter photochemical reactions directly, we applied photochemistry and allowed the 1D model to determine whether to include a photochemical reaction based on the presence of the involved species. The full model with photochemistry consists of 1906 “normal” reactions plus 52 photochemical reactions, totaling 1958 reactions (as shown in Table 1). The best scheme obtained by DARWEN, shown in Fig. 3, includes 48 molecules and 722 “normal” reactions (37% of the total number of normal reactions), and 34 photochemical reactions corresponding to the species present in the chemical scheme (65% of the total number of photochemical reactions).

As explained earlier, we ran DARWEN to optimize the test planets' atmospheres simultaneously. The resulting scheme attained accuracies of 0.16 and 0.18 for key molecules, with computation times of just 33 and 43 seconds for HD 209458b and HD 189733b, respectively. These computation times represent a substantial improvement over the full model, highlighting the practical advantages of our reduced schemes. Interestingly, HCN emerged as a major difference molecule again, alongside $CH_4$ and $CO_2$. The latter two molecules exhibited significant differences in the upper atmosphere, indicating a strong sensitivity to photochemistry, while HCN showed minor discrepancies in the lower layers, consistent with our previous low-cost schemes.

While the accuracy of the photoscheme was not as high as that of the validation scheme – being roughly two to three times lower – the abundances remain qualitatively well reproduced, as illustrated in Fig. 3. The discrepancies observed are still
<figure><img src="image_4.png" /><figcaption>Fig. 3. Same as Fig. 2 but for models in the presence of photochemistry, for HD 209458b (a) and HD 189733b (b).</figcaption></figure>

within acceptable ranges, especially considering that they may be smaller than current uncertainties arising from factors such as temperature, metallicity, and mixing rates, as well as uncertainties from unknown reactions, Arrhenius expression fitting of rate constants, and the eddy diffusion parameterization. This suggests that our reduced model with photochemistry is robust and reliable for practical applications.

Integrating photochemistry into the reduction process and achieving a consistent reduced scheme marks a significant advancement in this field, demonstrating DARWEN's ability to handle the complexities introduced by photochemical reactions. While the promising results indicate effectiveness, the observed lower accuracy highlights areas for further refinement, such as identifying and including critical missing reactions. Our data-driven approach is designed for iterative improvement, and future enhancements, including the incorporation of intrinsic model uncertainties, could lead to more reliable decision-making and models. Overall, this integration showcases DARWEN's strong potential, paving the way for more accurate and efficient chemical modeling of exoplanetary atmospheres.

## 6. Conclusions

We have developed DARWEN, a GA designed to efficiently reduce complex chemical networks for exoplanetary atmospheres. This method could potentially be applied to other atmospheric chemical contexts, such as those found in Solar System planets. We have demonstrated its ability to achieve accuracies similar to previous reduced models while significantly enhancing computational efficiency.

We applied DARWEN to reduce the full chemical model from Venot et al. (2020), V20, and compared our results to their associated reduced model, R20. We generated three types of reduced schemes: a validation scheme, a low-cost scheme, and a photoscheme, each with a different aim. The validation scheme focused on demonstrating DARWEN's capacity to enhance the accuracy of key molecules. The low-cost scheme aimed to improve computational efficiency by relaxing precision requirements. The photoscheme was designed to be the first reduced scheme in the field that includes photochemical reactions.

The validation scheme has a precision comparable to that of R20, with superior accuracy for the major species (maximum

discrepancy of 58% compared to 98% for R20). The low-cost scheme used the validation scheme as the initial solution and focused on minimizing the number of molecules while maintaining similar accuracy levels. This scheme ran 2.5 times faster than R20 and 10 times faster than V20, demonstrating DARWEN's capability to create more efficient schemes. Although some discrepancies were observed in HCN (with a maximum $\Delta$ of 33%), this is still within our acceptable uncertainty constraints. Our objective of obtaining a fast new scheme for a particular set of planets was achieved.

Our most challenging task was using DARWEN to develop the photoscheme, which involved incorporating photochemistry into the model reduction process. Although we achieved lower accuracies compared to our validation scheme and R20 – which cannot be applied when photochemistry is involved – the photoscheme represents a pioneering step in exoplanet research because it is the first time a chemical network that includes photochemical reactions has been reduced. The resulting schemes reached discrepancies of up to 18% in $\Delta$ for key molecules, with computation times of 33 and 43 seconds for HD 209458b and HD 189733b, respectively – more than 20 times faster than the full model. The challenges in optimizing the scheme for multiple planets, especially when incorporating photochemistry, highlight areas for further refinement.

Additionally, incorporating intrinsic model uncertainties into the decision-making process for evaluating discrepancies between the full and reduced models could enhance our approach. Considering that all full chemical schemes possess inherent limitations regarding the accuracy of their rate constants and other parameters derived from the literature, understanding the intrinsic uncertainties of the model is essential to assessing whether a discrepancy of 18% (or 33% in the case of the low-cost scheme) is significant. Nevertheless, we are optimistic about these results, as they mark a significant advancement in incorporating photochemistry into reduced chemical networks at a time when photochemical processes are increasingly important in the characterization of exoplanets.

While DARWEN has proven effective in reducing chemical networks for the specific exoplanets we have studied, applying it to other planets currently requires conducting a new run of the algorithm with the necessary planetary data for each case. This means that for detailed models – such as those needed for 3D atmospheric simulations – we must input the relevant data into
DARWEN to generate a new, tailored chemical scheme for each planet or group of planets of interest. In future work, we plan to expand this approach by applying DARWEN to a selection of exoplanets observed or expected to be observed with JWST and other facilities. By examining how the resulting chemistry of these planets varies, we aim to gain deeper insight into the diversity of exoplanetary atmospheres and, hopefully, construct a more comprehensive chemical scheme that can be applied more broadly.

Overall, DARWEN is a powerful tool for optimizing chemical network reductions, effectively balancing precision and computational efficiency. Its application to exoplanetary atmospheres opens new avenues for exploring and understanding the complex chemistry of these distant worlds. By developing these reduced networks, we enable the inclusion of detailed kinetic networks in complex models such as GCMs, thereby enhancing our understanding of exoplanets.

*Acknowledgements.* O.V. and R.V. acknowledge funding from Agence Nationale de la Recherche (ANR), project "EXACT" (ANR-21-CE49-0008-01). In addition, O.V. acknowledges funding from the Centre National d'Études Spatiales (CNES). A.L.B., J.N.H., T.K., R.B. and L.D. acknowledge KULeuven funding of the ID-N ESCHER project.

## References

Agúndez, M., Parmentier, V., Venot, O., Hersant, F., & Selsis, F. 2014, A&A, 564, A73

Alderson, L., Wakeford, H. R., Alam, M. K., et al. 2023, Nature, 614, 664

Astudillo-Defru, N., & Rojo, P. 2013, A&A, 557, A56

Baeyens, R., Decin, L., Carone, L., et al. 2021, MNRAS, 505, 5603

Baeyens, R., Konings, T., Venot, O., Carone, L., & Decin, L. 2022, MNRAS, 512, 4877

Baeyens, R., Désert, J.-M., Petrignani, A., Carone, L., & Schneider, A. D. 2024, A&A, 686, A24

Basilicata, M., Giacobbe, P., Bonomo, A., et al. 2024, A&A, 686, A127

Bauer, S. J. 1973, Physics of Planetary Ionospheres (New York: Springer-Verlag), 6

Bell, T. J., Welbanks, L., Schlawin, E., et al. 2023, Nature, 623, 709

Blain, D., Sánchez-López, A., & Mollière, P. 2024, AJ, 167, 179

Brogi, M., Emeka-Okafor, V., Line, M. R., et al. 2023, AJ, 165, 91

Cabot, S. H., Madhusudhan, N., Hawker, G. A., & Gandhi, S. 2019, MNRAS, 482, 4422

Carone, L., Mollière, P., Zhou, Y., et al. 2021, A&A, 646, A168

Changeat, Q., Edwards, B., Al-Refaie, A. F., et al. 2022, ApJS, 260, 3

Claire, M. W., Sheets, J., Cohen, M., et al. 2012, ApJ, 757, 6

Crossfield, I. J. M. 2023, ApJ, 952, L18

Diamond-Lowe, H., Stevenson, K., Bean, J., Line, M., & Fortney, J. 2014, ApJ, 796, 66

Dobrijevic, M., Hébrard, E., Plessis, S., et al. 2010, Adv. Space Res., 45, 77

Dobrijevic, M., Hébrard, E., Loison, J., & Hickson, K. 2014, Icarus, 228, 324

Drummond, B., Hébrard, E., Mayne, N. J., et al. 2020, A&A, 636, A68

Dyrek, A., Min, M., Decin, L., et al. 2024, Nature, 625, 51

Edwards, K., Edgar, T. F., & Manousiouthakis, V. I. 1998, Comp. Chem. Engng., 22, 239

Edwards, B., Changeat, Q., Tsiaras, A., et al. 2023, ApJS, 269, 31

Elliott, L., Ingham, D. B., Kyne, A. G., et al. 2004, Prog. Energy Combus. Sci., 30, 297

Elliott, L., Ingham, D. B., Kyne, A. G., et al. 2005, Industrial Eng. Chem. Res., 44, 658

Esparza-Borges, E., López-Morales, M., Redai, J. I. A., et al. 2023, ApJ, 955, L19

Etangs, A. L. D., Bourrier, V., Wheatley, P. J., et al. 2012, A&A, 543, L4

Finnerty, L., Xuan, J. W., Xin, Y., et al. 2024, AJ, 167, 43

France, K., Loyd, R. O. P., Youngblood, A., et al. 2016, ApJ, 820, 89

Giacobbe, P., Brogi, M., Gandhi, S., et al. 2021, Nature, 592, 205

Goldberg, D. E. 1989, Genetic Algorithms in Search, Optimization, and Machine Learning (New York: Addison-Wesley)

Grassi, T., Bovino, S., Gianturco, F. A., Baiocchi, P., & Merlin, E. 2012, MNRAS, 425, 1332

Grassi, T., Bovino, S., Schleicher, D., & Gianturco, F. A. 2013, MNRAS, 431, 1659

Guilluy, G., Sozzetti, A., Brogi, M., et al. 2019, A&A, 625, A107

Hindmarsh, A. C. 1983, Sci. Comput., 1, 55

Hoeijmakers, J., Ehrenreich, D., Kitzmann, D., et al. 2019, in EPSC-DPS Joint Meeting 2019, 2019, EPSC-DPS2019-1405

Holland, J. H. 1992, Sci. Am., 267, 66

Hébrard, E., Dobrijevic, M., Loison, J. C., et al. 2013, A&A, 552, A132

Kama, M., Shorttle, O., Jermyn, A. S., et al. 2019, ApJ, 885, 114

Kreidberg, L., Bean, J. L., Désert, J. M., et al. 2014, ApJ, 793, L27

Lampón, M., López-Puertas, M., Lara, L. M., et al. 2020, A&A, 636, A13

Landman, R., Sánchez-López, A., Mollière, P., et al. 2021, A&A, 656, A119

Lebedev, A. V., Okun, M. V., Chorkov, V. A., Tokar, P. M., & Strelkova, M. 2013, J. Math. Chem., 51, 824

Lines, S., Mayne, N. J., Boutle, I. A., et al. 2018, A&A, 615, A97

Lira-Barria, A., Rojo, P. M., & Mendez, R. A. 2022, A&A, 657, A36

Loyd, R. O. P., France, K., Youngblood, A., et al. 2016, ApJ, 824, 102

Madhusudhan, N. 2012, ApJ, 758, 36

Madhusudhan, N., Agúndez, M., Moses, J. I., & Hu, Y. 2016, Space Sci. Rev., 205, 285

Moses, J. I., Visscher, C., Fortney, J. J., et al. 2011, ApJ, 737, 37

Mounzer, D., Lovis, C., Seidel, J. V., et al. 2022, A&A, 668, A1

Nugroho, S. K., Kawahara, H., Gibson, N. P., et al. 2021, ApJ, 910, L9

Nurislamova, L. F., & Gubaydullin, I. M. 2017, J. Math. Chem., 55, 1779

Ohno, K., & Fortney, J. J. 2023, ApJ, 956, 125

Perumal, T. M., Krishna, S. M., Tallam, S. S., & Gunawan, R. 2013, Comp. Chem. Eng., 56, 37

Prinoth, B., Hoeijmakers, H. J., Pelletier, S., et al. 2023, A&A, 678, A182

Redfield, S., Endl, M., Cochran, W. D., & Koesterke, L. 2008, ApJ, 673, L87

Ridgway, R. J., Zamyatina, M., Mayne, N. J., et al. 2023, MNRAS, 518, 2472

Rimmer, P. B., & Rugheimer, S. 2019, Icarus, 329, 124

Rustamkulov, Z., Sing, D. K., Mukherjee, S., et al. 2023, Nature, 614, 659

Saltelli, A., Ratto, M., Tarantola, S., & Campolongo, F. 2005, Chem. Rev., 105, 2811

Sikalo, N., Hasemann, O., Schulz, C., Kempf, A., & Wlokas, I. 2014, Int. J. J. Kinetics, 46, 41

Sánchez-López, A., Alonso-Floriano, F. J., López-Puertas, M., et al. 2019, A&A, 630, A53

Till, Z., Varga, T., Sója, J., Miskolczi, N., & Chován, T. 2019, Chem. Eng. J., 375, 121920

Tsai, S. M., Lee, E. K., Powell, D., et al. 2023, Nature, 617, 483

Tsai, S.-M., Innes, H., Wogan, N. F., & Schwierterman, E. W. 2024, ApJ, 966, L24

Turrini, D., Schisano, E., Fonte, S., et al. 2021, ApJ, 909, 40

Veillet, R., Venot, O., Sirjean, B., et al. 2024, A&A, 682, A52

Venot, O., Hébrard, E., Agúndez, M., et al. 2012, A&A, 546, A43

Venot, O., Hébrard, E., Agúndez, M., Decin, L., & Bounaceur, R. 2015, A&A, 577, A33

Venot, O., Bounaceur, R., Dobrijevic, M., et al. 2019, A&A, 624, A58

Venot, O., Cavalié, T., Bounaceur, R., et al. 2020, A&A, 634, A78

Welbanks, L., Madhusudhan, N., Allard, N. F., et al. 2019, ApJ, 887, L20

Xu, R., Bai, X.-N., Öberg, K., & Zhang, H. 2019, ApJ, 872, 107

Xue, Q., Bean, J. L., Zhang, M., et al. 2024, ApJ, 963, L5

Youngblood, A., France, K., Loyd, R. O. P., et al. 2016, ApJ, 824, 101

Zamyatina, M., Hebrard, E., Drummond, B., et al. 2023, MNRAS, 519, 3129
# Appendix A: Sensitivity analysis and rate constant variations

The code receives coefficients to build the rate constants according to the modified Arrhenius equation:

$$k_f(T) = \alpha T^{\beta} \exp\left(\frac{-\gamma}{T}\right), \qquad (A.1)$$

where $k_f$ is the forward rate constant, $T$ is the temperature, $\gamma$ is the activation energy divided by the Boltzmann constant, and $\alpha$ and $\beta$ are coefficients. To obtain the reverse reaction rate constant $k_r$ of a given reaction, we used the fact that the constant of equilibrium, $K_{eq}$, is related to the standard Gibbs energy, $\Delta G^0$, and standard pressure, $P^0$,

$$K_{eq} = \frac{k_f}{k_r} = \exp\left(\frac{-\Delta G^0}{RT}\right) \left(\frac{P^0}{k_B T}\right)^{\Delta\nu}, \qquad (A.2)$$

where $k_B$ is the Boltzmann constant and $R$ is the Rydberg constant.

To assess the sensitivity of reactions that can be modeled as Eq. A.1 we varied their alpha coefficient, $\alpha$, such that $\alpha^* = 1.1\alpha$. $\beta$ and $\gamma$ are unchanged. The choice of $\Delta\alpha = 0.1\alpha$ strikes a balance between ensuring the alteration is sufficiently small to approximate a derivative, yet large enough to noticeably affect species abundances given the achievable steady-state convergence threshold. Also, to maintain the thermodynamics unvaried, for each forward reaction varied we varied the corresponding reverse reaction rate constant parameter by the same factor. We note that for some reactions, the rate constant as a function of temperature is expressed in a different way (see Venot et al. 2012). For all reactions, we used the same multiplicative perturbation approach so that $k_f$ is multiplied by 1.1 at each temperature.

A previous study advises against using variations of more than 5% for “local sensitivity analysis” (Saltelli et al. 2005; Dobrijevic et al. 2010) for related sensitivity analyses, but we have found only modest changes in terms of predicted model reduction properties between our chosen 10% variations and a smaller 5% change, with the molecules and reactions selected in the PCA-based reduced scheme being only slightly changed. Further changes following the GA component of the algorithm will still further reduce these differences.

# Appendix B: Genetic algorithm details

The iterative process at each generation involves the following steps:

1. **Selection:** At each generation, $n_{fit}$ schemes are chosen as progenitors. This can be done deterministically by selecting the $n_{fit}$ fittest schemes or using a weighted roulette rule, where selection probability is $\varphi_s / \sum_s \varphi_s$. In the first iteration, $n_{fit} = 1$.

2. **Crossover:** Each progenitor has a probability $p_c$ of crossover. We applied a one-point crossover by selecting a random position in the binary gene sequence and swapping genes beyond this point with another progenitor.

3. **Self-mutation:** Each progenitor has a mutation probability $p_m$. About $p_m \times n_{pop}$ schemes mutate per generation, flipping a random number of genes (up to $p_m \times N_{reacs}$). Mutations may introduce new reactions, adding the necessary molecules to the species list, increasing diversity.

4. **Exploring the species pool:** Each progenitor generates $n_{child}$ new schemes. Reaction inclusion is guided by mutation rate $\rho_m$ and restricted reaction pool, with the killing rate $\rho_\dagger$ controlling how many reactions are excluded. Our code sets species number limits, giving each new individual an equal chance to fall within this range. Molecules are randomly removed if fewer species are generated, or added if more are needed. Key species are always retained to balance computational cost and accuracy.

5. **Performance evaluation:** All individuals are evaluated according to Eq. 5.

6. **Elitism:** Every $n_{elit}$ cycles, if performance drops, the $n_{fit}$ fittest individuals from previous generations are selected, preserving valuable information.

7. **Curiosity mechanism:** To explore neglected regions of the solution space, we rewarded the inclusion of under-explored reactions (i.e., those added fewer than $N_{cur}$ times to fit schemes). Optionally, $n_{cur}$ such reactions can be added to the initial scheme, increasing variability and promoting thorough exploration. This optional step, original to our work, helps increase the variability of schemata and ensures the exploration of neglected solution spaces.

# Appendix C: Major species discrepancies

In the main text, we mention that our criterion for choosing between different promising schemes identified by our GA-based approach relates to major species discrepancies. Here, we present the species that exhibit a discrepancy of over 10% at least at one height of the model, as shown in Fig. C.1. This criterion for our final decision is optional; we could have chosen the scheme with the smallest $\Delta$. However, $\Delta$ did not seem to correlate well with major species abundances beyond a certain extent. For example, Fig. C.1 shows that some major species in R20 have discrepancies of almost 100%. The reasons why these discrepancies do not correlate and how they compare to the uncertainties in the model require follow-up studies.
<figure><img src="image_5.png" /><figcaption>**Fig. A.1.** Sensitivity analysis of HD 209458b without photochemistry: *Top*: Distribution of sensitivities for all molecules at all heights for each reaction ID (j), showing only those above the expected ODE solver noise level. Horizontal lines indicate the noise level and 10 % variation, the maximum expected if the response is linear. *Bottom*: Count of sensitivities above the noise level for each reaction.</figcaption></figure>

<table><caption>Table B.1: Glossary of terms related to our GA.</caption><thead><tr><td>Term</td><td>Value</td><td>Definition</td></tr></thead><tbody><tr><td>N<sub>reacs</sub></td><td>958</td><td>Total number of forward reactions in our tested network</td></tr><tr><td>n<sub>reacs</sub></td><td>~ 150 - 600</td><td>Total number of forward and reverse reactions tested per scheme</td></tr><tr><td>n<sub>fit</sub></td><td>15-20</td><td>Number of progenitors selected per generation</td></tr><tr><td>n<sub>child</sub></td><td>15-20</td><td>Number of children generated from each progenitor</td></tr><tr><td>n<sub>pop</sub></td><td>n<sub>fit</sub> × n<sub>child</sub></td><td>Total number of individuals per generation</td></tr><tr><td>p<sub>c</sub></td><td>0.5</td><td>Probability of crossover</td></tr><tr><td>p<sub>m</sub></td><td>0.05</td><td>Probability of self-mutation</td></tr><tr><td>ρ<sub>m</sub></td><td>~0.15</td><td>Mutation rate, determining how many reactions from the pool are added.</td></tr><tr><td>ρ<sub>‡</sub></td><td>~ 0.1</td><td>Killing rate, determining how many reactions from the individual are excluded.</td></tr><tr><td>n<sub>mlc</sub></td><td>~ 35-55</td><td>Number of molecules per scheme</td></tr><tr><td>n<sub>elit</sub></td><td>~ 7</td><td>Frequency of elitism application (generations)</td></tr><tr><td>n<sub>cur</sub></td><td>≤ 5</td><td>Number of curious reactions added per scheme</td></tr><tr><td>N<sub>cur</sub></td><td>100</td><td>Upper threshold for defining an under-explored reaction</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Fig. C.1.** Discrepancies on selected major species. Major species with discrepancies of more than 10% compared to the full model at any height are shown for (a) HD 209458b (without photochemistry), (b) HD 189733b (without photochemistry), (c) HD 209458b (with photochemistry), and (d) HD 189733b (with photochemistry). **Top panels:** Absolute abundances. **Bottom left:** Vertical profile of $\Delta_{\max}$ for each subset of species. **Bottom right:** Maximum discrepancy for those species at each height. Note that in HD 189733b with photochemistry (d), the maximum discrepancy does not always map the species present because $N_2H_3$ is defined as a major species only at certain pressure levels.</figcaption></figure>