Mechanism and Energetics of Hydrogen Sulfide Thermolysis from Reactive Molecular Dynamics: Cutoff-Radius Effects, Thermochemically Validated Energy Costs, and the Elementary Reaction Network

1. Introduction

The transition toward sustainable energy systems requires simultaneously reducing pollutant emissions and developing clean energy vectors. Hydrogen sulfide (H₂S) occupies a singular position in this context: it is at once one of the most problematic gaseous pollutants of the oil-and-gas industry and a potential source of hydrogen. Large volumes of H₂S are produced during the sweetening of sour natural gas, in petroleum refining, particularly in hydrodesulfurization [1], in the anaerobic digestion of wastewater and organic waste [2], and in geothermal fields [3]. Its acute toxicity, comparable to that of cyanide, together with its corrosivity, makes it both an occupational and an environmental hazard. The dominant treatment, the Claus process, partially oxidizes H₂S to recover elemental sulfur but converts the hydrogen of the molecule into water, discarding its energy value [4].

Above roughly 700 °C, H₂S undergoes thermal decomposition—thermolysis—an endothermic, reversible reaction, H₂S(g) ⇌ H₂(g) + (1/n) Sₙ(g), with S₂ the dominant sulfur vapor species [5]. The reaction is strongly endothermic and thermodynamically limited: at atmospheric pressure the equilibrium conversion does not exceed ~4.2% at 1000 °C and 1 atm, it approaches 25% only above ~1350 °C [5]. Thermolysis thus offers a circular-economy alternative to the Claus process, recovering both sulfur and usable molecular hydrogen, but its viability depends on operating at elevated temperatures and on understanding the kinetics, the mechanism, and the energy costs of the reaction.

The kinetics of H₂S thermolysis have been debated for decades. In the low-to-medium temperature regime (800–1250 °C), tubular flow-reactor studies agree on a first-order dependence on H₂S concentration , whereas in the high-temperature regime (2700–3800 K), shock-tube studies in the presence of an inert third body M report a second-order primary step, H₂S + M → SH + H + M [6,7,8]. Karan and co-workers reconciled the high- and low-temperature data into a unified rate expression [7]. In thermal-plasma reactors, thermal dissociation dominates, with near-complete conversion reported around ~2400 K [9,10]. The identity of the initiation step has itself been controversial, with homolytic S–H fission [6,11] competing in the literature with a spin-forbidden molecular elimination to S(³P) + H₂ [12,13]. Yet the molecular sequence of elementary steps—where exactly the H₂ comes from, and which step limits conversion—cannot be resolved by macroscopic kinetics alone.

Reactive molecular dynamics (RMD) simulations with ReaxFF-type force fields allow the formation and breaking of chemical bonds to be observed atom-by-atom with femtosecond resolution [14,15,16], making them an ideal tool to address these questions where direct experimentation is complex. Their application to dilute gaseous systems, however, requires reconsidering the interaction cutoff radius of the potential, originally parametrized for condensed phases. There is thus a twofold knowledge gap—methodological, concerning the appropriate simulation parameters, and applied, concerning the mechanism and energy costs of H₂S thermolysis. The present work addresses both within a single, internally consistent set of simulations: we first establish the cutoff radius and temperature window required for reliable dilute-gas RMD, then quantify the energy costs of dissociation and validate them against independent thermochemistry, and finally reconstruct the complete elementary reaction network with forward/reverse event counts for every step.

2. Materials and Methods

2.1. Reactive Molecular Dynamics with ReaxFF

Reactive molecular dynamics integrates Newton's equations of motion for a system of N atoms, with forces obtained from an interatomic potential in which chemical connectivity emerges dynamically from the instantaneous atomic positions. This work employs the ReaxFF potential [15] in the parametrization of Zhang and van Duin, tuned to describe weak interactions between hydrocarbons and water [17]. ReaxFF describes the energy through the bond order, a continuous function of interatomic distance that tends smoothly to zero as atoms separate, so that all connectivity-dependent energy terms (bond, angle, torsion, over/under-coordination) decay continuously as a bond breaks, while van der Waals and Coulomb interactions act between all atom pairs. Partial charges are recomputed at each step by charge equalization (QEq) [18]. Simulations were run in LAMMPS [19] in the isothermal–isobaric (NPT) ensemble using the equations of motion of Shinoda and co-workers [20], with an integration time step of 0.1 fs. Each system consisted of 1000 H₂S molecules (3000 atoms) at 1 atm.

