© The Authors 2026

1 Introduction

Over the past several decades, and continuing today, the chemical characterization of exoplanet atmospheres has relied on a sustained community effort. This difficult and delicate task has generated many research programs on the acquisition and treatment of observational data for these celestial bodies, but also on the modeling of these data.

The commissioning of the James Webb Space Telescope (JWST) in July 2022 has flooded the exoplanet community with unprecedented, high-quality data, enabling rapid progress in our understanding of these worlds.

Forthcoming observational programs, including space-based missions such as ARIEL, scheduled for launch in 2029 (Tinetti et al. 2018), and ground-based instruments such as the Extremely Large Telescope, scheduled for first light in 2028, will intensify the need for precise modeling to ensure that these data are used to their full potential. The present study is motivated by recent results acquired by the JWST, which reveal the presence of SO₂ in WASP-39 b and WASP-107 b (Tsai et al. 2023; Dyrek et al. 2024; Powell et al. 2024; Welbanks et al. 2024).

The implications of these detections are twofold. Firstly, SO₂ is the first sulfur-bearing species ever detected in an exoplanet atmosphere, which confirms the importance of sulfur chemistry for exoplanets, as it demonstrates that these species can affect the spectrum. This opens the door to its quantification and potentially to retrieving information on the planet’s formation through the ratios with other elements (Cridland et al. 2019).

Secondly, SO₂ is thermodynamically unstable in these atmospheres, as it is a highly oxidized form of sulfur. In the hydrogen-dominated atmospheres of hot Jupiters such as WASP-39 b and WASP-107 b, the thermodynamically stable form of sulfur is H₂S, which, unlike SO₂, has not been unambiguously detected. This indicates that photochemical processes oxidize H₂S to SO₂ (Tsai et al. 2023), making it necessary to include these species in models to estimate their chemical abundances and understand their impact on synthetic spectra. This is generally neglected in most retrieval analyses, which typically assume thermochemical equilibrium or constant abundance profiles (Wakeford et al. 2017; Kreidberg et al. 2014). This is because this modeling is computationally intensive as it requires integrating the chemical kinetics occurring inside the atmosphere until the steady state is reached. These chemical processes are represented as a kinetic network, the size of which increases rapidly with its comprehensiveness and can become prohibitive for 3D modeling or retrievals. However, for these applications, this network size can be reduced using different reduction approaches (Lira-Barria et al. 2022; Venot et al. 2019; Tsai et al. 2022; Yang & Hu 2024) and used to perform retrievals with TauREx thanks to its new FRECKLL plugin (Al-Refaie et al. 2024). Another constraint on the modeling of chemical kinetics is its reliance on two types of crucial data.

The first type is accurate photolysis data, which are available for most major species at room temperature, but can deviate drastically under the high-temperature conditions of exoplanet atmospheres (Venot et al. 2013; Venot et al. 2018; Fleury et al. 2025). The second is kinetic data for the reactions in the network, whose rates are parameterized with rate constants that depend heavily on temperature and also pressure for some reaction types. While data for earth-like conditions (300 K, N₂–O₂ atmosphere) are available for simple and abundant species, higher-temperature data for hydrogen-dominated atmospheres are more difficult to gather. Several approaches have been adopted to address this problem. The most common approach is to gather multiple rate constants from various databases, such as KIDA (Wakelam et al. 2012) or the NIST chemical kinetics database (Manion et al. 2008). On the one hand, this approach can lead to small network sizes, which lower computational cost by enabling the inclusion of only relevant species. On the other hand, it is biased toward including only reactions available in the literature that directly involve the main species. This can lead to internal inconsistencies in the network and the neglect of crucial species and reactions to accurately model the kinetics of the species of interest.

To solve this issue, networks directly derived from those used in combustion modeling were developed, as these usually undergo extensive validations¹ against experiments in ideal reactors (Venot et al. 2012; Venot et al. 2015; Venot et al. 2020). Our previous work also focused on deriving such a network from the most recent advances in combustion and pyrolysis kinetics for C/H/O/N compounds and separately performing validation against a wide range of combustion experiments (Veillet et al. 2024).

We take the same approach in this work, starting from available sulfur networks in the combustion and pyrolysis literature that are used to model the thermal decomposition of hetero-atomic organic compounds in many combustion applications. These compounds may contain oxygen, sulfur, nitrogen, phosphorus, or halogens, and are the building blocks of much combustion-relevant matter, such as materials (Xu & Hadjichristidis 2023), pesticides (Casida & Durkin 2017), toxic agents (Sirjean et al. 2017), pharmaceuticals, or fuels, either derived from organic matter or fossil resources, such as biogas or coal-derived syngas (Mathieu et al. 2015). The combustion of these hetero-atomic compounds is studied both for applications under controlled conditions, such as the incineration of waste and its energy recovery, but also to prevent accidental combustion such as fires and explosions. Forest fires also fall into this category, in which biomass containing hetero-atoms is burned. These phenomena involve both solids and gases, and, in most cases, the first step required for combustion to take place is an initial pyrolysis stage (Stauffer et al. 2007). For the combustion of gases with impurities, the literature shows that the presence of hetero-atoms such as sulfur or nitrogen perturbs the reactivity of organic compounds and the type and quantity of pollutants emitted (Mathieu et al. 2014; Honorien et al. 2023; Glarborg et al. 2018). In combustion processes, the rapid decomposition of fuel is followed by fragment oxidation, which is kinetically limiting in the fuel oxidation process (Wang et al. 2018). Robust reaction bases, which have been validated over a wide range of pressure and temperature conditions, are therefore also required in combustion research on hetero-atomic organic compounds, where the temperatures encountered are similar to those found in hot exoplanets. Although these core reaction mechanisms are widely available for hydrocarbon combustion, they remain scarce for describing the combustion chemistry of hetero-atomic compounds containing both sulfur and nitrogen. Nitrogen chemistry sub-mechanisms were originally included along with the combustion chemistry of small hydrocarbons for NO_x formation. Recently, some of the nitrogen sub-mechanisms also account for NO_x reduction reactions with ammonia (Glarborg et al. 2018). Kinetic models incorporating the combustion chemistry of sulfur compounds were developed mainly for H₂S, to account for its presence as an impurity in some types of gas and, to a lesser extent, for CH₃SH (Stagni et al. 2022; Mathieu et al. 2014; Zhou et al. 2013; Cooper et al. 2022). In this work, we used experimental and theoretical combustion and pyrolysis data, along with our previous work on a C/H/O/N chemical network (Veillet et al. 2024) and CH₃SH kinetics (Veillet et al. in rev), to develop our network for sulfur modeling in exoplanets, with fully coupled C/H/O/N/S interactions, especially between carbon and sulfur. This C/H/O/N kinetic network has been statistically validated against a database of 1618 experimental combustion and pyrolysis datasets (Veillet et al. 2024), featuring CH₄, CH₃OH, C₂H₅OH, CO, H₂, H₂O, NO, HCN, NH₃ and their mixtures. We used a similar approach to select a H₂S kinetic model from the literature, based on extensive comparisons with a newly developed sulfur combustion database (SO_x, H₂S, and CS₂) composed of 1151 experimental measurements. For methanethiol (CH₃SH), the main link in the carbon-sulfur chemistry, we propose a new detailed pyrolysis and combustion mechanism based on previous theoretical calculations and validate it against recent experimental data (Colom-Díaz et al. 2021; Alzueta et al. 2019). We then used this kinetic network to model both hot Jupiters and warm Neptunes, and to investigate the effect of sulfur chemistry on other species and their synthetic transit spectra.

