Automated Chemical Reaction Network Generation and Its Application to Exoplanet Atmospheres

With the advent of JWST and the spectroscopic characterization of exoplanet atmospheres in unprecedented detail, there is a demand for more complete pictures of chemical and photochemical reactions and their impacts on atmospheric composition. Traditionally, building reaction networks for (exo)planetary atmospheres involves manually tracking relevant species and reactions, a time-consuming and error-prone process. This approach’s applicability is also often limited to specific conditions, making it less versatile for different planetary types (i.e., photochemical networks for Jupiters may not be directly applicable to water-rich exoplanets). We introduce an automated approach using a computer-aided chemical reaction network generator, combined with a 1D photochemical kinetic-transport model, offering significant advantages. This approach automatically selects reaction rates through a rate-based iterative algorithm and multiple refinement steps, enhancing model reliability. Also, this approach allows for the efficient simulation of diverse chemical environments, from hydrogen to water, carbon dioxide, and nitrogen-dominated atmospheres. Using WASP-39b and WASP-80b as examples, we demonstrate our approach’s effectiveness, showing good agreement with recent JWST data. Our WASP-39b model aligns with prior studies and JWST observations, capturing photochemically produced sulfur dioxide. The WASP-80b model reveals an atmosphere influenced by deep-interior thermochemistry and vertical mixing, consistent with JWST NIRCam observations. Furthermore, our model identifies a novel initial step for the N2–NH3–HCN pathway that enhances the efficiency of the conversion in high-temperature/high-pressure environments. This automated chemical network generation offers a novel, efficient, and precise framework for studying exoplanetary atmospheres, marking a significant advancement over traditional modeling techniques.


INTRODUCTION
Our knowledge of other stellar systems and their accompanying planets has been expanding significantly since the exoplanet surveys of Kepler, K2, and Transiting Exoplanet Survey Satellite (TESS) satellites.Over 5500 exoplanets have been confirmed (NASA's Exoplanet Archive 2023).Adding to this, the recent launch of JWST has provided us with a deluge of high-quality Corresponding author: Jeehyun Yang; Renyu Hu jeehyun.yang@jpl.nasa.gov;renyu.hu@jpl.nasa.govspectroscopic data.This allows for the characterization of exoplanet atmospheres with unprecedented detail, exemplified by the detection of SO 2 in the hot-Jupiter WASP-39 b's atmosphere, which indicates active photochemistry (Tsai et al. 2023).Another example is the detection of CO 2 in the temperate sub-Neptune K2-18 b's atmosphere (Madhusudhan et al. 2023), which supports the hypothesis of a water-rich interior (Madhusudhan et al. 2021;Hu et al. 2021), but can also be potentially explained by a high-metallicity atmosphere (Yu et al. 2021;Hu et al. 2021;Tsai et al. 2021;Wogan et al. 2024).
As shown above, JWST enables detailed atmospheric measurements of diverse types of exoplanets from Jupiter-sized to Earth-sized, from cool to hot atmospheres, on which one may expect diverse atmospheric Yang & Hu composition and redox conditions.1-D photochemical atmospheric modeling is crucial for interpreting the JWST observations as well as guiding future observations.Thus, enhancing 1-D photochemical atmospheric modeling with more comprehensive chemical and photochemical reaction networks enables more precise characterization of exoplanet atmospheres, guiding future observations, and advancing our understanding of exoplanets.
Multiple photochemical reaction networks have already been developed for exoplanet atmospheric modeling studies (e.g., Moses et al. 2011;Hu et al. 2012;Venot et al. 2012;Tsai et al. 2017;Rimmer & Rugheimer 2019).However, there still have been several limitations to constructing these chemical reaction networks applicable to exoplanet atmospheric conditions.One of the major issues is the way photochemical networks are constructed.Most of these existing chemical networks that describe various exoplanet atmospheres are constructed by hand, adopting reaction rate and thermodynamic parameters by carefully keeping track of all possible species and reactions relevant to the target system (e.g., Moses et al. 2011;Hu et al. 2012;Venot et al. 2012;Tsai et al. 2017;Rimmer & Rugheimer 2019;Tsai et al. 2021).This process is very time-consuming and error-prone, and the resulting model significantly depends on the chemistry knowledge of the builder who manually chooses these parameters from various sources (i.e., laboratory measurements, ab initio calculations, estimations, etc.).For this reason, the chance that the number of missing and dubious reactions are included in the model increases as the model size grows, eventually leading to a failure in precisely predicting and interpreting important reaction species and pathways.
A recent study by Veillet et al. (2023) constructed C-H-O-N chemical networks based on an extensive amount of combustion experimental data gathered over recent decades.This provides a relatively robust chemical network for describing the atmospheres of hot Jupiter exoplanets, predominantly composed of H 2 with an insignificant amount of sulfur species (since the network doesn't contain sulfur-bearing species).However, the applicability of such a model to other types of planets is limited, often necessitating significant time and effort to develop another model for another system.This limitation arises because the relevance of specific chemical species and reactions is intrinsically tied to the system's conditions, such as temperature, pressure, and the dominant gas species.As a result, the chemical network built by Veillet et al. (2023), while built with an intent to model H 2 -dominated atmospheres only, is not applicable for, e.g., Venus-like exoplanet atmospheres, given the differences in temperature and pressure profiles, as well as in the dominant atmospheric gas composition (e.g., CO 2 -dominated atmosphere with sulfurbearing species).While some might suggest including all known species and reactions studied so far, doing so is impractical.The computational time required for large chemical kinetics simulations scales approximately linearly with the number of chemical reactions and approximately quadratically with the number of chemical species, N (Schwer et al. 2002).Given the vast amount of spectroscopic data expected from the JWST and future observational missions, a fundamentally new approach to photochemical reaction network construction is essential.
Over the past decade, advancements in computational chemical engineering have paved the way for automated chemical reaction network generation.These automation techniques can be categorized based on their approaches to species and reaction selection, as well as parameter generation.One common method is to define reaction families to find possible reactions, which allows for the expansion of the network starting from an initial set of molecules (e.g., Sarathy et al. 2012).Another approach is a rule-based method, as adopted by packages like Genesys (e.g., Vandewiele et al. 2012).The inclusion of species and reactions in the network is determined by a set of user-defined constraints.Although computationally efficient, these constraints, often dependent on the user's chemistry knowledge, can potentially bias network generation.In contrast, the Reaction Mechanism Generator (RMG, Gao et al. (2016); Liu et al. (2021); Johnson et al. (2022)) employs a ratebased method, where the importance of a species or reaction is determined based on iterative simulations of the chemical system.For this reason, this rate-based approach is more objective (i.e., independent of the user's chemistry knowledge) and can provide a more comprehensive chemical network than other methods.The only downside of this rate-based method is that it is computationally more expensive than other methods.These RMG-generated networks have been actively utilized in various chemical engineering fields, such as a computergenerated acetylene pyrolysis model by Liu et al. (2020), which successfully described the previous experiment by Norinaga et al. (2008), butyl acetate pyrolysis and combustion model by Dong et al. (2023), and methyl propyl ether pyrolysis and oxidation model by Johnson et al. (2021) showing excellent agreement with most of the shock tube and rapid compression machine data.Notably, this automation has recently been used to simulate and successfully rationalize the previous laboratory photochemical studies by Fleury et al. (2019Fleury et al. ( , 2020) ) that simulated hot Jupiter-like atmospheric conditions (Yang et al. 2023).Such applications underscore the reliability and vast potential of automated reaction network generators, particularly when used together with existing photochemical reaction networks, offering solutions to challenges associated with manual methodologies.
Given the challenges and potential of recent advancements in computational chemical engineering, we develop the atmospheric chemistry module of EPACRIS (ExoPlanet Atmospheric Chemistry & Radiative Interaction Simulator), an innovative atmospheric simulation framework for exoplanets.The atmospheric radiative transfer module of EPARCRIS will be described in a separate paper (Scheucher & Hu, in prep).The EPACRIS atmospheric chemistry module integrates a cutting-edge automated chemical reaction network generation by RMG with a general-purpose one-dimensional photochemical kinetic-transport atmospheric simulation, originally developed by Hu et al. (2012), and since then expanded and upgraded substantially (Hu et al. 2013;Hu & Seager 2014;Hu 2019Hu , 2021)).This integration facilitates the fast and reliable construction of tailored reaction networks for specific exoplanet atmospheres.This paper details our methodology and demonstrates its effectiveness using the well-characterized atmosphere of WASP-39 b and the atmosphere of WASP-80 b as case studies for two different types of H 2 -dominated hot or warm Jupiters, compared with the recent JWST observations and photochemical modeling studies (Alderson et al. 2023;Ahrer et al. 2023;Feinstein et al. 2023;JWST Early Release Science Team et al. 2023;Rustamkulov et al. 2023;Tsai et al. 2023;Powell et al. 2024;Bell et al. 2023).

