A Predictive Model for Catalytic Methane Pyrolysis

Methane pyrolysis provides a scalable alternative to conventional hydrogen production methods, avoiding greenhouse gas emissions. However, high operating temperatures limit economic feasibility on an industrial scale. A major scientific goal is, therefore, to find a catalyst material that lowers operating temperatures, making methane pyrolysis economically viable. In this work, we derive a model that provides a qualitative comparison of possible catalyst materials. The model is based on calculations of adsorption energies using density functional theory. Thirty different elements were considered. Adsorption energies of intermediate molecules in the methane pyrolysis reaction correlate linearly with the adsorption energy of carbon. Moreover, the adsorption energy increases in magnitude with decreasing group number in the $d$-block of the periodic table. For a temperature range between 600 and 1200 K and a normalized partial pressure range for $\mathrm{H_2}$ between $10^{-1}$ and $10^{-5}$, a total of eighteen different materials were found to be optimal catalysts at least once. This indicates that catalyst selection and reactor operating conditions should be well-matched. The present work establishes the foundation for future large-scale studies of surfaces, alloy compositions, and material classes using machine learning algorithms.


Introduction
Hydrogen is an important energy carrier and chemical feedstock and has the potential to replace fossil hydrocarbons in industry and transportation.However, hydrogen production based on fossil fuels accounted for almost the entire production in 2021, with associated emissions of over 900 Mt of CO 2 .Although emission-free hydrogen production from water electrolysis has accelerated over the past few years, it covers only 0.1% of global demand. 1 Thus, an economically and ecologically feasible bridging technology for hydrogen production is needed.
As suggested by Sánchez-Bastardo et al., a promising option is methane pyrolysis, i.e., thermal decomposition of methane in the absence of oxygen to form solid carbon and hydrogen. 2 4 (g) → 2 H 2 (g) + C (s) (1)   While methane pyrolysis still relies on the extraction of natural gas, the process is CO 2 emission-free.From an energetic point of view, only 37.5 kJ of energy is necessary to produce 1 mol of H 2 via methane pyrolysis, compared to 286 kJ required for water electrolysis. 3Even though a significant reaction yield can be achieved only above 1200 °C, 4 this figure can be substantially reduced by using a catalyst.3][14][15][16] In particular, pyrolysis in liquid bubble column reactors holds a big promise for sustainable process which through its design circumvents catalyst deactivation as the carbon that accumulates on top of the metal bath can be continuously removed from the reactor.
Ab initio methods provide a useful tool for high-throughput screening of material choices.
In this work, a straightforward and flexible model for qualitative comparison of possible catalysts is derived.The reaction will be divided into its elementary steps, and a combination of the microkinetic model and Sabatier analysis will be used to describe the reaction rate as a single function of temperature T , pressure p and adsorption energy of carbon ∆E Ads,C which is applied to thirty different elements.

Model Methane pyrolysis
The methane pyrolysis reaction and its individual reaction steps have been widely studied in the past.
Chemisorbed methane dissociates into a methyl radical (CH 3 ) and a hydrogen atom.
The methyl radical dissociates further into a methylene (CH 2 ), and then into a methine (CH) radical, followed by the last dissociation step to form adsorbed elemental carbon and hydrogen.
While carbon desorbs to form solid carbon, the hydrogen atoms recombine to produce hydrogen gas.
While this reaction path is generally accepted, it is still disputed which reaction step is the slowest and, therefore, rate-determining for the overall reaction.In the 1960s, Kozlov and Knorre concluded that Eq. ( 4) is the rate-limiting step, 18 whereas in 1972, Baker considered it to be carbon diffusion. 19More recently, surface transport phenomena are also suggested. 20mputational studies using Density Functional Theory calculations found Eq. ( 3) to be rate-determining are combined since the energy barrier for initial adsorption of methane is negligibly small. 22ence, the first reaction step is reducing the number of steps to six.Secondly, the role of carbon desorption will be neglected.
Carbon deposition on a catalyst surface is a complex phenomenon that can induce the formation of a wide variety of structures, such as carbon nanotubes, among many others. 23,24us, for practical reasons, the role of carbon in this context is greatly simplified by omitting Eq. ( 7) and setting the activity of carbon a C to 1.