Section 2 discusses the choices and available data used to build the C/H/O/N/S kinetic network, both theoretical and experimental. In Sect. 3, this network is used to model the atmospheres of six exoplanets: GJ 436 b, GJ 1214 b, HD 189733 b, HD 209458 b, WASP-39 b, and WASP-107 b. We also perform a detailed chemical analysis on the main species of interest to identify the main contributing reactions and formation pathways for these species, and to investigate how the coupling between sulfur, carbon, and nitrogen chemistry impacts their abundances. The corresponding synthetic spectra are then simulated to evaluate the impacts on observables. Finally, we conclude in Sect. 4 and discuss potential future improvements to this work.

2 Detailed combustion network selection

2.1 Considered combustion networks

The starting point of this work is a recently developed C₀–C₂ C/H/O/N core mechanism, built based on a systematic comparison of the performance of several available mechanisms from the literature for C/H/O and C/H/O/N chemistry against experimental data (Veillet et al. 2024). The kinetic model includes oxygenated or nitrogen compounds such as methanol, ethanol, ammonia and hydrogen cyanide, and small hydrocarbons (C₂H₆, C₂H₄, C₂H₂, and CH₄). We gathered experimental data on combustion and pyrolysis of these species from the literature, covering a wide range of temperatures (800–2400 K) and pressures (0.2–50 bar). The experimental measurements included in this dataset consist of ignition delay times and species profiles as a function of temperature, time, and pressure in shock tubes, plug flow, jet-stirred, and closed reactors. To develop the present C/H/O/N/S mechanism, we built on this previous work and added sulfur chemistry following the same methodology. We gathered experimental measurements of the combustion and pyrolysis of four different species: H₂S, CH₃SH, OCS, and CS₂. This dataset consisted of species profiles as a function of temperature in flow reactors and jet-stirred reactors. The resulting dataset, composed of 1606 experimental measurements, is summarized in Table 1. We then compared four kinetic models from the literature on this dataset: Stagni et al. (2022), Cooper et al. (2022), Zeng et al. (2021), and Colom-Díaz et al. (2021), which we briefly describe below.

Stagni 2022: a mechanism for H₂S combustion and pyrolysis from Stagni et al. (2022), built and validated on a large review of H₂S experimental data. It is one of the state-of-the-art model for H₂S combustion and pyrolysis kinetics. No carbon chemistry is included.

Colom-Diaz 2021: a mechanism for CH₃SH combustion and pyrolysis from Colom-Díaz et al. (2021), which iterates on previous work on atmospheric CH₃SH kinetics.

Cooper 2022: a state-of-the art mechanism for H₂S, built on the work of Haynes and co-workers and Glarborg and coworkers, which includes its oxidation by N₂O from Cooper et al. (2022). This mechanism is one of the few models that attempt to describe N/S coupling.

Zeng 2019: a mechanism for CS₂ combustion from Zeng et al. (2019), based on the CS₂ kinetics from a previous mechanism proposed in Glarborg et al. (2014).

2.2 Comparison with experiments

To test the accuracy of these models on H/O/S chemistry, we first compared them with our H₂S combustion and pyrolysis data using the Cantera software (Goodwin et al. 2023). Figs. 1 and 2 show two examples of these comparisons. In these comparisons, the model from Stagni et al. (2022) shows the best agreement with the experimental data. This is especially true for H₂S pyrolysis, since the other three models show an offset of approximately 100 K in their H₂S consumption and H₂ production. For the combustion data, the Cooper 2022 model also performs well, with a temperature shift of around 50 K in H₂S consumption and SO₂ production. By contrast, the Zeng 2019 and Colom-Diaz 2021 models are in poorer agreement with the experimental data. The Zeng 2019 model shows a temperature shift of around 150 K, while the Colom–Diaz model exhibits strong variations (likely due to numerical instabilities) in H₂S consumption and SO₂ production between 700 and 850 K. We observe these variations in all simulations of the H₂S combustion dataset for this model with different convergence parameters; therefore, they are not due to convergence issues.

Consequently, to implement H₂S chemistry in the present work, we merged the Stagni 2022 model with our previous C/H/O/N mechanism (Veillet et al. 2024). As Stagni 2022 already contains an H/O reaction subset, we replaced it with the one from our previous study and verified that this did not affect its performances on our H₂S experimental dataset. Because the Stagni 2022 model only described H/O/S chemistry, C/S coupling had to be added to create our final mechanism. To fully cover the C/H/O/S chemistry, we first included the Colom-Diaz 2021 and Zeng 2019 models in our mechanism, excluding their H₂S chemistry subset, which was already included from Stagni 2022. As the Colom-Diaz 2021 model mainly described CH₃SH combustion and pyrolysis, it also contains a CS₂ reaction subset, which we replaced with the one from Zeng 2019, which specializes in CS₂ and OCS combustion. However, as shown in Fig. 1, the Colom-Diaz 2021 model exhibits numerical issues, which led us to develop a new kinetic model for CH₃SH. Most parts of this new sub-mechanism were estimated by analogies with close analog oxygenated species sharing the same molecular geometry (i.e., CH₃OH for CH₃SH, CH₂O for CH₂S). Activation energies were corrected for differences in reaction enthalpy when necessary, and a total of more than 120 reactions were included in this CH₃SH. The rate coefficients of the most sensitive reactions for pyrolysis simulations and unknown thermochemistry data were calculated using ab initio calculations at the CCSD(T)-F12/CBS level. Finally, to account for the N/S coupling, we added the corresponding reaction subset from Cooper 2022 to the present work, which included 40 additional reactions. This work is fully detailed in Veillet et al. (in rev) and is not described in more detail here.

Fig. 1 Comparison between models from the literature for the combustion of 500 ppm of H2S at atmospheric pressure and an air-to-fuel equivalence ratio of four in a flow reactor (no. 1 in Table 1).

2.3 Thermal decomposition of H₂S