2.2. The Cutoff Radius for Dilute Gases

Both the connectivity-dependent and the non-bonded energy terms are truncated at a maximum interatomic distance, the cutoff radius (rC), with a seventh-order taper function ensuring that the energy and its first derivatives vanish continuously at rC. In condensed systems a cutoff of 10 Å, the value for which the ReaxFF potentials were originally parametrized [17,21,22] are sufficient to describe structural properties of condensed phases, but realistic calculations require the consideration of long-ranged electrostatic contributions [23]. Classical molecular dynamics simulations have shown that the accountability of long-ranged contributions are important to describe interfacial properties of systems interaction through Van der Waals and electrostatic + Van der Waals interactions [24,25]. In a dilute gas, the mean intermolecular distance can considerably exceed this radius, so a 10 Å cutoff may systematically exclude non-bonded interactions relevant to the thermodynamics and dynamics. Cutoff radii of 10, 15, 18, 30 and 40 Å were therefore evaluated, first at 298.15 K against the experimental density [26], and then at 3000 K against the dissociation kinetics. All production runs used rC = 40 Å at 2500, 3000 and 3500 K, temperatures chosen to span the interval over which thermolysis is progressively activated; each was followed for 130–167 ns of effective chemistry after a 10 ns isothermal equilibration at 298.15 K and 1 atm; times are reported relative to the start of the reactive window.

2.3. Reconstruction of the Reaction Network by Co-Occurrence Signatures

Extracting the mechanism from the trajectory requires identifying every elementary reaction event from atom-resolved species-transition records. Each sulfur and hydrogen atom was assigned a species code equal to 10·nS + nH, where nS and nH are the numbers of sulfur and hydrogen atoms in its molecule (12 = H₂S, 11 = HS•, 1 = H•, 2 = H₂, 10 = S, 13 = [H₃S]•, 22 = HSSH). An elementary reaction modifies several atoms in the same time step, and this temporal co-occurrence constitutes its signature; recognizing each signature allows the forward and reverse events of every step to be counted separately, and degenerate exchanges that produce no net chemistry to be excluded. Channels that form a bound pair (the S₂H₂ complex/HSSH, or H₂ from H•+H•) were counted as reactive events only when the product persisted for at least 10 ps; more fleeting encounters were classified as transient collision complexes and excluded, a criterion applied consistently throughout (including to the first appearance of sulfur chains in Section 3.6). As quality control, mass balances were propagated for every atom, and exact closure of the sulfur and hydrogen inventories confirmed that no transitions were missed. The growth of polysulfur chains was followed with a dedicated cluster-detection program based on a union–find algorithm [27,28,29], which reconstructs frame by frame which atoms belong to the same sulfur aggregate.

2.4. Energy Costs and Thermochemical Validation

The energy cost of dissociation was quantified as the slope of the total system enthalpy against the number of H₂S molecules consumed, dH/dN, evaluated over the chemically controlled window (excluding the first 2–3 ns of thermal expansion after the instantaneous heating). The slope was obtained by block averaging (18 blocks) with bootstrap error estimation [30,31], a procedure that is robust to the choice of kinetic-regime boundaries. For independent validation, the complete-dissociation enthalpy was estimated by Kirchhoff's law using the experimental H–SH bond dissociation energy of 376.24 ± 0.05 kJ·mol⁻¹ [32], corroborated by high-level ab initio calculations [33], together with the heat-capacity polynomials of the species [34], yielding 4.11 – 4.12 eV per H₂S molecule over the simulated temperature range.

3. Results and Discussion

3.1. The Cutoff Radius Is Critical for Dilute-Gas RMD