Sabatier analysis
From the Sabatier principle, it follows that stronger adsorption of a reactant, i.e., a higher magnitude of the negative adsorption energies ∆E Ads , correlates with a higher reaction rate constant k.This is valid up to a point where too strong adsorption inhibits product desorption, resulting in a lack of free surface sites θ * . 25Together with the Brønsted relationship that links thermodynamics and kinetics of a reaction, 26 both concurring mechanisms caused by adsorption are responsible for the characteristic volcano-shaped plots of catalytic reaction rates that were first reported in 1969. 27In addition, the findings of Abild-Pedersen et al. 28 showing that the adsorption energy of CH x (x = 1, 2, 3) scales approximately with the adsorption energy of a C atom are applied to simplify the model, resulting in the reaction rate as a function of only the adsorption energy of carbon ∆E Ads, C .
As a result, the system of elementary steps that was established above contains five reaction steps, with the fourth step (Eq.( 6)) being rate-determining based on findings by Fan et al.. 22 Four equilibrium constants can be derived: from Eq. ( 9) : from Eq. ( 4) : from Eq. ( 5) : from Eq. ( 8) : where K i is the equilibrium constant of reaction i, θ j is the surface coverage of species j, p CH 4 and p H 2 are partial pressures of CH 4 and H 2 , respectively.Partial pressures are normalized with respect to the total system pressure.K 4 is missing because it is in a non-equilibrium state.Equilibrium constant and Gibbs free energy are connected via For any intermediate steps where no gas-phase molecules are being adsorbed or desorbed, namely (11) and (12), the entropy contributions are ignored as the changes in entropy are negligible, i.e., Reaction steps (10) and ( 13) require a different approach.As shown by Norskov et al., to a good approximation, adsorbed molecules can be expected to lose all their entropy upon adsorption: 29 Consequently, the same must apply in reverse for desorption.
For Eq. (10) to Eq. ( 13) there are now five unknown variables, i.e. all coverages θ j , with only four equations.Rearranging the system of equations above yields four coverages, which are all functions of θ * .
From Eq. ( 13): Eq. ( 10) with Eq. ( 17): Eq. ( 11) with Eq. ( 18): and Eq. ( 12) with Eq. ( 19): The site conservation rule yields a fifth equation with its general form of where a fraction of j-covered sites can be defined as λ j = θ j /θ * .This leads to and the coverage of free sites can then be expressed as Inserting all terms from above yields the coverage of free surface sites for the system as Now, the only non-equilibrated reaction step (Eq.( 6)) is being considered.For this, a 7 reaction rate expression can be set up as with a C = 1.It is not possible to define an equilibrium constant for Eq. ( 25).However, it is convenient to define an "approach to equilibrium" γ. 29 This approach to equilibrium is a positive number that provides information about whether the reaction proceeds in a forward or backward direction.
γ < 1 : The reaction proceeds in the forward direction.
γ = 1 : The reaction is in equilibrium.
γ > 1 : The reaction proceeds in the backward direction.
Hence, the approach to equilibrium γ 4 is defined as leading to a simplified reaction rate of with k 4 defined as Having only one rate-determining reaction step allows for the important simplification of where γ is the approach to equilibrium for the overall reaction rate, i.e., and p CH 4 and p H 2 are normalized partial pressures of CH 4 and H 2 , respectively, K eq is the equilibrium constant, and ∆G is the Gibbs free energy of the overall reaction.Using these simplifications in Eq. ( 27) leads to the reaction rate as with ∆S A 4 and ∆E A 4 as the reaction entropy and energy of non-equilibrated reaction step 4. While ∆E A 4 follows the transition-state scaling relation that is derived from adsorption energy trends, 28 ∆S A 4 is neglected for step 4 (Eq.( 6)), i.e. ∆S A 4 = 0 because reactants and products stay in an adsorbed state.This leads to the final reaction rate equation Impressively, this model describes the reaction rate as a combination of microscopic quantities, such as adsorption and activation energy, and macroscopic properties, such as temperature and pressure.Importantly, the microscopic quantities can be obtained using ab initio calculations.From R, the volcano-shaped curve is obtained, the maximum position of which corresponds to adsorption energy, which in turn is characteristic of a particular ideal catalyst material.

Methodology
Adsorption, reaction, and activation energies for catalytic methane pyrolysis were calculated using density functional theory (DFT) calculations. 30Values for standard enthalpy and entropy of the total reaction were taken from the National Institute of Standards and Technology (NIST). 31The Vienna Ab initio Simulation Package (VASP) was used for all DFT calculations in this work. 32Exchange-correlation functional parametrized by Perdew, Burke, and Ernzerhof (PBE-GGA) 33 and revised by Hammer et al. 34 to improve chemisorption energetics of atoms and molecules on transition-metal surfaces was employed.
A projector augmented wave (PAW) method was used for electron-ion interactions. 35Magnetic moments for ferro-or paramagnetic substrates were taken from the Materials Project (www.materialsproject.org), an open dataset for properties of inorganic materials. 36A supercell of 3 × 3 × 4 (2 × 2 × 4) primitive cells was used to generate the cubic (hexagonal) slabs for the surface calculations, thus leading to (111), (110), and (0001) surfaces for fcc, bcc, and hcp metals, respectively.These represent the closely packed plane for the considered crystal systems.The additional vacuum of more than 20 Å was inserted to separate the periodic images of the free surfaces.Γ-centered k-mesh of 6 × 6 × 1 was chosen.Default plane wave cutoff energy from the PAW pseudopotentials 37 was used.More detailed information on the calculations carried out can be found in the supporting material.