METHODS
A schematic diagram of the methodology and overall workflow adopted in this study is reported in Figure 1.

Automatic chemical reaction network generation using the Reaction Mechanism Generator (RMG) software
A detailed chemical reaction network for modeling the H 2 -dominated atmospheres of warm and hot Jupiters with equilibrium temperatures of 800-1500 K (T eq of WASP-39 b and WASP-80 b are within this range) was constructed automatically by RMG (Gao et al. 2016;Johnson et al. 2022).RMG is a Python-based opensource software and has been extensively used in the chemical engineering community to automatically generate chemical networks to simulate numerous pyrolysis and combustion chemistry successfully (Class et al. 2016;Dana et al. 2018;Chu et al. 2019;Keçeli et al. 2019;Liu et al. 2020).RMG is previously described in detail in Gao et al. (2016); Liu et al. (2021); Johnson et al. (2022) and only briefly described here along with Figure 1.
In given reactor conditions (i.e., temperature, pressure, reaction time, and initial mixing ratio of gas species), RMG will first place the initial species in the reaction system into the 'core' of the model and then find all the possible reactions based on these 'core' species (i.e., indicated as (1) on the left side of Figure 1).Chemical reaction rates depend on the species concentrations at the previous timestep.Thus, determining this initial set of core species is crucial when applying automatic reaction mechanism generation to (exo)planetary atmospheres.Unlike laboratory experiments, where initial concentrations are well-controlled, such information is often not fully available for exoplanet atmospheres.In response to this challenge, there are several ways to address this issue.
One approach to setting initial conditions is based on existing observational constraints.For instance, the recent JWST observation of K2-18 b has constrained concentrations of certain chemical species, such as methane (CH 4 ) and carbon dioxide (CO 2 ), from atmospheric retrievals (Madhusudhan et al. 2023).Building on this, researchers can construct the chemical network of interest by assuming specific compositions-e.g., approximately 50% H 2 O and 50% H 2 , along with constraints on CH 4 , CO 2 , and (CH 3 ) 2 S), which is tailored to explore the "Hycean Worlds" scenario of K2-18 b proposed by Madhusudhan et al. (2021).
Alternatively, one can define a set of grids, such as varying solar metallicity (Lodders 2020) from 1×to 100× solar metallicity as inputs for Reaction Mechanism Generator (RMG) and build a chemical network for different solar metallicities.Benchmarking these networks against existing JWST observations can check the feasibility of each chemical network.In the context of planetary atmospheres, constructing a chemical network tailored for a CO 2 -dominated atmosphere with a trace amount of sulfur species can provide insights into the atmospheres of Venus or Venus-like exoplanets.Such applications offer invaluable insights into (exo)planetary atmospheres.
After the reactor condition is provided, the next step for RMG (indicated as (2) in Figure 1) is to simulate the possible reactions using its database (maintained and updated with the latest data sources by the developers of RMG ( 2023)) of reaction parameters from previous experiments, ab initio calculations, or estimation methods (e.g., Benson group additivity, Benson & Buss (2004)), which will generate a list of possible product species  (i.e., 'edge' species).It should be noted that RMG relies on a chemical kinetics database compiled from various sources, each with its inherent errors.However, a well-maintained database represents our best knowledge at any given time.RMG initializes simulation at t=0 (indicated as (3) in Figure 1), followed by the next steps (indicated as (4-6) in Figure 1) that determine if these 'edge' species are important enough to be added to the 'core' species.'Edge' species i are included into the 'core' species if where R i is the production and loss flux of 'edge' species i, defined as an infinitesimal change in the concentration of 'edge' species (i.e., dC i ) in an infinitesimal time (i.e., dt), ϵ is the user-specified error tolerance, and R core is the characteristic flux of the reaction system, defined by A typically recommended range for this user-specified error tolerance, ϵ, is between 0.01 and 0.05 for users seeking a larger and more comprehensive model, despite the higher computational cost.Consequently, ϵ was set to 0.1 in this work (as specified by toleranceMoveToCore=0.1 in the RMG input file, available in the supplementary materials).As shown in (7-8) in Figure 1, the reaction generation and integration steps continue until they meet the termination criteria for reaction time, t term or the concentration of a specific species X term .This process results in the completed chemical network, encompassing all 'core' species and reactions with significant fluxes at the given reactor conditions.
In this work, temperatures from 700 to 2000 K and pressure between 10 -3 and 10 2 bar were sampled to generate chemical networks relevant within these T and P ranges using the ranged reactors setting in RMG (Liu et al. 2021), later then combined.We used an initial molecular mixing ratio of 10×solar metallicity, following the previous models of WASP-39 b (Tsai et al. 2023) and WASP-80 b (Bell et al. 2023), and automatically gener- The pressure dependence feature of RMG was enabled to automatically construct pressure-dependent networks for species with up to 10 atoms (i.e., constraining a total number of atoms).The resulting chemical network contained 105 species and 2337 reactions (forward-reverse reaction pairs), which can be found in the supplementary information as the CHEMKIN input file.Among these 2337 generated reactions, 2271 reactions didn't violate their respective collision limits, k coll (i.e., any bimolecular reaction rate coefficient doesn't exceed its Lennard-Jones collision rate constant), which were then incorporated into our 1D kinetic-transport model (Section 2.4) after the adaptation described in Sections 2.2 and 2.3.