Simulations at 298.15 K and 1 atm, started from an out-of-equilibrium configuration, showed that the cell density first overshoots and then relaxes toward the experimental value for all cutoff radii, but the overshoot is much larger (a peak of ≈4× the equilibrium density, versus ≈2×) and the relaxation slower (~1.7 ns versus ~0.4–0.6 ns) with the original 10 Å cutoff than with 15 or 18 Å (Figure 1). The energetics display the same pattern: the 15 and 18 Å models reach a steady total enthalpy within ~0.2 ns and agree with each other to within 1 eV, whereas the 10 Å model requires ~13 ns and settles ≈23 eV above the converged value (about 0.023 eV per molecule)—that is, it misses roughly 23 eV of attractive non-bonded interaction energy—with instantaneous fluctuations ~1.5 times larger (standard deviation 21 versus 14.5 eV; Figure 2). A significant fraction of the energetically relevant non-bonded interactions therefore occurs beyond 10 Å in the dilute gas.

The consequences for reactivity are more severe. Upon instantaneous heating to 3000 K, no dissociation was observed with rC of 10 or 15 Å within 20 ns, whereas radii of 18, 30 and 40 Å captured the thermolysis, with a dissociation rate that increases with rC but converges: at 130 ns the conversion rises by 2.5 percentage points from 18 to 30 Å but by only 1.0 point from 30 to 40 Å (Figure 3). The cell volume proved essentially insensitive to rC, indicating that the volumetric properties are governed by short-ranged interactions [35], whereas the enthalpy tracked the rC-dependent reaction progress. All subsequent results therefore employ rC = 40 Å.

3.2. Temperature Threshold and Global Kinetics

With rC = 40 Å, no dissociation was detected between 1000 and 2000 K within 20 ns, consistent with the seconds-scale residence times required experimentally at low temperature [7] and with the intrinsic time-scale limitation of RMD. At 2500, 3000 and 3500 K dissociation proceeded clearly (Figure 4). The decay of H₂S and the growth of H₂ and HS• follow exponential asymptotic behavior consistent with first-order kinetics over the conversions reached, in agreement with the flow-reactor literature [7,8]. By the end of the runs (130–167 ns), the net conversion referred to the initial 1000 molecules was 8.7% at 2500 K, 26.7% at 3000 K and 46.3% at 3500 K—an increase of more than a factor of five across the 1000 K interval. The product species (H•, HS•, H₂, S, with HSSH and HS₂• as intermediates) coincide with those of the homolytic mechanisms reported for thermal plasma at comparable temperatures [9], supporting the physical validity of the reactive potential.

3.3. Energy Costs Rise with Temperature and Approach the Kirchhoff Limit

Figure 5a shows the change of the total system enthalpy plotted against the number of H₂S molecules consumed, after excluding the initial 2–3 ns thermal-expansion transient. The relationship is strikingly linear at all three temperatures (R² = 0.990, 0.995 and 0.988), demonstrating that the enthalpy cost per molecule is a single well-defined quantity at each temperature. The block-averaged slopes are 2.41 ± 0.07, 3.03 ± 0.06 and 3.89 ± 0.18 eV per H₂S consumed at 2500, 3000 and 3500 K. The cost increases approximately linearly with temperature, at ≈1.36 eV per 1000 K (Figure 5b), and approaches the complete-dissociation reference of 3.90 eV at 0 K [32], obtained independently from Kirchhoff's law: at 3500 K the simulated cost (3.89 ± 0.18 eV) is statistically indistinguishable from the 0 K value, but slightly shorter than Kirchhoff’s law predictions at this temperature. An extrapolation of a linear regression of the 3 values matches the Kirchhoff’s law prediction at ≈ 3781 K, and beyond this temperature it is expected an asymptotic behavior.