The current understanding of the thermal decomposition of H₂S remains debated and has required us to make important choices that we detail here. The literature review of Raj et al. (2020) on sulfur combustion, particularly H₂S kinetics, further details the conflicts between modeling and experiment, also noted in Stagni et al. (2022). In essence, the two reactions in question are the possible pathways for thermal decomposition of H₂S: H₂S → SH + H and H₂S → S + H₂. The first pathway was the most implemented in early kinetic networks, including H₂S pyrolysis (Higashihara et al. 1976; Bowman & Dodge 1977; Roth et al. 1982), as it is a simple bond-breaking reaction and a standard, well-known type of initiation reaction. The second pathway is a spin-forbidden reaction, asH₂S and S are triplets and H₂ a singlet. Consequently, this reaction must proceed through an intersystem crossing from the singlet potential energy surface to the triplet energy surface. Because this process is typically slow and difficult to model theoretically, it was initially ignored in models until experimental evidence for this decomposition channel was reported by Woiki & Roth (1994) and Olschewski et al. (1994). A subsequent theoretical study by Shiina et al. (1996) analyzed these potential energy surfaces and the crossing point energy, but the derived rate constant was an order of magnitude lower than those measured by Woiki & Roth (1994) and Olschewski et al. (1994). Furthermore, additional experimental measurements by Shiina et al. (1998) and Karan et al. (1999) confirmed the previous experimental values, indicating that this theoretical value was likely underestimated. Shiina et al. (1998) also experimentally measured the corresponding reverse reaction, H₂ + S → H₂S and compared it to the channel H₂ + S → SH + H, showing that the latter is favored at 100 bar at 900 K. However, these experimental measurements are also subject to caution, as Zhou (2009) demonstrated that even quartz surfaces can catalyze H₂S decomposition. Nonetheless, the activation energy deduced for H₂S → S + H₂ from these independent experiments is in direct conflict with ab initio calculations and thermochemical data, because the corresponding activation energy is lower than the reaction enthalpy. This results in a negative activation energy for the reverse reaction H₂ + S → H₂S, which favors this pathway over H₂ + S → SH + H reaction, in direct conflict with the measurements of Shiina et al. (1998). Such an activation energy is unphysical, but it is the only one that reproduces both high- and low-temperature experimental data. To reconcile this activation energy with the correct reverse pathway, we used an intermediate species corresponding to the first excited state of the sulfur atom S(¹D). This species was included in our mechanism by splitting the reaction H₂S → S + H₂ into two reactions: H₂S → S(¹D) + H₂ and S(¹D) S. The former was written with the same parameters as H₂S → S + H₂, while the latter used the same parameters as the analog reaction O(¹D) → O(³P). For S(¹D) thermochemical data, we used the formation enthalpy of the S(³P) ground state (denoted simply as S in this work) and added the energy difference with S(¹D). For the entropy, we used the Sackur–Tetrode equation (Sackur 1913; Tetrode 1912). Because the S(¹D) → S reaction is much faster than H₂S → S(¹D) + H₂, the total reaction rate for the forward reaction remains unchanged. The reverse reaction is prevented by the step S → S(¹D), which is very slow due to the high energy barrier. This approach effectively reconciles the observed activation energy with the preferred reverse pathway, although it does not resolve the issue of the non-physical activation energy for H₂S → S(¹D) + H₂. Further theoretical and experimental studies on H₂S pyrolysis are required to fully resolve this conflict, but this solution provides the best compromise given the current state of the literature.

Fig. 2 Comparison between models from the literature for H2S pyrolysis (500 ppm at atmospheric pressure in a flow reactor, no. 2 in Table 1).

3 Application to exoplanetary atmospheres

3.1 Models and data sources

Throughout the development of this experimentally validated C/H/O/N/S kinetic network, we focused on the coupling of sulfur chemistry with the rest of the network. We primarily aimed to model the impact of sulfur chemistry in hot exoplanets, both on sulfur species and on the main C/H/O/N species (such as H₂O, CH₄, CO, CO₂, NH₃, HCN, C₂H₂, andC₂H₄). In this section, we use this kinetic network to model the atmosphere of six exoplanets: three warm Neptunes (GJ 436 b, GJ 1214 b, WASP-107b) and three hot Jupiters (HD 189733 b, HD 209458 b, WASP-39 b).

WASP-39 b and WASP-107 b were chosen specifically because of recent observations of SO₂ in their transit spectra by the JWST (Tsai et al. 2023; Dyrek et al. 2024), while HD 189733 b, HD 209458 b, GJ 436 b, and GJ 1214 b are planets that have been extensively studied and are commonly used to represent their respective exoplanet classes (Venot et al. 2012, 2020; Moses et al. 2011). The parameters used to model these planets (pressure-temperature profile, UV stellar flux, eddy diffusion coefficient, and metallicity) are the same as in Veillet et al. (2024) for GJ 436 b, GJ 1214 b, HD 189733 b, and HD 209458 b, while for WASP-39 b and WASP-107 b we used the same input parameters as Tsai et al. (2023) and Dyrek et al. (2024).respectively. Table 2 summarizes these parameters. We used the 1D model FRECKLL (Al-Refaie et al. 2024) to run the simulations with both our newly developed kinetic network and the sulfur kinetic network from VULCAN used in Tsai et al. (2023) to simulate WASP-39 b, hereafter referred to as V25 and Tsai 2023. The pressure-temperature profile was discretized in 130 layers, and the photodissociation data used for the sulfur species are listed in Table A.1. Fig. 3 shows the obtained abundance profiles. To study the effect of the coupling on the C/H/O/N species, we also compared our results with those obtained using the scheme developed in Veillet et al. (2024) (without sulfur but featuring the same C/H/O/N chemistry), hereafter referred to as V23, is shown in Fig. 4. The corresponding synthetic transmission spectra at a resolution of 50 were also generated using TauREx (Al-Refaie et al. 2021). We used opacity data from Exomol (Tennyson & Yurchenko 2012) for H₂O, CH₄, CO, CO₂, NH₃, HCN, C₂H₂, C₂H₄, H₂CO, OH, CH₃, CN, NS, SH, SO, SO₂, SO₃, H₂S, CS, and OCS, and HITRAN (Gordon et al. 2022), downloaded from the petitRADTRANS website (Mollière et al. 2019), for CS₂. Collision-induced absorption was included for H₂ H₂ and H₂ He, using data from HITRAN (Gordon et al. 2022). The calculations also included the default Rayleigh scattering from TauREx. The resulting spectra are shown in Fig. 5.