Adsorption energies
Calculated adsorption energies are presented below.A striking correlation between adsorption energies and position in the d-block of the periodic table can be seen in Fig. 2, where each element is color-coded according to its group on the example of the ∆E Ads,CH vs. ∆E Ads,C presented in Fig. 1.Elements that are not included in the d-block are grayed out.One explanation could be due to the fact that transition metals are mainly characterized by their partially filled d-subshells.Due to the directional character of d-orbitals, the nucleus is weakly shielded, interactions between d-electrons are weak, and the nucleus does not only strongly attract d-but also s-electrons from the next higher s-orbital.This results in relatively high but slowly increasing ionization sidering the enthalpy of hydration, which also decreases in magnitude with a higher group number. 39

Activation energies
Norskov et al. 29 proposed that the adsorption and activation energies should be correlated since the same fundamental physics governs them.Abild-Pedersen et al. 28 showed that in the case of methane, this correlation is exclusively related to the ∆E Ads,C term.The activation corresponding to the dissociation of the CH molecule, the rate-determining step, Eq. ( 6), is needed.The explicit activation energy calculations were carried out for a (111)-surface of fcc-Cu.Details can be found in the supporting material.The resulting relationship is

Reaction energies
The remaining ingredient in the free site coverage, θ * , depends on the equilibrium constants K 1 , K 2 , K 3 , and K 5 .For their evaluation through Eq. ( 14), the Gibbs free energy of a reaction is needed.As discussed in Section "Sabatier analysis", different treatment is applied depending on whether the reactions involve gas-phase molecules.
Similar to the activation energies, reaction energies can also be correlated with adsorption energies.The reaction calculations were again carried out for a (111) surface of fcc-Cu and are detailed in the supporting material.The resulting relationships (Eq.( 15)) are and Reaction energies for CH 4 (g) → CH 3 (g) + H(g) and 2H(g) → H 2 (g) were calculated as

Reaction rate
The reaction rate, R, is now calculated from the adsorption energy trends, reaction, and ac-  Figure 4 shows the elements exhibiting the highest reaction rate according to the Fig. 3 and color-coded according to their global market price as of February 2024 41 (in US Dollar ($) per kg).For economic upscaling of catalytic methane pyrolysis to industrial application, a price consideration will be of great importance.Hence, price-informed volcano plots offer a straightforward way to make strategic decisions for a given set of calculated catalyst surfaces with user-defined operating conditions.
Figure 5 shows R as a function of different temperatures and pressures.The upper plot demonstrates that the slope of log(R) (as a function of ∆E Ads,C ) decreases with increasing 16 temperature while the maximum is shifted to lower magnitudes of adsorption energy.Overall, the rate increases with temperature for otherwise the same conditions, as intuitively expected.However, we note that increasing temperature to increase the yield has its practical and economical limits.The lower graph shows that the pressure also serves as a tool to tune the maximum reaction rate (optimal operating conditions) to other elements.Furthermore, the model enables the investigation of the best-performing catalyst material in the T -p space for a given database of adsorption energies on a set of pure fcc, bcc and hcp metals.For a temperature range between 600 and 1200 K and a partial pressure range for H 2 between 10 −1 and 10 −5 , a total of 18 different metals were found to be optimal catalyst materials at least once.
Finally, we point out that our study demonstrates that the type of catalyst and specific reactor operating conditions should be matched.

Fig. 1
depicts adsorption energies of CH, CH 2 , CH 3 , and H, respectively, as a function of carbon adsorption energy for a given substrate element.In the lower right corners of each graph, the Pearson correlation coefficient and the linear fit function, which will be used for the Sabatier model later, are shown.The

Figure 1 :
Figure 1: Adsorption energies of CH, CH 2 , CH 3 , and H as a function of ∆E Ads,C for fcc, bcc, and hcp substrate structures with close-packed surface planes.

Figure 2 :
Figure 2: ∆E Ads,CH as a function of ∆E Ads,C for different substrates.For d-block elements, the correlation between group number and their adsorption energy is clearly evident.
tivation energies.Since R = R(T, p, ∆E Ads,C ) (cf.Eq. (32)), representative thermodynamic input parameters are chosen: T = 1000 K, p CH 4 = 0.99, and p H 2 = 0.01.The resulting reaction rate is shown in Fig. 3 as a function of ∆E Ads,C , which can be linked with various base metals via the data presented in Fig. 1.

Figure 3 :Figure 4 :
Figure 3: Logarithm of the reaction rate, R, as a function of ∆E Ads,C for T = 1000 K, p CH 4 = 0.99, and p H 2 = 0.01.

Figure 5 :
Figure 5: Top: log(R) as a function of ∆E Ads,C for different temperatures (p CH 4 = 0.99 and p H 2 = 0.01).Bottom: log(R) as a function of ∆E Ads,C for different pressures (T = 1000 K).A selection of elements is displayed on the top axis to illustrate their respective positions.18

Figure 6 :
Figure 6: Arrhenius plot of the ideal catalyst as a function of 1/T and log(p H 2 ) = log(1 − p CH 4 ).Note that the colors have no meaning other than to help distinguish areas with different optimal catalysts.