The physical interpretation is direct. At 2500 K a large fraction of the H₂S consumed proceeds only as far as HS• (a single S–H bond broken), so the net cost per molecule lies well below the complete-dissociation value. As the temperature rises, dissociation proceeds further per molecule consumed—more H₂ is released and sulfur begins to condense—and the cost per molecule converges toward the full thermochemical value. The linear extrapolation of the trend reaches the Kirchhoff limit near 3781 K; given that only three temperatures support the fit, this figure should be read as indicative, and the expected behavior beyond the simulated range is a saturation at the thermochemical limit rather than a crossing. This quantitative agreement between the simulated slopes and the independent Kirchhoff estimate validates the RMD energetics once the thermal-expansion contribution is excluded (first 2-3 ns of simulation).

3.4. The Elementary Reaction Network Is Temperature-Invariant

The co-occurrence signature analysis reconstructed the complete elementary reaction network (Figure 6). The fundamental structural finding is that the same network operates at the three temperatures: no new elementary step appears at high temperature that was not already present, if only marginally, at low temperature. What changes is the relative frequency of the higher-barrier channels, which are progressively activated as the temperature increases. Table 1 summarizes the forward/reverse event counts.

Three trends stand out. First, homolysis dominates in number at all temperatures and lies very close to equilibrium, with almost identical forward and reverse counts—settling, for these conditions, the historical question of the initiation step in favor of homolytic S–H fission [8,11]; the molecular-elimination channel does appear, but since ReaxFF evolves on a single effective potential-energy surface and does not track spin, its rate cannot be compared quantitatively with the spin-forbidden channel discussed in the shock-tube literature [12,13]. Second, the channels of elimination, thiyl recombination and disproportionation, marginal at 2500 K, grow by roughly an order of magnitude at 3000 and 3500 K. The branching pathway S + H₂S ⇌ 2 HS• proceeds through the bound S₂H₂ complex: 43, 482 and 822 complete forward passages were observed at the three temperatures, contained within the recombination/disproportionation counts of Table 1. Third, each forward/reverse pair is balanced to within a small percentage—the signature of partial equilibrium at high temperature.

3.5. Hydrogen Originates from Abstraction, Not from Recombination

The data are unambiguous regarding the source of molecular hydrogen. Persistent recombination H• + H• → H₂ (product surviving ≥10 ps) occurred only 1, 13 and 17 times at 2500, 3000 and 3500 K—below 0.5% of the abstraction count at every temperature. Fleeting H–H contacts that redissociated within a few picoseconds were two orders of magnitude more frequent, underscoring the third-body requirement. Essentially all the H₂ is produced by hydrogen abstraction, H• + H₂S → H₂ + HS•, both directly and through a transient hypervalent adduct [H₃S]•, in quantitative agreement with the canonical propagation step of the H/S system [36,37,38,39]. The reason is mechanistic: H+H recombination is a three-body process that cannot compete when each hydrogen atom is surrounded by an enormous excess of H₂S, whereas abstraction is bimolecular, first-order in H• and H₂S, and requires no third body because the HS• fragment carries away the excess energy. Because H₂ arises from abstraction, appreciable H₂ is produced even when the population of free H• is small, as observed at low temperature.

The temperature dependence of the rate-controlling structure was resolved into three successive stages: initiation (homolysis far from equilibrium, radical reservoir accumulating), propagation (abstraction activates and rapidly equilibrates, H₂ rises to a plateau), and partial equilibrium (all reversible channels balanced, net conversion sustained only by the drainage of sulfur into chains). At 2500 K the mechanism reduces to a two-step chain: the population of free H• peaked at 11 atoms and ended at 5—a fraction of 10⁻³ of the 2000 hydrogen atoms—behaving as a quasi-steady-state intermediate, and the products satisfied the radical balance N(HS•) = 2N(H₂) + N(H•) exactly (87 = 2×41 + 5 at the end; Figure 7, Table 2). At 3500 K the radical reservoir thermalizes (H• ≈ 226), the quasi-steady-state approximation ceases to hold, and the system reaches a thermodynamically controlled partial equilibrium.

Figure 8 compares the main species at the three temperatures on a common scale. The HS• radical closely tracks the consumed H₂S at all temperatures, reflecting its role as chain carrier, while atomic sulfur remains practically absent at 2500 K and accumulates appreciably only at 3500 K—the visual evidence that the sulfur sink opens only at the highest temperatures.