In the following, we discuss the main differences observed in the synthetic spectra. Detailed contributions are shown in Fig. B.1. We investigated the causes of these differences using the abundance profiles. We computed the formation pathways responsible for changes in the abundance profiles by solving the minimum flow problem from graph theory using the Python package NetworkX. Nodes were chosen as species in the network while edges were determined from reactions, with a capacity equal to their rate. We determined the nodes linked by an edge for each reaction by maximum matching on a bipartite graph formed from the set of reactants and the set of products. We based the maximum match chosen for each reaction on a greedy algorithm (Chvatal 1979) using a cost matrix calculated from the number of identical and different atoms for each match. This choice favors links between the product that most closely resembles one of the reactants, with a bias toward reactants that share the highest number of identical atoms. This is optimal for simple abstraction bimolecular reactions, for example, but is limited for more complex cases such as association reactions. For this reason, when identifying pathways between species, we allow the user to specify which atoms must be included in all intermediate species. Therefore, when generating the graph to calculate the total flow, we only include species that contain these atoms and exclude those that did not fulfill this condition. This results in a directional graph that represents the chemical pathways available from one species to another, together with the corresponding flows. This approach enables quantitative characterization of both primary and secondary formation pathways. In the following, a chemical pathway (with its corresponding flow) refers to the list of intermediate species from the source to the sink. Each link between intermediate species is composed of the sum of the rates of all the reactions that link these species. Two important properties of these “steps” in the chemical pathway are therefore the contribution of each reaction to the total capacity of the link and the ratio of this capacity to the global flow between the source and sink, which we call saturation in the following. This saturation provides information on the sensitivity of the total flow of the chemical pathway to the reaction rates and is by definition 100% for the limiting step in the pathway. The global flow of this pathway will therefore be less sensitive to variations in the rates of reactions with low saturation.

Fig. 3 Abundance profiles of GJ 436 b, GJ 1214 b, HD 189733 b, HD 209458 b, WASP-39 b, and WASP-107 b simulated with FRECKLL. Solid lines show results from the present work, while dashed lines show those from Tsai 2023. We omit H2 and He for clarity to focus on the main species.

Fig. 4 Abundance profiles of GJ 436 b, GJ 1214 b, HD 189733 b, HD 209458 b, WASP-39 b, and WASP-107 b simulated with FRECKLL. Solid lines show results from the present work, while dashed lines show results from V23.

3.2 Results for GJ 436 b and GJ 1214 b

For GJ 436 b, we modeled the planet using a 10× solar metallicity (from Lodders (2010)), as exact values are uncertain but likely at least one order of magnitude above solar (Madhusudhan & Seager 2011; Moses et al. 2013). Recent data from the JWST (Mukherjee et al. 2025b) suggest that its metallicity could be closer to 100× solar, indicating that this value might be underestimated. In general, some planet parameters, such as metallicity, may deviate from current estimates or be revised by future data. When no exact literature value is available, we used general estimates that represent one of the possible cases. For GJ 1214 b, we chose a higher metallicity of 100× solar, as suggested by recent observations of JWST (Nixon et al. 2024; Kempton et al. 2023). For GJ 436 b, a constant eddy diffusion coefficient of 10⁸ cm²/s was used in the absence of a precise profile from a General Circulation Model (GCM) calculations, while for GJ 1214 b, we adopted the same profile as in Veillet et al. (2024). The resulting abundance profiles, shown in Fig. 3, reveal several differences between our sulfur network (solid lines) and Tsai 2023 (dashed lines). Firstly, with our network, the quenching point of CO is slightly higher in the atmosphere of GJ 436 b, resulting in increased CO and CO₂ above 10 bar. The sole pathway converting CO to CO₂ is the single elementary reaction CO + OH → CO₂ + H. Secondly, the abundance of CS₂ is clearly higher with our network, especially between 1 and 10⁻⁴ bar, with a difference of two to three orders of magnitude for GJ 436 b and one to two orders of magnitude for GJ 1214 b. This discrepancy arises from differences in the CS₂ formation mechanisms between the two networks. For V25, CS₂ forms from CH₄ through CH₃, CH₂S, CHS, and CS via the following steps: $$C{H}_4+\text{H}\to C{H}_3+{\text{H}}_2$$(1) $$C{H}_3+\text{S}\to C{H}_2\text{S}+\text{H}$$(2) $${\text{CH}}_2\text{S}+\text{H}\to \text{CHS}+{\text{H}}_2$$(3) $$\text{CHS}+{\text{S}}_2\to \text{CS}+{\text{HS}}_2$$(4) $$\text{CS}+\text{SH}\to {\text{CS}}_2+\text{H}.$$(5)

This pathway represents 99.4% of the flow from CH₄ to CS₂, resulting in around 3.8 × 10⁵ molecule cm⁻³s⁻¹. The flows between each intermediate species involved in these steps are greater than 80% attributable to these reactions, except for the second, which contributes equally with the reaction CH₃ + S₂ → CH₂S + SH to the formation of CH₂S. The sulfur atom and disulfur species originate from the decomposition of the thermodynamically stable species H₂S through the reactions H₂S + H → SH + H₂, SH + H → S + H₂, and CS₂ + S → CS + S₂ as CS₂ becomes abundant. The limiting step is CHS + S₂ → CS + HS₂, but the second and third reactions in this pathway, CH₃ + S → CH₂S + H or CH₃ + S₂ → CH₂S + SH and CH₂S + H → CHS + H₂, are also very close to saturation (i.e., they are rate-limiting, with 99.2 and 92.3%, respectively). The primary sink of CS₂ is its advection to lower atmospheric layers and destruction via oxidation to OCS and subsequently CO through the steps CS₂ + S → CS + S₂, CS + SO → OCS + S, and OCS + H → CO + SH, with the second reaction as the limiting step. This formation pathway is absent from the Tsai 2023 network because CH₂S is not included. This results in a very reduced flow of 1.4 × 10⁻⁶ molecule cm⁻³s⁻¹, which is equivalent to the flows from other carbon-bearing species such as CO, CO₂, HCN, and OCS. We are confident that this pathway is accurate because it arises through reactions that were investigated in our preliminary study on CH₃SH pyrolysis (Veillet et al. in rev). The corresponding experimental data (no. 3 in Table 1), shown in Fig. D.1, were compared to the ideal reactor simulations with the present work and Tsai 2023. These conditions clearly show that the Tsai 2023 kinetic network lacks a formation mechanism to explain the experimental CS₂ abundance observed during the pyrolysis of CH₃SH. In particular, the formation pathway and sensitivity analysis in this study demonstrate that CH₂S is the key species that controls the formation of CS₂. Precise ab initio calculations have been performed on this species, and the reaction CH₃ + S₂ → CH₂S + SH is calculated in Veillet et al. (in rev) with full pressure dependence by solving the master equation for the full potential energy surface. However, many analogies were also used in this work given the sheer number of possible reactions, and further thorough studies are required to fully describe the kinetics, particularly under oxidation conditions. This will be the focus of a follow-up study aiming to improve this oxidation mechanism (Veillet et al., in prep). This species has already been observed as a product of photochemistry experiments (He et al. 2020), and the detection of CS₂ in TOI-270 d (Holmberg & Madhusudhan 2024; Benneke et al. 2024) further supports the likelihood that high concentrations of CS₂ are induced by photochemistry. In addition, CS₂ has been highlighted in the exoplanet photochemistry modeling literature (Mukherjee et al. 2025a; Zahnle et al. 2009), particularly for cooler exoplanets, a trend also observed in the planets studied here. Recent work (Moses et al. 2024) also highlights the same formation pathway for CS₂, but notes that OCS chemistry, which is closely linked to CS₂ oxidation, remains uncertain. This suggests that greater effort is needed to validate networks against the widest possible range of available data to highlight these blind spots in our models.