Sorting out newly generated species and reactions by RMG
Besides the RMG-generated reactions and species outlined in Section 2.1, our initial kinetic-transport atmosphere model already possesses a reaction library that includes photodissociation and associated reactions.This library comprises 111 species (see Ta-ble 1) and 914 reactions, broken down as follows: 657 bi-molecular reactions, 91 ter-molecular reactions, 93 thermo-dissociation reactions, and 71 photochemistry reactions (Hu et al. 2012;Hu & Seager 2014;Hu 2021).The rates of photochemical reactions are calculated according to Equation ( 12), as detailed in Section 2.3 of Hu et al. (2012).Because RMG does not generate photochemistry-driven reactions, it is necessary to combine the original reaction network and the network generated by RMG, and the first step is to identify any overlapping species and reactions to prevent duplicates.We have annotated the 111 original species using RMG's 'adjacency lists' methodology (Gao et al. 2016;Johnson et al. 2022), which allows the RMG-EPACRIS adapter to compare the reactants and products, including reverse reactions, and ensure no duplications.The 'adjacency lists' method uses a graph-based structure to illustrate molecules, identifying atoms as vertices and their connecting bonds as edges in the list.For example, the adjacency list for hydroxymethylene (HCOH) and methoxy radical (CH 3 O) are shown in Figure 2. The structure of the 'adjacency lists' method is defined as follows: The first column specifies the atom index, the second column specifies the atom element, and the third column, prefixed by the lowercase letter 'u' for "unpaired", specifies the count of unpaired electrons for each atom.The fourth column, prefixed by the lowercase letter 'p' for "pairs", specifies the number of lone electron pairs.The fifth column, prefixed by the lowercase letter 'c' for "charge," specifies the formal charge on the atom.Bracketed values specify the presence of a bond, with the first value (i.e., number) indicating the index of the atom to which the current atom is bonded, and the second value (i.e., the uppercase letter) denoting the bond order: 'S' for single, 'D' for double, 'T' for triple, or 'B' for benzene-type bonds.If the molecule has an overall spin multiplicity (i.e., the degeneracy of the electronic ground state) larger than 1, it will be defined above the adjacency list (e.g., see the methoxy radical case in Figure 2).In the adjacency list of the methoxy radical molecule shown on the right side of Figure 2, the oxygen atom has a single unpaired electron (thus having an overall spin multiplicity of 2) and 1 single bond to the carbon atom that has 3 single bonds to hydrogen atoms, forming a methoxy radical.After sorting out newly generated species and reactions by RMG using the adjacency lists method, we found that 65 species (indicated with the footnote b in Table 1) were both included in the RMG-generated network and the original EPACRIS library, with 40 new species generated by RMG as shown in Table 2.We included nonreactive species, helium and Neon, and 20 additional

Yang & Hu
Note-a Simplified Molecular-Input Line-Entry System b 65 chemical species that were included in the original EPACRIS library and were also picked by RMG to be important for describing atmospheric conditions of WASP-39 b and WASP-80 b. c Despite appearing as other chemical species or different spin states in SMILES representation, these species are singlets in the 'adjacency lists' representation, indicating all electrons are paired.d 20 chemical species additionally included in the photochemical network to fully account for photochemistry and aerosol chemistry that might be important in the atmospheres of WASP-39 b and WASP-80 b.
chemical species (indicated with the footnote d in Table 1) as well as their other relevant thermochemical reactions imported from the original EPACRIS reaction library in the photochemical network.We found that the species that were not included in the RMG-generated list (i.e., thermochemically not important) are mostly associated with photodissociation (e.g., O( 1 D)) and the chemical network enabled by photodissociation.Consequently, we included them in our analysis to account for the impacts of photochemistry and aerosol chemistry, which might be significant in the atmospheres of WASP-39 b and WASP-80 b, except for the molecules that have more carbon atoms than C 3 hydrocarbons.Although C 3 and larger species are observed in Titan's atmosphere and mainly formed through photochemistry (Yung et al. 1984), we omitted them here due to the unfavorable physical (high temperature) and chemical conditions (H 2 -dominated) in hot Jupiters, and the uncertainties in relevant photodissociation rates.As a result, the final photochemical network contained 126 species and 2578 reactions (693 original reactions = 71 photo- Note-a Simplified Molecular-Input Line-Entry System Except for the 71 photochemistry reactions, the other 2484 reactions are forward-reverse reaction pairs.In the future, we will incorporate only the thermal-driven reaction rate coefficients that are generated by RMG, and only the photochemical reaction rate coefficients that are stored in the original reaction list into 1D kinetictransport modeling.

Converting CHEMKIN format to EPACRIS-readable format
After automatic chemical network generation by RMG as described in Section 2.1, the thermodynamic parameters (i.e., heat capacity, enthalpy, and entropy) and rate coefficients are provided in both CHEMKIN and Cantera (Goodwin et al. 2021) formats.The RMG-EPACRIS adapter then imports this information and converts it into a format adopted by EPACRIS to become inputs to 1D photochemical kinetic-transport atmospheric models, as described in the following Sections 2.3.1 and 2.3.2.

Thermodynamic parameters
RMG uses the NASA polynomial representation (McBride 1992) to calculate the relevant thermodynamic parameters.The NASA polynomial representation was originally developed by scientists at NASA to express temperature-dependent thermodynamic parameters such as the heat capacity C p (T ), enthalpy H (T ), and entropy S (T ) using seven or nine coefficients (McBride 1992).In this representation, the following thermodynamic parameters are given by nine polynomial coefficients a = [a −2 , a −1 , a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 ] (a −2 = a −1 = 0 in the seven-coefficient version): Then, the RMG-EPACRIS adapter obtains each species' Gibbs free energy, G, by the following equation: Yang & Hu The Gibbs free energy of species is used to calculate the reverse reaction rates in 1D photochemical kinetictransport models using the methods outlined in Hu & Seager (2014).