3.6. The Rate-Limiting Step: Sulfur Condensation

S–H bond events are counted in the tens of thousands, whereas S–S bond-forming events number only in the hundreds. Every sulfur atom that becomes fixed in a growing Sₓ chain irreversibly loses its capacity to produce net H₂, so the net production of hydrogen follows the slow condensation of sulfur and not the fast, equilibrated abstraction step. Figure 9 and Table 3 show the first appearance of each sulfur-chain class that persisted for at least 10 ps—a criterion that excludes fleeting collision complexes. The species appear in a strictly sequential order—HSSH, then S₂, then S₃, then S₄—evidencing atom-by-atom chain growth, and every stage accelerates with temperature. Persistent S₃ forms only at 3000 K and above; S₄ was observed only at 3500 K, and only transiently (near 87 ns, surviving less than 10 ps). At 2500 K no sulfur species beyond S₂ ever persisted, which explains the stalling of the conversion at that temperature. The energetics of Section 3.3 tell the same story from the enthalpy side: the cost per molecule approaches the complete-dissociation limit precisely as the sulfur sink opens.

This result has an immediate reading for process engineering. Because the conversion advances only insofar as sulfur is removed from the quasi-equilibrated reservoir into chains, and because this drainage is slow and thermally activated, raising the temperature progressively opens the sulfur sink at the price of an energy cost that converges to the thermochemical limit. Below a certain temperature the process is effectively blocked at the sulfur link, however active the radical chemistry. The engineering lever to improve the yield is therefore to facilitate the removal or condensation of sulfur, rather than to optimize the already-equilibrated hydrogen chemistry.

3.7. Limitations

Three limitations bound the quantitative scope of these results. First, ReaxFF evolves on a single effective potential-energy surface and does not track electronic spin, so the molecular-elimination channel to S(³P) + H₂ appears in the simulations but its rate should not be compared quantitatively with experiment [12,13]. Second, classical nuclear dynamics neglects tunneling in the hydrogen-abstraction step; at 2500–3500 K the system is far above the crossover temperature and the effect is minor, but it introduces uncertainty in absolute abstraction rates at lower temperatures [38]. Third, the global charge-equalization scheme leaves small residual charges on separated radicals; their electrostatic effect is much smaller than kT and does not alter the identity or ordering of the reaction channels, which rest on bond topology. Finally, the accessible time scale (hundreds of ns) restricts the study to T ≥ 2500 K; the experimental low-temperature regime remains beyond direct reach of RMD.

4. Conclusions

A single, internally consistent set of reactive molecular dynamics simulations resolved the methodology, the energetics and the mechanism of H₂S thermolysis. Methodologically, the interaction cutoff radius proved critical for dilute gases: the original 10 Å parametrization misses ~23 eV of non-bonded interaction energy and fails to capture dissociation at 3000 K, whereas 30–40 Å converges; a thermal threshold for observable dissociation lies between 2000 and 2500 K on RMD time scales. Energetically, the system enthalpy rises linearly with the number of H₂S molecules consumed, defining costs of 2.41 ± 0.07, 3.03 ± 0.06 and 3.89 ± 0.18 eV per molecule at 2500, 3000 and 3500 K that increase by ≈1.36 eV per 1000 K, matching at 3500 K the complete-dissociation limit of 3.90 eV obtained independently from Kirchhoff's law—a quantitative, thermochemically anchored validation of the approach. Mechanistically, thermolysis proceeds by a free-radical chain with temperature-invariant elementary steps: homolytic S–H initiation, hydrogen abstraction as the essentially exclusive source of H₂ (direct H•+H• recombination is negligible), and a rate-limiting sulfur-condensation stage that determines the net conversions, with control shifting from a low-temperature quasi-steady-state chain to a high-temperature partial equilibrium. Beyond its mechanistic value, the analysis indicates that the removal of sulfur—not the hydrogen chemistry—is the bottleneck to be addressed when engineering the valorization of residual H₂S as a source of green hydrogen.