Fig. 4 shows the effect of coupling sulfur chemistry to C/H/O/N species by comparing the present sulfur network (solid lines) with the previous C/H/O/N network, V23 (dashed lines). For highly abundant species such as H₂O, CH₄, CO, andCO₂, we observe very few consequences. For C₂H₂ and C₂H₄, however, we observe a higher production when including sulfur species in the network. At 10⁻³ bar, for example, C₂H₂ is more than two orders of magnitude more abundant in GJ 436 b, and C₂H₄ is between one and two orders of magnitude more abundant when using the new sulfur network. For GJ 1214 b, the magnitude impact on C₂H₄ is similar, whereas for C₂H₂ the effect is less than one order of magnitude. In the absence of sulfur, the normal formation pathway for C₂H₄ from CH₄ is as follows: $${\text{CH}}_4+\text{H}\to {\text{CH}}_3+{\text{H}}_2$$(6) $${\text{CH}}_3+{\text{CH}}_3\to {\text{C}}_2{\text{H}}_6$$(7) $${\text{C}}_2{\text{H}}_6+\text{H}\to {\text{C}}_2{\text{H}}_5+{\text{H}}_2$$(8) $${\text{C}}_2{\text{H}}_5\to {\text{C}}_2{\text{H}}_4+\text{H}.$$(9)

However, in the presence of sulfur, this pathway accounts for only 2.5% of the C₂H₄ produced. Most of the production skips the second step CH₃ + CH₃ →C₂H₆, which is the limiting reaction, and proceeds directly to C₂H₄ through the reaction CH₂S + CH₃ → C₂H₄ + SH. This reaction was also computed in our CH₃SH pyrolysis study (Veillet et al. in rev), with the full pressure dependence calculated using the same master equation method. This result highlights the benefits of cross-validation between exoplanet modeling and combustion or pyrolysis experiments. The pathway to C₂H₂ with sulfur involves the reaction CH₂S + CHS → C₂H₂ + HS₂, which was an analogy with the equivalent addition CH₂O + HCO → CH₂OCHO. As the accuracy of this analogy is uncertain, further studies are required to confirm or disprove the reality of this reaction.

These abundance differences have very little impact on the synthetic transmission spectrum of GJ 436 b between 1 and 20 µm. However, around 25 µm, a new and significant CS₂ feature of approximately 200 ppm is visible. Fig. B.1 further highlights the CS₂ contribution to the spectra, showing that the features around 6.5 µm and 11 µm are also due to CS₂, although they remain too small to significantly outweigh the NH₃ features. For GJ 1214 b, abundance changes due to the chemical scheme have a more significant impact on the spectrum, although much of the spectral differences are due to an increase in the continuum, probably due to the large fraction of heavy elements in the atmosphere. The CS₂ feature at 25 µm is also present, and a new CS₂ feature at 11 µm can be seen next to the C₂H₄ feature.

Fig. 5 Synthetic transmission spectra of GJ 436 b, GJ 1214 b, HD 189733 b, HD 209458 b, WASP-39 b, and WASP-107 b simulated with TauREx at a resolution of 50 for the abundance profiles shown in Figs. 3 and 4.

3.3 Results for HD 189733 b and HD 209458 b

To model HD 189733 b and HD 209458 b, we adopted solar metallicity, as this appears to be a reasonable approximation based on known values (Xue et al. 2024; Finnerty et al. 2024). We also used the same eddy diffusion coefficient profiles as in Veillet et al. (2024). The corresponding abundance profiles of Fig. 3 show similar differences to those observed for GJ 436 b and GJ 1214 b in Sect. 3.2. CS₂ is more abundant with our sulfur kinetic network than with Tsai 2023 by one to two orders of magnitude, due to the formation pathway not present in Tsai 2023. However, while this species reaches more than 1 ppm in HD 189733 b, its abundance remains low in HD 209458 b, peaking at 10⁻⁹ between 10⁻³ and 10⁻² bars. Another notable difference is the CH₄ abundance in HD 189733 b, which is lower in our network than in Tsai 2023. This is because in this layer of the atmosphere, the equilibrium between CS and CS₂ is shifted toward CS. This causes a decrease in the CH₄ abundance around 10⁻⁴ bar, which is not seen in the Tsai 2023 network due to the low abundances of these species. Other differences among the major species are mostly minor, although in the upper atmosphere of HD 209458 b H₂S, CO₂, and H₂O abundances are slightly higher with our sulfur network than in Tsai 2023. Fig. 4 compares the abundance profiles obtained with our new sulfur scheme with those obtained with our network without sulfur. We observe several differences for HD 209458 b, mainly in NH₃, CH₄, C₂H₂, and HCN. However, these differences remain smaller than one order of magnitude. In contrast, for HD 189733 b, multiple C/H/O/N species are affected by the coupling with sulfur chemistry. The CH₄ difference has already been discussed and is linked to the abundance of CS at this altitude. However, this directly impacts the HCN chemistry, which is formed from CH₄ by the sequence of reactions: $${\text{CH}}_4+\text{H}\to {\text{CH}}_3+{\text{H}}_2$$(10) $${\text{CH}}_3+{\text{N(}}^4\text{S)}\to {\text{H}}_2\text{CN}+\text{H}$$(11) $${\text{H}}_2\text{CN}\to \text{HCN}+\text{H}.$$(12)

Remarkably, the reaction CH₃ + N(⁴S) → H₂CN + H resembles the previous reaction CH₃ + S → CH₂S + H. However, since the product is a radical, unlike in the CH₂S case, it rapidly undergoes a β-scission to HCN. This close coupling between CH₄ and HCN reduces the HCN abundance by approximately two orders of magnitude when accounting for sulfur, due to its impact on the CH₄ abundance. The differences in the NH₃ abundance profile are directly linked to this formation of HCN, as it also requires the nitrogen atoms that are primarily derived from NH₃. At 10⁻⁴ bar in HD 189733 b, another notable coupling effect of sulfur chemistry is the main pathway of HCN formation from ammonia, which proceeds through NS, SNO, NO, HCNO, and HCN as follows: $${\text{NH}}_3+\text{H}\to {\text{NH}}_2+{\text{H}}_2$$(13) $${\text{NH}}_2+\text{H}\to \text{NH}+{\text{H}}_2$$(14) $$\text{NH}+\text{H}\to {\text{N(}}^4\text{S)}+{\text{H}}_2$$(15) $$\text{SH}+{\text{N(}}^4\text{S)}\to \text{NS}+\text{H}$$(16) $$\text{NS}+\text{OH}\to \text{SNO}+\text{H}$$(17) $$\text{SNO}+\text{H}\to \text{SH}+\text{NO}$$(18) $${}^3{\text{CH}}_2+\text{NO}\to \text{HCNO}+\text{H}$$(19) $$\text{HCNO}+\text{H}\to \text{HCN}+\text{OH}$$(20)