Rate-coefficient expressions
RMG adopts eight types of expressions for reaction rate constants.The RMG-EPACRIS adapter converts these into the formats adopted by EPACRIS (available in supplementary materials), enabling the importation of the rate constants, k, for 1D kinetic-transport atmospheric modeling.Consequently, EPACRIS implements the same eight rate-coefficient expressions, as elaborated below.
[1] Arrhenius-type expression Type 1 is the Arrhenius-type expression whose temperature-dependent rate-coefficient, k(T ), follows the Arrhenius equation as shown in Equation 7.
In the Arrhenius equation, A represents the preexponential factor, T 0 the reference temperature in kelvins [K], n the temperature exponent, and E 0 the activation energy in joules per mole [J • mol −1 ].Here, T denotes temperature [K], and R is the ideal gas constant, 8.314 The unit of A depends on the reaction order-[s −1 ] for first-order (i.e., unimolecular reaction, an elementary reaction in which the rearrangement of a single reactant produces one or more products), [m 3 • mol −1 • s −1 ] for second-order (i.e., bimolecular reaction, involving the simultaneous collision of any combination of two reactants), and [m 6 • mol −2 • s −1 ] for third-order (i.e., termolecular reaction, an elementary reaction involving the simultaneous collision of any combination of three reactants) reactions.
[2] Multi-Arrhenius-type expression Type 2 is the Multi-Arrhenius-type expression whose temperature-dependent rate-coefficient, k(T ), follows a set of Arrhenius equations summed to obtain the overall rate-coefficient, as shown in Equation 8.
In the Multi-Arrhenius equation, the index i refers to the ith set of Arrhenius parameters, which are consistent with those outlined in equation 7.
The rate coefficients are then determined by log scale interpolation between these expressions at each pressure.For example, the rate at an intermediate pressure In the Pdep-Arrhenius equation, P refers to pressure [P a].A(P ) represents the pressure-dependent preexponential factor, n(P ) the pressure-dependent temperature exponent, and E 0 (P ) the pressure-dependent activation energy in joules per mole [J • mol −1 ].The unit of A(P ) depends on the reaction order-[s −1 ] for first-order, [m 3 • mol −1 • s −1 ] for second-order, and [m 6 • mol −2 • s −1 ] for third-order reactions.The index i refers to the ith set of Arrhenius parameters.For pressures beyond the specified range, that is, P ≤ P lowest or P ≥ P highest , the rate-coefficient is determined by the Arrhenius equation at P lowest for lower, and at P highest for higher pressure values.
The rate coefficients are then determined by log scale interpolation between these expressions at each pressure as described in Equation ( 10).In the Multi-Pdep-Arrhenius equation, all terms are consistent with those defined in the Pdep-Arrhenius expression (Type 3).Unlike the Pdep-Arrhenius-type expression, which formulates the rate coefficient exclusively as either Type 1 or Type 2, the Multi-Pdep-Arrhenius-type expression offers a less consistent approach.In this equation, the rate coefficient expression can vary among Type 1 and Type 2, thus presenting a mixed format.
[5] Third-Body-type expression Type 5 is the Third-Body-type expression whose kinetics simply introduces an inert third body to the rate expression as shown in Equation ( 12).
In the Third-Body equation, k 0 (T ) refers to the lowpressure limit temperature-dependent rate coefficient.Its unit depends on the reaction order-[ for first-order, and [m 6 • mol −2 • s −1 ] for second-order reactions.
[M ] is the concentration of the bath gas [mol • m −3 ].
[6] Lindemann-type expression Type 6 is the Lindemann-type expression which qualitatively models the falloff behavior of pressuredependent reactions as shown in Equation (13).
k(T, P ) =k ∞ (T ) In the Lindemann equation, the Arrhenius expressions (i.e., Type 1) k 0 and k ∞ represent the low-pressure and high-pressure limit kinetics, respectively.The units of k 0 and k ∞ vary with reaction order: for first-order reactions, they are [m 3 • mol −1 • s −1 ] and [s −1 ], and respectively.
[7] Troe-type expression Type 7 is the Troe-type expression which quantitatively models the falloff behavior of pressure-dependent reactions by introducing a broadening factor F to the Lindemann equation, as shown in Equation 14.
The broadening factor F is computed via following Equation 15 In the Troe equation, the Arrhenius expressions (i.e., Type 1) k 0 and k ∞ represent the low-pressure and highpressure limit kinetics, respectively.The units of k 0 and k ∞ vary with reaction order: for first-order reactions, they are [m 3 • mol −1 • s −1 ] and [s −1 ], and for second- respectively.Four parameters (i.e., α, T 1 , T 2 , and T 3 ) are provided to calculate the broadening factor F .
[8] Chebyshev-type expression Type 8 is the Chebyshev-type expression which adopts the Chebyshev polynomial formulation as a means of fitting a wide range of complex k(T, P ) behavior as

Yang & Hu
shown in Equation 16.
logk(T, P ) = In the Chebyshev equation, α tp are the constants defining the rate coefficient, and ϕ n (x) is the Chebyshev polynomial of the first kind of degree n evaluated at x.The first few Chebyshev polynomials of the first kind are described in Equation 17.
T and P represent the reduced temperature and reduced pressures, respectively, mapping the ranges (T min , T max ) and (P min , P max ) to the interval (-1, 1).The Chebyshev rate expression is defined by the coefficient matrix α, comprising α tp , specific temperature and pressure ranges, typically involving 6 values for temperature (i.e., N T =6) and 4 for pressure (i.e., N P =4).It is important to note that Chebyshev polynomials are only defined within the interval (-1,1).Therefore, extrapolating rates beyond the defined temperature and pressure ranges is strongly discouraged, as the polynomials do not provide valid results outside these limits.

1-D photochemical kinetic-transport atmospheric modeling using EPACRIS
After generating the chemical network by RMG for the conditions relevant to the H 2 -dominated atmospheres of warm and hot Jupiters whose equilibrium temperature is 800-1500 K (Section 2.1) and adapting it for EPACRIS using the RMG-EPACRIS adapter (Sections 2.2 and 2.3), we performed 1D photochemical kinetic-transport atmospheric modeling with EPACRIS to simulate the steady-state mixing ratio of chemical species in the atmospheres of WASP-39 b and WASP-80 b.The photochemical kinetic-transport module of EPACRIS was employed to calculate the steady-state chemical composition of WASP-39 b's atmosphere of each morning and evening terminator (following Tsai et al. 2023) and that of WASP-80 b's atmosphere (following Bell et al. 2023), considering thermochemical equilibrium, vertical transport, and photochemical processes.We assumed cloud-free conditions and zero-flux boundary conditions.The temperature-pressure profiles (Figure 3a), the eddy diffusion coefficient profiles (Figure 3b), and the stellar spectra (Figure 3 c) are adopted from Tsai et al. (2023); Bell et al. (2023).In the case of WASP-80 b, we used the stellar flux of HD 85512 (K6V) at the 1 AU distance, adopted from the MUSCLES survey III (Loyd et al. 2016).It should be noted that the stellar spectrum significantly influences photolysis rates.From this, one can intuitively infer that the atmospheric chemistry of WASP-39 b is more significantly impacted by photochemistry compared to that of WASP-80 b based on Figure 3c, which shows that the stellar flux on WASP-39b is 10 to 100 times stronger than that on WASP-80 b.We assumed atmospheric abundances of 10×solar metallicity (Lodders 2020) for both WASP-39 b (Rustamkulov et al. 2023;Tsai et al. 2023) and WASP-80 b (Bell et al. 2023).These choices facilitate a direct comparison between the EPACRIS-simulated WASP-39 b and WASP-80 b atmospheres and published results.After the model has converged and reached the steady state, we computed the synthetic transmission spectra of WASP-39 b and WASP-80 b based on the molecular mixing ratio profiles (Figure 4 and 8), using the transmission spectra generation module of EPACRIS (Hu et al. 2013), and compared the resulting transmission spectra with JWST observations (Rustamkulov et al. 2023;Powell et al. 2024;Bell et al. 2023).Figure 4 shows the comparison between the previously reported vertical molecular mixing ratio profile of major species and those simulated from the current work using EPACRIS.As shown in Figure 4, the vertical mixing ratios of all species at pressures higher than ∼ 10 2 mbar are consistent across all five models, including EPACRIS and four others in Tsai et al. (2023).This consistency suggests that deep atmospheric chemistry of WASP-39b is primarily governed by thermal chemistry, which excludes photochemistry, and aligns with thermochemical equilibrium as depicted in Figure A1.The minor variations in the behaviors of SH and S 2 are primarily attributed to differences in the Gibbs free energy of these species.For example, VULCAN (Tsai et al. 2017(Tsai et al. , 2021)) 2000), and for S 2 from the NIST-JANAF table (Chase 1985).It is evident that each parameter in photochemical modeling inherently contains a certain amount of uncertainties.Among these, in general, rate coefficients and thermodynamic parameters are the most significant contributors to the uncertainties of chemical kinetic models.Therefore, assessing the model's sensitivity to thermodynamic parameters as well as rate coefficients is of paramount importance when it comes to the preciseness of atmospheric photochemical modeling.
At pressure less than ∼ 10 2 mbar (i.e., where photochemistry and vertical mixing become more significant), the EPACRIS simulations (solid lines in Figure 4) of the main sulfur species (H 2 S, SH, S 2 , S, SO, SO 2 ) align closely with the models in Tsai et al. (2023).The peak mixing ratios of these species are generally within an order of magnitude of each other for both morning and evening terminators, as illustrated in Figures 4a and  b.Although these models represent the atmospheric steady state, making it challenging to track the timedependent chemical evolution, some insights can still be gleaned.For instance, in the deep atmosphere (around ∼ 10 4 mbar), thermochemically favored H 2 S (hydrogen sulfide) is the dominant sulfur-bearing species up to a pressure of about 10 mbar (at temperatures above 900 K), for both terminators.In addition to H 2 S, sulfur monohydride (SH) is the second most abundant sulfur species in the deep atmosphere.Between 5-8 mbar, H 2 S rapidly transitions to S 2 and S, with SO and SO 2 also present.Above 1 mbar, atomic sulfur (S) becomes the dominant sulfur species, with SO and SO 2 peaking the maximum mixing ratio at ∼ 10 −1 − 10 −2 mbar for both the morning ([SO] max = ∼ 42 ppm, [SO 2 ] max = ∼68 ppm) and evening ([SO] max = ∼40 ppm, [SO 2 ] max = ∼28 ppm) terminators.One notable difference between the current model (i.e., EPACRIS) and the models previously reported in Tsai et al. (2023) lies in the amount of SO 2 levels predicted between 0.5 to 10 mbar for both the morning and evening terminators of WASP-39b (see Figure 4a and  b).In the morning terminator, Tsai et al. (2023) predicts SO 2 mixing ratios above 1 ppm (i.e., 10 −6 ) starting from around 10 mbar.In contrast, EPACRIS indicates this level is reached at around 2 mbar, a higher altitude (Figure 4a).Conversely, for the evening terminator, Tsai et al. (2023) suggests SO 2 exceeds 1 ppm starting from 0.5 mbar, whereas EPACRIS shows this occurring at around 3 mbar, a lower altitude (Figure 4b).This increased formation of SO 2 on the evening side compared to the morning side is explained in Section 3.1.2.It has to be noted that the mixing ratio prediction at around 1 mbar is potentially significant for JWST observations in transmission, especially considering the NIRSpec mode's primary probing range of 0.1 to 2 mbar (Rustamkulov et al. 2023).