The oxidation of N(⁴S) to NO can also occur through via the reaction N(⁴S) + SO → NO + S. In the Tsai 2023 network, N(⁴S) reacts with CS to form HCN via two steps: N(⁴S) + CS → S + CN and H₂ + CN → H + HCN. Although the effect is small, these reactions provide a clear example of the coupling between sulfur and nitrogen chemistry. The abundance profiles C₂H₂ and C₂H₄ are also affected by sulfur in two ways, due to the new formation pathway discussed in Sect. 3.2. Sulfur increases the abundances of these species near 10⁻² bar, while the reduced CH₄ abundance decreases them, dominating above 10⁻⁴ bar. For pressures lower than this, sulfur is mostly in the atomic form on both planets; however, there is a narrow pressure range around 10⁻³ bar where S₂ dominates for HD 189733 b, while SH radicals dominate for HD 209733 b.

Looking at the synthetic transmission spectra of HD 209458 b and HD 189733 b (Fig. 5), we see that for HD 209458 b there is no significant change, except for slight deviations, particularly visible in the continuum above 10 µm. This can be understood from the spectra contributions (Fig. B.1), as the spectrum of HD 209458 b is entirely dominated by water. For HD 189733 b, however, several differences are notable. First, the CS₂ feature at 25 µm is visible in the spectrum obtained with the abundances from our sulfur kinetic network. The second CS₂ feature at 11 µm is also visible, with an amplitude of approximately 100 ppm. Second, the HCN feature at 13 µm disappears from the spectrum with our network compared to the Tsai 2023 network. Third, two amplitude changes occur at 3 and 7 µm, respectively. These correspond to CH₄ features, which show higher amplitudes in the Tsai 2023 and V23 networks. This difference is due to the lower CH₄ abundance obtained with our C/H/O/N/S network due to carbon-sulfur coupling, which causes fainter CH₄ features.

3.4 Results for WASP-39 b and WASP-107 b

In order to stay as close to the original simulations as possible, we chose to model WASP-39 b and WASP-107 b with 10× solar metallicities, the eddy diffusion profile used in Tsai et al. (2023) for WASP-39 b, and a constant eddy coefficient of 10¹⁰ cm ²/s for WASP-107 b, as in (Dyrek et al. 2024). These two planets are the first in which the JWST has detected a sulfur species, namely SO₂. This raises the question of how the differences between our network and Tsai 2023 impact the SO₂ abundances and the corresponding synthetic spectra. For WASP-39 b, the abundance profiles in Fig. 3 show that the SO₂ abundance profiles are very similar. This is because, in both kinetic networks, the formation mechanism of SO₂ is the same, $${\text{H}}_2\text{S}+\text{H}\to \text{SH}+{\text{H}}_2$$(21) $$\text{SH}+\text{H}\to \text{S}+{\text{H}}_2$$(22) $$\text{S}+\text{OH}\to \text{SO}+\text{H}$$(23) $$\text{SO}+\text{OH}\to {\text{SO}}_2+\text{H,}$$(24) with the limiting step being the last one. This corresponds to the mechanism presented in Tsai et al. (2023), although the reaction rates are not exactly from the same source, which causes the observed differences. This high concentration of SO₂ arises from the dissociation of H₂S, which mainly forms S₂ between 10⁻³ and 10⁻⁴ bar, and sulfur atoms below 10⁻⁴ bar. SO₂ is not the only species that differs between the networks, we also observe that the CS₂ abundance is higher in the present network than in Tsai 2023, as previously explained for all planets. This high concentration of sulfur atoms also impacts both NH₃ and CH₄. Indeed, N(⁴S) atoms are produced via the same mechanism as described when discussing the formation of HCN in Sect. 3.3, and they can then be oxidized by SH radicals through SH + N(⁴S) → NS + H. This reaction accounts for the differences in NH₃ between the present sulfur network and Tsai 2023 seen around 10⁻⁴ bar for WASP-39 b in Fig. 4. For CH₄, the abundance of SO oxidizes the carbon to CO₂ via the following mechanism: $${\text{CH}}_4+\text{H}\to {\text{CH}}_3+{\text{H}}_2$$(25) $${\text{CH}}_3+\text{S}\to {\text{CH}}_2\text{S}+\text{H}$$(26) $${\text{CH}}_2\text{S}+\text{H}\to \text{CHS}+{\text{HS}}_2$$(27) $$\text{CHS}+{\text{S}}_2\to \text{CS}+{\text{HS}}_2$$(28) $$\text{CS}+\text{SO}\to \text{OCS}+\text{S}$$(29) $$\text{OCS}+\text{H}\to \text{CO}+\text{SH}$$(30) $$\text{CO}+\text{OH}\to {\text{CO}}_2+\text{H}.$$(31)

Both reactions CHS + S₂ → CS + HS₂ and CH₃ + S → CH₂S + H are the limiting steps. The effects of this pathway are more evident when comparing the CH₄ abundances from the C/H/O/N kinetic network (dashed lines) with those obtained using the sulfur network (solid lines). This effect is also observed with Tsai 2023, although it does not include the CH₂S species, because this oxidation can proceed through the following alternate pathway from CH₃ to CS: $${\text{CH}}_3+\text{H}{\to }^1{\text{CH}}_2+{\text{H}}_2$$(32) $${}^1{\text{CH}}_2{\text{+H}}_2{\to }^3{\text{CH}}_2{\text{+H}}_2$$(33) $${}^3{\text{CH}}_2\text{+S}\to \text{CHS+H}$$(34) $$\text{CHS+H}\to {\text{CS+H}}_2.$$(35)

The oxidation of CS can then proceed through the same steps.

For WASP-107 b, the situation is very different, as the planet is much cooler. This allows for a huge production of CS₂ with our sulfur kinetic network, which predicts CS₂ abundances that are over seven orders of magnitude higher than the Tsai 2023 network. Other small differences are observed for SO₂ and CH₄, for the same reasons discussed above. The abundance of H₂S is also affected by variations in the abundance of CS₂. The production of C₂H₂ and C₂H₄ is strongly favored from the CH₄ species through the following pathway: $${\text{CH}}_4\text{+H}\to {\text{CH}}_3{\text{+H}}_2$$(36) $${\text{CH}}_2{\text{S+CH}}_3\to {\text{C}}_2{\text{H}}_4\text{+SH}$$(37) $${\text{C}}_2{\text{H}}_4\text{+H}\to {\text{C}}_2{\text{H}}_3{\text{+H}}_2$$(38) $${\text{C}}_2{\text{H}}_3\to {\text{C}}_2{\text{H}}_2\text{+H,}$$(39) with the last step being the limiting one. It is also worth noting that the reaction CH₂S + CH₃ → C₂H₄ + SH is only approximately 7.3% saturated in this pathway, so it should not be very sensitive to the value of the rate constant; however, its inclusion in the kinetic network is important to allow bypassing the formation of the C₂H₆ intermediate species, and proceeding directly to the C₂H₄ product. At 10⁻³ bar, this represents about 55.7% of C₂H₂ formation, whereas the standard pathway without sulfur accounts for 14%. A third formation pathway, which contributes 12.1% is noteworthy because it directly involves the species CH₃SH, $${\text{CH}}_4\text{+H}\to {\text{CH}}_3{\text{+H}}_2$$(40) $${\text{OCS+CH}}_3\to {\text{CO+CH}}_3\text{S}$$(41) $${\text{CH}}_3{\text{S+H}}_2\text{S}\to {\text{CH}}_3\text{SH+SH}$$(42) $${\text{CH}}_3\text{SH+SH}\to {\text{CH}}_2{\text{SH+H}}_2\text{S}$$(43) $${\text{CH}}_2\text{SH+SH}\to {\text{H}}_2{\text{S+CH}}_2\text{S}$$(44) $${\text{CH}}_2\text{S+CHS}\to {\text{C}}_2{\text{H}}_3{\text{+S}}_2$$(45) $${\text{C}}_2{\text{H}}_3\to {\text{C}}_2{\text{H}}_2\text{+H}.$$(46)

The limiting step of this pathway is the reaction CH₂S + CHS → C₂H₃ + S₂, which is an analogy made in the network with the addition CH₂O + HCO → CH₂OCHO;therefore, further work on this reaction is needed to confirm this C₂H₂ formation mechanism. However, it involves two more reactions that we calculated in Veillet et al. (in rev), which are CH₃SH + SH → CH₂SH + H₂S and CH₃S + H₂S → CH₃SH + SH. These reactions are also key in CH₃SH pyrolysis, and our extensive work and validation on experimental data covering these conditions provide high confidence in their rates. For other species such as CH₄, the same mechanisms observed for previous planets apply, which are responsible for its oxidation to CO₂, the difference between the C/H/O/N and the C/H/O/N/S network in HCN abundance around 10⁻⁵ bar, and the difference in NH₃ around the same pressure.

The transmission spectra of these planets are shown in Fig. 5. For WASP-39 b, the impact of sulfur coupling on the transmission spectra is negligible. The two features not predicted with the C/H/O/N network are the H₂S feature at 4 µm and the SO₂ feature at 7.3 µm. Despite the few differences between the networks, both predict a very similar spectrum for this planet.

For WASP-107 b, however, many different features are visible. These features include CS₂ at 25 and 11 µm, as well as a different amplitude of about 2000 ppm on the CH₄ feature at 3.3 µm, due to the coupling between sulfur and carbon chemistry, which reduces its abundance. The SO₂ feature at 7.3 µm also changes the amplitude between models by about 4000 ppm. This indicates that the transit depth is very sensitive to SO₂ abundance, as the SO₂ differences in abundance are relatively minor between V25 and Tsai 2023. The differences in NH₃ abundance also cause the disappearance of the NH₃ feature at 20 µm, which is not present in the spectrum corresponding to our sulfur kinetic network but has an amplitude of about 3000 ppm with Tsai 2023. For both planets, we compared our synthetic spectra with observation data, primarily from the JWST. For WASP-39 b, our synthetic spectra appear to correctly reproduce these observational data, but the SO₂ feature at 4.1 µm is not observed. This difference is not related to the kinetic network, as we do not reproduce this feature with the Tsai 2023 kinetic network. Since we also used the same planet parameters as in their original work, for which this feature is the main highlight, the absence of this feature is likely linked to differences in the physical modeling of the planet. These differences may arise from both the use of the 1D model FRECKLL and the radiative transfer model TauREx. If due to the 1D model, differences might lie in the simulation parameters, such as the UV flux incident angle or differences in the altitude profile, which is not recomputed from hydrostatic equilibrium while solving for the steady state in the current version of FRECKLL. However, comparison between different radiative transfer models of synthetic spectra generated from the same FRECKLL output shows that some are able to reproduce this SO₂ feature at 4.1 µm, strongly suggesting that the differences mainly stem from the radiative transfer modeling. However, for WASP-107 b, using a clear-sky model as presented here cannot reproduce the observed spectrum, which is unsurprising given that silicate clouds have already been proposed to explain observations (Dyrek et al. 2024). We tested simple, optically thick cloud models with TauREx, but identifying suitable cloud pressures and photospheric radii proved too challenging to achieve a satisfactory match between our modeled synthetic spectra and observations. This difficulty also appears related to the puffiness of WASP-107 b, for which the synthetic spectrum was very sensitive to these parameters compared to WASP-39 b. More detailed comparisons with observational data would require careful modeling and a search for an optimal value of these parameters, beyond the scope of this paper, which focuses on the modeling of the chemical kinetics. A more comprehensive parameter study may be conducted in a follow-up paper.

4 Conclusions