Yang & Hu
Another notable difference is the amount of CH 4 levels predicted in the atmosphere above P ∼10 mbar for both the morning and evening terminators (see Figure 4a and b).According to the thermochemical equilibrium vertical mixing ratio profile of CH 4 as shown in Figure A1), CH 4 mixing ratio at around 10 mbar should be more than 100 ppm in the morning terminator, and more than 1 ppm in the evening terminator, but this is not the case since the quenching kinetics starts to happen around 100 mbar for both the morning and evening terminator, followed by photodissociation of CH 4 at higher altitude (P ∼1 mbar), as shown in Figure 4.However, on top of this, the current model shows a more depleted CH 4 level compared to previous modeling studies in Tsai et al. (2023), which indicates additional scavenging reactions for CH 4 sources (e.g., CH 3 ).According to the rate analysis, the reactions  abundance at around 1 mbar matches with the rapid appearance of S and S 2 originated from the dissociation of H 2 S (see Figure 4).Generally, sulfur dioxide (SO 2 ) at high altitudes (P ∼ 10 −2 − 0.5 mbar) is more prevalent at the cooler morning terminator of WASP-39b, whereas at lower altitudes (P ∼ 0.5 − 10 mbar), SO 2 is more abundant at the hotter evening terminator.To understand this distribution pattern, it is crucial to track the origin of oxidizers (i.e., OH and H radicals), since SO 2 in the atmosphere of WASP-39 b is mainly produced by successive oxidation of sulfur species originating from the deep atmosphere hydrogen sulfide (H 2 S) as pointed out in the previous modeling work of Tsai et al. (2023).
As illustrated in Figure 5, the vertical molecular mixing ratios of both OH and H radicals display similar patterns, largely due to the dissociation of H 2 O into H and OH.This also indicates that H 2 O is a main source for OH as well as H.However, their magnitudes differ, with additional sources of H such as H 2 and H 2 S contributing to these variations.The morning and evening vertical mixing ratio profiles of these species (i.e., OH and H radicals), as shown by the solid (morning) and dashed (evening) lines in Figure 5, cannot be fully explained by thermochemical equilibrium (dash-dotted line for the morning and dotted line for the evening) alone.This discrepancy indicates that a combination of thermal chemistry, photochemistry, and vertical mixing influences these behaviors.Further analysis, as shown in Figure 6a, indicates that the origins of OH radicals change with altitude (or pressure), suggesting a complex interplay of atmospheric processes at different levels.As shown in Figure 6a, from the upper atmosphere at around 10 −7 mbar down to 10 −3 mbar, the reactions serve as the major OH source and sink, respectively.
(rates ∼ ±2 × 10 5 [molecules/cm 3 /s]) Then in between 10 −3 − 10 −1 mbar, the reactions serve as the major OH source and sink, respectively.(rates ∼ ±3 × 10 7 [molecules/cm 3 /s]).And as shown in Figure 6a, at the pressure between 10 −1 − 10 1 mbar, the reactions HSO + H → OH + SH (23) serve as the major OH source and sink, respectively (rates ∼ ±6 × 10 9 [molecules/cm 3 /s]).The reaction rates get larger with decreasing altitude since molecular number density [molecules/cm 3 ] gets larger with decreasing altitude (i.e., increasing pressure ∝ number density).This interconversion of OH and H radicals is rapid, leading to the formation of a combined H + OH chemical group whose relative ratio remains constant under specific atmospheric conditions, as depicted in Figure 5.As highlighted previously, water vapor (H 2 O) is a primary source of these radicals.Rate analysis involving H 2 O (Figure 6 b) shows that the formation of OH and H, crucial for SO 2 production, primarily occurs through two distinct reactions in different atmospheric regimes.In the upper atmosphere, at pressures below 10 −2 mbar, H 2 O photolysis is the predominant reaction (see Figure 6a and b.Conversely, in the middle atmosphere, Yang & Hu within the pressure range of 1 to 10 mbar, the interaction between H 2 O and sulfur radicals (originated from H 2 S) becomes increasingly significant (as shown in Figure 6b).These interactions can be summarized by the following reactions: This summarized scheme is consistent with the rapid increase in H, OH around 10 mbar, as shown in Figure 5, and the rapid depletion of H 2 S along with a rapid increase in S at the same pressure level, as shown in Figure 4a and b.It is also noteworthy that the thermochemical equilibrium mixing ratio of H 2 O at pressures between 1-10 mbar is slightly higher than its vertical mixing ratio of H 2 O (see Figure 5).The transport rate (vertical mixing) in this region (i.e., P ∼1-10 mbar) was at least 2 orders of magnitude slower than the total loss rates.This suggests that vertical transport is not the predominant factor for the straight feature of H 2 O vertical mixing ratio profile at this region.Instead, this straight feature indicates that the decreased amount of H 2 O in this region is being converted into H and OH, corroborating the above reaction scheme.
If we combine the T − P profiles (see Figure 3a) with OH vertical mixing ratio of the morning and evening terminator (solid and dashed lines, respectively, in Figure 5), we can see the positive correlation between the OH (and H) radical mixing ratios and temperature within the 10 −1 − 10 1 mbar pressure range.Notably, in this range (i.e., 10 −1 − 10 1 mbar), the morning terminator consistently exhibits temperatures ∼200 K lower than the evening terminator.As described earlier, the origin of the lower atmospheric OH is more thermally driven chemistry while less H 2 O photolysis-driven, and thus sensitive to temperatures.As a result, within the 10 −1 − 10 1 mbar pressure range, OH radicals are more than an order of magnitude more abundant in the hotter evening terminator compared to the cooler morning terminator, as illustrated in Figure 5.This increased OH abundance leads to more SO 2 formation at lower altitudes (P ∼ 0.5 − 10 mbar) through OH-aided successive oxidation.
3.1.3.Theoretical transmission spectra of the atmosphere of WASP-39 b generated by EPACRIS Figure 7 shows the comparison between the averaged theoretical transmission spectra generated by EPACRIS and JWST observations of WASP-39 b (Alderson et al. 2023;Feinstein et al. 2023;Powell et al. 2024).The EPACRIS-predicted transmission spectra are broadly consistent with the near-and mid-infrared observations.The model accurately captures the features of H 2 O (blue), CO 2 (red), and SO 2 (magenta).While the CO (green) feature at 4.5-5 µm appears to be overpredicted, considering the uncertainties, the overall agreement between the model predictions and observations is still considered decent.The predicted transmission spectra generated by the previous photochemical networks, assuming 10× solar metallicity as discussed in Tsai et al. (2023), align well with the NIRSpec/G395H spectra (Feinstein et al. 2023).However, they overpredicted the transit depth in the 7-8 µm wavelength range, which corresponds to the dominant SO 2 absorption feature when compared to the MIRI data reported by Powell et al. (2024).Consequently, using the same photochemical networks resulted in the lowest χ 2 value of 2.51 when assuming 7.5× solar metallicity, while the χ 2 for 10× solar metallicity was 2.91 (Powell et al. 2024).In contrast, the predicted transmission spectra using the photochemical network from the current study showed an even lower χ 2 of 2.10.This network matches well with both NIRISS/SOSS, NIRSpec/G395H, and MIRI/LRS spectra without the need to vary solar metallicity.As highlighted in the previous Section 3.1.1,the vertical mixing ratios of species in the 0.1-2 mbar pressure range are probed by JWST observations in transmission.Differences in the photochemical network can lead to significantly varied results in the vertical mixing ratio profiles, thus influencing the predicted transmission spectra.

WASP-80 b
3.2.1.Upper atmospheric chemistry affected by the deep interior thermochemistry and quenching kinetics WASP-80 b's equilibrium temperature (T eq = 825K, Triaud et al. (2015)) is approximately 300 K cooler than WASP-39 b's (T eq = 1116 K, Faedi et al. (2011)).This suggests that transport-induced quenching, where the lifetime of chemical species becomes longer relative to the vertical mixing timescale, could play a more significant role in the cooler atmosphere of WASP-80 b compared to the hotter WASP-39 b. Figure 8 illustrates that all major species, including H 2 O, CO, CH 4 , NH 3 , and HCN, originate from the deep interior (the quenching point is at around P = 10 3 mbar) and are transported to the upper atmosphere, where some species such as CH 4 and NH 3 undergo photodissociation.Notably, the model-predicted H 2 O volume mixing ratio aligns well with the JWST observations (Bell et al. 2023)  as NH 3 , HCN, CO 2 , CO, and SO 2 were not constrained from the observations, all these model-predicted mixing ratios fall within the upper limit values derived from emission and transmission spectra except an SO 2 upper limit determined from emission spectra (see magenta left-pointing triangle symbol in Figure 8).As expected, the CH 4 mixing ratio as well as other species (e.g., NH 3 and HCN) formed from deep interior thermochemistry are sensitive to quenching kinetics.As depicted in Figure 8, the black solid and dashed lines represent the CH 4 volume mixing ratio using an eddy diffusion coefficient (K zz ) profiles that are 2 times slower and 5 times slower than the K zz profile adopted from Bell et al. (2023), respectively (see Figure 3b).When using these slower K zz profiles, the predicted CH 4 mixing ratio becomes more consistent with the observational constraints while not changing the predicted H 2 O mixing ratio, and NH 3 and HCN mixing ratios shift to the lower mixing ratios, indicating a shift of the deep interior quenching point toward lower pressure (see Figure 8).Consequently, more detailed constraints on other species such as CO 2 , NH 3 , CO, and SO 2 are required to precisely describe the WASP-80 b's atmospheric chemistry.Despite uncertainties in metallicity and the K zz profile, the current model aligns reasonably well with the observational data, providing valuable insights into the atmospheric behavior of this exoplanet.

Detailed chemistry and newly suggested deep-interior nitrogen incorporation pathway
Yang & Hu  (Powell et al. 2024) (grey square symbol points with uncertainties).The uncertainties are 1σ standard deviations.The reduced χ 2 value for the near-infrared region was calculated against the NIRISS/SOSS and NIRSpec/G395H data, while the reduced χ 2 value for the mid-infrared region was calculated against the MIRI/LRS data.Each color represents a spectrum generated by excluding specific species: blue for no H 2 O, red for no CO 2 , green for no CO, magenta for no SO 2 , and black for all species included.As shown in Figure 10a, the deep-interior CO-CH 4 conversion is (26) with M representing any third-body molecule.This scheme is identical to the scheme (1) in Moses et al. (2016).The SO 2 formation mechanism in the upper atmosphere predicted in the model was similar to the SO 2 formation in the atmosphere of WASP-39 b previously described in Section 3.1.2.The CO 2 formation mechanism was the combination of deep interior CO 2 formation and the additional oxidation of CO through OH radicals at the atmosphere above P ∼ 1 mbar.
Figure 10b visualizes the N 2 to NH 3 to HCN conversion pathway in the deep interior.This pathway mostly resembles those detailed in Moses et al. (2016), with a key difference in the initial incorporation of nitrogen from N 2 into species like NH 3 and HCN.The well-known N 2 →NH 3 route, as outlined in Moses et al. (2016), involves multiple hydrogenation steps starting with N 2 activation by a hydrogen atom to form NNH, ultimately yielding NH 3 as following: Despite this scheme being included in the RMGgenerated chemical network for hot Jupiter exoplanet atmospheres described in Section 2.1, RMG suggests a different dominant deep-interior pathway for N 2 →NH 3 conversion that initiates with N 2 directly interacting with two H 2 molecules to form N 2 H 2 , eventually leading to NH 3 (highlighted in Figure 10b) as following scheme: In this scheme, N 2 is activated by two H 2 molecules, forming cis-N 2 H 2 (see Figure 9).The transition state for this reaction was first identified by Asatryan et al. (2010) using the CBS-QB3 level of theory.In contrast to the simple bimolecular addition of N 2 +H 2 with an initial barrier of 125 kcal/mol (calculated at G2M(MP2)//MP2/6-31G** level of theory by Hwang & Mebel (2003)), this dihydrogen-activated nitrogen fixation has a relatively much lower barrier of about 77 kcal/mol (Asatryan et al. 2010).While the initial barrier in Scheme 27 is lower (17.1 kcal/mol, Hwang & Mebel (2003)), subsequent steps, such as NNH+H 2 →N 2 H 2 +H, face much higher barriers (42 kcal/mol, Hwang & Mebel (2003)), and Scheme 27 has to go through one more additional elementary reaction compared to Scheme 28.We tested the sensitivity of the photochemical model to the newly suggested reaction by intentionally disabling it.We found that under the conditions of WASP-80 b's atmosphere, the overall mixing ratio did not show significant changes.However, under certain favorable deep-atmospheric conditions (e.g., hot and high-pressure conditions), this newly suggested N 2 incorporation step could lead to significantly more formation of NH 3 or HCN species.
In this study, the rate coefficient for the N 2 +2 H 2 →N 2 H 2 +H 2 is considered as a high-pressure limit, calculated via conventional transition state theory.However, as this termolecular reaction involves three actual reactants (unlike the usual third body [M ] considered as an unreacted agent), it is inherently pressure-dependent and entropically less favorable.Despite these, such reactions are likely viable under the high pressures characteristic of the deep-interior chemistry in hot Jupiter atmospheres, underscoring the importance of its inclusion for accurate nitrogen incorporation modeling.This case highlights the substantial advantages of systematic, computer-aided automatic chemical network generation, which can reveal previously overlooked chemical pathways and provide detailed insights into exoplanet atmospheric chemistry.
As shown in Figure 10b, the NH 3 →HCN conversion scheme at pressures between 10 3 − 10 4 mbar is This pathway mostly resembles detailed in Moses et al. (2016) (Bell et al. 2023).The EPACRIS prediction aligns well with the NIRCam data (χ 2 =1.748 for Eureka! and 1.407 for tshirt), particularly in capturing H 2 O (blue) and CH 4 (green) features identified in Bell et al. (2023).Due to NIRCam's wavelength range limitations for CO 2 and SO 2 detection, comparing these species with model predictions is challenging.However, our model indicates that future JWST NIRSpec/G395H observations could potentially confirm the presence of CO 2 and SO 2 .Additionally, the spectral feature near 3 µm could signify NH 3 or HCN presence, both anticipated to exceed 1 ppm at pressures below 1 mbar (Figure 8).This observation underscores the need for more detailed exploration in this wavelength range, potentially through repeated observations.Overall, similar to the WASP-39 b case, there is a decent agreement between the model predictions and the observational data on the WASP-80b atmosphere.

Future applications and limitations of the current study
The current framework-automatic reaction mechanism generation coupled with one-dimensional (1-D) photochemical kinetic-transport modeling-has numerous potential applications for tools used in studying (exo)planetary atmospheres.One application is in the realm of climate modeling.Accounting for disequilibrium chemistry is essential in climate modeling.Consequently, reducing the size of the photochemical network becomes crucial, especially since two-dimensional (2-D) or three-dimensional (3-D) climate modeling is computationally intensive.This challenge amplifies in the context of general circulation modeling (GCM).
To address this, it is important to retain major chemical species that significantly impact climate structure while pruning less significant species from the chemical network to enhance computational efficiency.The current framework can offer substantial benefits to climate modeling and GCMs by eliminating unimportant NIRCam t-shirt (grey square symbol points with uncertainties) and Eureka!(grey circle symbol points with uncertainties) reductions by Bell et al. (2023).The uncertainties are 1σ standard deviations.The reduced χ 2 value of 1.407 was calculated against the NIRCam tshirt!reduction data.Although not shown, the reduced χ 2 value against the NIRCam Eureka! reduction data was 1.748.Each color represents a spectrum generated by excluding specific species: blue for no H 2 O, red for no CO 2 , green for no CH 4 , magenta for no SO 2 , and black for all species included.species and reactions.Pruning can be achieved by adjusting features in the Reaction Mechanism Generator (RMG), such as increasing the user-specified error tolerance, ϵ, or limiting the total number of atoms.However, this process involves a trade-off between minimizing the network size and maintaining the precision of the chemical network, necessitating a balance between these two aspects.
Moreover, several areas within the current framework could be improved.For instance, since RMG is primarily developed for simulating combustion chemistry, it lacks photochemical reactions in its library and does not account for vertical mixing processes (e.g., molecular diffusion and eddy diffusion).This omission might lead to significant gaps in identifying atmospherically important chemical species and reactions involving photons and vertical mixing.Therefore, incorporating photochemistry and physical processes in reaction mechanism generation is a critical area for future study.

CONCLUSIONS
In this study, we have developed a new framework for exoplanet atmospheric photochemical modeling.This framework, for the first time, integrates a rate-based automatic chemical network generator (RMG) with a 1-D photochemical kinetic-transport atmospheric simulator, forming the chemistry module of EPACRIS.We first constructed the reaction network specifically tailored for the atmosphere of H 2 -dominated hot Jupiter whose equilibrium temperature is 800-1500 K, and then incorporated this chemical network into EPACRIS for one-dimensional photochemical kinetic-transport modeling.Our model results generally align with previous studies of WASP-39 b by Tsai et al. (2023), particularly in capturing the photochemical production of OH radicals from H 2 O photolysis in the upper atmosphere and the formation of SO 2 through successive oxidation by OH and H radicals.A key difference between our study and previous models is the predicted SO 2 abundance in the middle atmosphere (pressure range of 0.5-10 mbar).Our results indicate a higher SO 2 formation at the warmer evening terminator compared to the cooler morning terminator.This discrepancy is attributed to an increased presence of sulfur species-oxidizers (OH and H radicals), predominantly generated from thermally-driven reactions between sulfur-bearing species (such as S or S 2 ) and H 2 O within this middle atmosphere.The predicted transmission spectrum of the atmosphere of WASP-39 b based on our model was compared to the JWST NIRISS/SOSS, NIRSpec/G395H, and MIRI/LRS observations (Feinstein et al. 2023;Alderson et al. 2023;Powell et al. 2024), showing good consistency and capturing H 2 O, CO 2 , and SO 2 spectral features for a 10×Solar metallicity atmosphere.Our model result of the WASP-80 b shows that the deep interior chemistry and vertical mixing dominate the general atmospheric chemistry, with predicted concentrations of CH 4 , H 2 O, CO, NH 3 , HCN, and SO 2 exceeding 1 ppm at pressures below 1 mbar.Utilizing RMG, we identified a dominant, previously overlooked reaction for the initial nitrogen incorporation (N 2 +2H 2 →N 2 H 2 +H 2 ), significant in high-pressure environments like deep interior atmospheres.Such use of RMG can unveil new reactions within chemical networks, potentially leading to the discovery of novel species in (exo)planetary atmospheres.The predicted transmission spectrum of the atmosphere of WASP-80 b based on our model was compared to the JWST NIRCam observation reported by Bell et al. (2023), showing good consistency and capturing H 2 O and CH 4 spectral features.This new approach not only provides the 1-D photochemical kinetic-transport modeling of (exo)planetary atmospheres with unprecedented efficiency and preciseness but also the applicability to di-verse atmospheric conditions (e.g., from H 2 -dominated to H 2 O-dominated atmospheres), enabling us to more effectively and precisely predict and interpret the vast amounts of data from upcoming JWST observations of various exoplanet atmospheric conditions.
Mention of any commercial product, process, or service by name, trademark, or manufacturer does not imply endorsement by the U.S. Government or the Jet Propulsion Laboratory, California Institute of Technology

Figure 1 .
Figure 1.A schematic diagram and the flowchart describing the expansion of the chemical network during automated reaction network generation by RMG using the rate-based algorithm (left) and overall workflow of implementation into 1-D photochemical transport atmospheric modeling (right) in this work.RMG stands for the Python-written "Reaction Mechanism Generator" (Gao et al. 2016; Liu et al. 2021; Johnson et al. 2022), while EPACRIS stands for the overall "ExoPlanet Atmospheric Chemistry & Radiative Interaction Simulator" written in C. Each light blue-colored shaded box refers to the corresponding output after each method described in Section 2.
Overall behavior of main sulfur-bearing species in the atmosphere of WASP-39b

Figure 3 .
Figure 3. (a) Temperature-pressure profiles for the morning (red) and evening limbs (blue) for WASP-39 b, and for WASP-80 b (black) (b) eddy diffusion coefficient (Kzz) profile for WASP-39 b (lime), WASP-80 b (black solid line), 2× slower eddy diffusion coefficient profile for WASP-80 b (black dashed line), and 5× slower eddy diffusion coefficient profile for WASP-80 b (black dotted line), and (c) the stellar flux at the 1 AU distance for WASP-39 (purple) and WASP-80 (green).These input parameters are adopted from Tsai et al. (2023); Bell et al. (2023) except for the stellar spectra of WASP-80.The stellar flux of HD 85512 (K6V) at the 1 AU distance, adopted from the MUSCLES survey III Loyd et al. (2016), was used for WASP-80 b.
, one of the models used in Tsai et al. (2023), utilizes thermodynamic parameters for SH and S 2 from Goos et al. (2016) and McBride (2002).These sources calculated NASA polynomials based on the data from McBride (2002).In contrast, RMG uses thermodynamic parameters for SH whose NASA polynomials are calculated by Song et al. (2017), which are based on the data from Shiell et al. (

Figure 4 .
Figure 4. Comparison between the previously reported vertical molecular mixing ratio profile of major species simulated for the (a) morning and (b) evening terminators of WASP-39 b in Tsai et al. (2023) (color-shaded areas) and those simulated from the current work using EPACRIS (solid lines).Each color indicates the corresponding species (SO 2 : magenta, H 2 O: blue, CH 4 : green, CO: red, CO 2 : dark blue, H 2 S: brown, S: purple, S 2 : grey, SH: yellow, and SO: light blue), and the color-shaded areas indicate the span enclosed by the photochemical models presented in Tsai et al. (2023).

Figure 5 .
Figure 5.The vertical molecular mixing ratio profile of several species (H 2 , H 2 O, H, OH) simulated for the morning (solid lines) and evening (dashed lines) terminators of the WASP-39 b atmosphere.The thermochemical equilibrium vertical molecular mixing ratios are also indicated by the dash-dotted lines (morning terminator) and dotted lines (evening terminator).Although not shown in the figure, the thermochemical equilibrium vertical molecular mixing ratios for the morning and evening terminators of OH show similar patterns to those of H radicals with the ∼6 orders of magnitude smaller amplitude.Each color indicates the corresponding species (OH: light blue, H: grey, H 2 O: blue, and H 2 : black).
3.1.2.Various origins of H and OH radicals for the formation of SO 2 in the atmosphere of WASP-39 b

Figure 6 .
Figure 6.The rate-pressure profile of dominant reactions involving specific chemical species in the morning terminator: (a) rates of reactions involving OH radicals presented at pressure ranges between 10 −7 -10 3 mbar.Among 318 OH-involved reactions, only 8 dominant reactions are shown here for readability.Each color of the solid lines indicates a corresponding reaction.Negative values in the rate indicate that the reaction acts as a sink for OH species, while positive values indicate the reaction serves as a source for OH species.The rate-pressure profile for the evening terminator exhibited behavior similar to that of the morning terminator, with the primary difference being in the amplitude of the rates (the evening rates are in general slightly faster within a factor of 2-3).Due to this similarity, the evening profile is not separately illustrated; (b) reactions involving H 2 O, presented at pressure ranges between 10 −2 -10 mbar.Among 124 H 2 O-involved reactions, only 6 dominant reactions are shown here for readability.Each solid line color corresponds to a specific reaction, with the colors representing the six rates in descending order from largest to smallest.Negative values in the rate indicate the reaction consumes H 2 O, while positive values indicate the reaction that produces H 2 O.As mentioned in the main text, in the upper atmosphere, H 2 O photolysis (lime) serves as a major source for OH species, while H 2 O + S→OH + SH reaction (blue) in the middle atmosphere at a pressure of 1-10 mbar serves as a major source for OH.

Figure 7 .
Figure 7.Comparison between the terminator-averaged theoretical transmission spectra generated by EPACRIS (solid lines) and the JWST observations: (top) NIRISS/SOSS and NIRSpec/G395H data (Feinstein et al. 2023; Alderson et al. 2023) (grey circle symbol points with error-bars indicate NIRISS/SOSS data, while grey square symbol points with uncertainties indicate NIRSpec G395h) and (bottom) MIRI/LRS data(Powell et al. 2024) (grey square symbol points with uncertainties).The uncertainties are 1σ standard deviations.The reduced χ 2 value for the near-infrared region was calculated against the NIRISS/SOSS and NIRSpec/G395H data, while the reduced χ 2 value for the mid-infrared region was calculated against the MIRI/LRS data.Each color represents a spectrum generated by excluding specific species: blue for no H 2 O, red for no CO 2 , green for no CO, magenta for no SO 2 , and black for all species included.

Figure 8 .
Figure 8.Comparison of the vertical molecular mixing ratio profiles for WASP-80 b's atmosphere: previous JWST observations (Bell et al. (2023), symbols and dashed lines) versus current EPACRIS simulations (solid lines for 1D photochemical kinetic-transport modeling; dotted lines for thermochemical equilibrium).Square symbols with uncertainties indicate the emission data, and circle symbols with uncertainties indicate the transmission spectra.Left-pointing triangle symbols indicate the upper limit (U.L.) values determined from emission spectra, while down-pointing triangle symbols indicate the upper limit values determined from transmission spectra.Each color indicates the corresponding species (SO 2 : magenta, H 2 O: blue, CH 4 : green, CO: red, CO 2 : dark blue, NH 3 : yellow, and HCN: cyan).The dark solid and dashed lines are the predicted CH 4 volume mixing ratio when using 2× slower Kzz profile and 5× slower Kzz profile in Figure 3b, respectively.Although not shown, the H 2 O volume mixing ratio using 2× slower Kzz profile is almost identical to the H 2 O volume mixing ratio using the original Kzz profile of WASP-80 b.

Figure 9 .
Figure 9.A schematic diagram that visualizes the first reaction:N 2 +H 2 +H 2 →N 2 H 2 +H 2 in Equation28.From the left side, N 2 reacts with two H 2 molecules to form cis-N 2 H 2 and H 2 on the right side.

Figure 11
Figure 11  compares EPACRIS-generated theoretical transmission spectra with the JWST NIRCam observations of WASP-80 b(Bell et al. 2023).The EPACRIS prediction aligns well with the NIRCam data (χ 2 =1.748 for Eureka! and 1.407 for tshirt), particularly in capturing H 2 O (blue) and CH 4 (green) features identified inBell et al. (2023).Due to NIRCam's wavelength range limitations for CO 2 and SO 2 detection, comparing these species with model predictions is challenging.However, our model indicates that future JWST NIRSpec/G395H observations could potentially confirm the presence of CO 2 and SO 2 .Additionally, the spectral feature near 3 µm could signify NH 3 or HCN presence, both anticipated to exceed 1 ppm at pressures below 1 mbar (Figure8).This observation underscores the need for more detailed exploration in this wavelength range, potentially

Figure 10 .
Figure 10.Model-predicted major reaction pathways at deepinterior (P ∼ 3 × 10 4 mbar) for (a) CO-CH 4 conversion, and (b) N 2 -NH 3 -HCN conversion.Bold characters represent major species with a mixing ratio above 1 ppm.Dashed pathways indicate directions at least three times slower than other branching reactions (e.g., N 2 H 2 branching into N 2 H 3 is 3× slower than into N 2 H 4 ).The yellow highlighted region in panel (b) indicates a newly suggested initial nitrogen incorporation path in the RMGgenerated chemical network.

Figure 11 .
Figure 11.Comparison between the theoretical transmission spectra generated by EPACRIS (solid lines) and the JWST observations:

Table 1 .
111 molecular species originally included in the EPACRIS library

Table 2 .
40 newly included molecular species in the chemical network tailored for H 2 -dominated atmospheres by RMG , with additional species (i.e., CH 3 NH and CHNH are newly included species by RMG) and reactions (i.e., CH 3 NH 2 +H→CH 3 NH+H 2 , CH 3 NH