In this work, we developed the first sulfur kinetic network extensively validated on combustion and pyrolysis data to model exoplanet disequilibrium chemistry. The network builds on our previous C/H/O/N work (Veillet et al. 2024) and kinetic networks from the combustion literature, which we specifically expanded for CH₃SH pyrolysis kinetics with high-level ab initio calculations and analogies to explore a large number of new reactions and yield the most comprehensive sulfur network possible (Veillet et al. in rev). Special emphasis was placed on the coupling of sulfur with all other atoms, including carbon and nitrogen. The coupling between carbon and sulfur has received extensive care and study due to available experimental data for heteroatomic compounds containing sulfur, carbon, and hydrogen together. These kinds of kinetic experimental data are vital for improving our understanding of the reactions involved in these conditions. The species CS₂ has received particular attention because it forms during CH₃SH pyrolysis. However, existing networks and experiments in the literature mainly focus on CS₂ combustion and on pyrolysis, which does not involve hydrogen and therefore cannot be the only source of kinetic data for networks used to model hydrogen-dominated hot exoplanets. We modeled six planets: GJ 436 b, GJ 1214 b, HD 189733 b, HD 209458 b, WASP-39 b, and WASP-107 b, two of which were proven to have at least one sulfur signature species in their atmosphere distinctive of photochemistry, SO₂. The simulations highlight a substantial contribution of the CS₂ species to the abundance profiles, with abundances as high as 1000 ppm. We computed the corresponding synthetic transmission spectra, which show that this species can contribute to the spectra of the cooler planets. This is in good agreement with recent JWST observations of TOI-270 d (Holmberg & Madhusudhan 2024; Benneke et al. 2024), which indicate potential signatures of CS₂ in the spectrum. This finding is clearly relevant for the exoplanet community, as it could potentially be a photochemistry marker for cooler planets, analogous to SO₂, and be closely linked to its metallicity. However, some care should be taken, as the full picture of CS₂ in exoplanets is far from complete, and further validation should be performed by exploring CS₂ oxidation in a hydrogen-heavy medium. This species is not the only one of interest for a reliable and fully coupled C/H/O/N/S network. Other species of interest, such as C₂H₂ and C₂H₄, are also affected by the addition of sulfur chemistry, with an increase in abundance of up to two orders of magnitude. Other species, such as CH₄, NH₃, and HCN, are also found to be impacted by sulfur chemistry, mainly through the central species CH₂S, which is key in the coupling of carbon and sulfur chemistry. may contribute significantly to the chemistry through formation from the NS radical, indicating a possible role of nitrogen in the coupling with sulfur chemistry. However, additional work is required to develop a comprehensive network for nitrogen-sulfur chemistry. These results clearly highlight the need for more comprehensive kinetic networks and for extensive validation of our chemistry models in the JWST era. The best tools available for this task are combustion and pyrolysis experiments, which remain the best way to improve the reliability of our photochemistry and thermochemistry models of exoplanet atmospheres.

Data availability

The kinetic network is available on Zenodo (https://doi.org/10.5281/zenodo.17554652), on KIDA (https://kida.astrochem-tools.org/) and on the ANR EXACT website (https://www.anr-exact.cnrs.fr/fr/chemical-schemes/).

Acknowledgements

This project is fund by the ANR project ‘EXACT’ (ANR-21-CE49-0008-01). In addition, O.V. acknowledges funding from the Centre National d’Études Spatiales (CNES). High performance computing resources were provided by IDRIS under the allocation AD010812434R3 made by GENCI and also by the EXPLOR center hosted by the University of Lorraine. This work was partly supported by a Science and Technology Facilities Council Small Award [ST/Y00261X/1].

References

Appendix A Photodissociation data

Appendix B Transmission spectrum contributions

Fig. B.1 Contribution to the synthetic transmission spectra for all considered absorbing species for the abundances profiles calculated with the present work kinetic network. The plotted quantity is the relative difference to the total transit depth, $\frac{D_{tot}-D_{spec}}{D_{tot}}$, with Dtot being the total transit depth and Dspec being the contribution of the species to the total transit depth. For each wavelength, only the uppermost and close lines have significant impact on the spectrum.

Appendix C Comparisons on CH₃SH pyrolysis data

Fig. C.1 Comparison between the kinetic networks from present work (green), Colom-Díaz et al. (2021) (blue) and Tsai et al. (2023) on experimental CH3SH pyrolysis data (N°3 in Table 1) from Alzueta et al. (2019).

Appendix D Benchmark on the entire combustion and pyrolysis experiments dataset

Fig. D.1 Distribution of each kinetic network prediction errors on the entire experimental dataset, computed as detailed in Veillet et al. (2024). Data is grouped between different conditions of each reference publication in the bataset, based on the fuel studied. The resulting species groups are therefore H2S from Stagni et al. (2022) for conditions 1 to 2 (in cyan), CH3SH from Alzueta et al. (2019) for conditions 3 to 6 (in yellow), CS2 from Zeng et al. (2019) for conditions 7 to 9 (in red), OCS from Zeng et al. (2021) for conditions 10 to 12 (in green), CS2 and OCS with H2O from Abián et al. (2015) for conditions 13 to 15 (in blue) and 16 to 18 (in black). Each condition reference number can be found in Table 1. CH3SH data with water are included in the same group as dry conditions to reduce visual clutter.

All Tables

All Figures

	Fig. 1 Comparison between models from the literature for the combustion of 500 ppm of H2S at atmospheric pressure and an air-to-fuel equivalence ratio of four in a flow reactor (no. 1 in Table 1).
In the text

	Fig. 2 Comparison between models from the literature for H2S pyrolysis (500 ppm at atmospheric pressure in a flow reactor, no. 2 in Table 1).
In the text

	Fig. 3 Abundance profiles of GJ 436 b, GJ 1214 b, HD 189733 b, HD 209458 b, WASP-39 b, and WASP-107 b simulated with FRECKLL. Solid lines show results from the present work, while dashed lines show those from Tsai 2023. We omit H2 and He for clarity to focus on the main species.
In the text

	Fig. 4 Abundance profiles of GJ 436 b, GJ 1214 b, HD 189733 b, HD 209458 b, WASP-39 b, and WASP-107 b simulated with FRECKLL. Solid lines show results from the present work, while dashed lines show results from V23.
In the text

	Fig. 5 Synthetic transmission spectra of GJ 436 b, GJ 1214 b, HD 189733 b, HD 209458 b, WASP-39 b, and WASP-107 b simulated with TauREx at a resolution of 50 for the abundance profiles shown in Figs. 3 and 4.
In the text

Fig. B.1 Contribution to the synthetic transmission spectra for all considered absorbing species for the abundances profiles calculated with the present work kinetic network. The plotted quantity is the relative difference to the total transit depth, $\frac{D_{tot}-D_{spec}}{D_{tot}}$, with Dtot being the total transit depth and Dspec being the contribution of the species to the total transit depth. For each wavelength, only the uppermost and close lines have significant impact on the spectrum.

In the text

	Fig. C.1 Comparison between the kinetic networks from present work (green), Colom-Díaz et al. (2021) (blue) and Tsai et al. (2023) on experimental CH3SH pyrolysis data (N°3 in Table 1) from Alzueta et al. (2019).
In the text

Fig. D.1 Distribution of each kinetic network prediction errors on the entire experimental dataset, computed as detailed in Veillet et al. (2024). Data is grouped between different conditions of each reference publication in the bataset, based on the fuel studied. The resulting species groups are therefore H2S from Stagni et al. (2022) for conditions 1 to 2 (in cyan), CH3SH from Alzueta et al. (2019) for conditions 3 to 6 (in yellow), CS2 from Zeng et al. (2019) for conditions 7 to 9 (in red), OCS from Zeng et al. (2021) for conditions 10 to 12 (in green), CS2 and OCS with H2O from Abián et al. (2015) for conditions 13 to 15 (in blue) and 16 to 18 (in black). Each condition reference number can be found in Table 1. CH3SH data with water are included in the same group as dry conditions to reduce visual clutter.

In the text