Treatment of Thermal Non-Equilibrium Dissociation Rates: Application to $\text{H}_{2}$

Alex T. Carroll amoricar@caltech.edu Jacob Wolmer Guillaume Blanquart Department of Mechanical and Civil Engineering, California Institute of Technology, Pasadena, California 91125 Aaron M. Brandis NASA Ames Research Center, Mountain View, California 94035 Brett A. Cruden AMA, Inc., NASA Ames Research Center, Mountain View, California 94035

Abstract

This work presents a detailed description of the thermochemical non-equilibrium dissociation of diatomic molecules, and applies this theory to the case of $\rm H_{2}$ dissociation. The master equations are used to derive corresponding aggregate rate constant expressions that hold for any degree of thermochemical non-equilibrium. These general expressions are analyzed in three key limits/ regimes: the thermal equilibrium limit, the quasi-steady-state (QSS) regime, and the pre-QSS regime. Under several simplifying assumptions, an analytical source term expression that holds in all of these regimes, and is only a function of the translational temperature, $T_{\rm t}$ , and the fraction of dissociation, $\phi_{\rm A}$ , is proposed. This expression has two input parameters: the QSS dissociation rate constant in the absence of recombination, $k_{\rm d,nr}(T_{\rm t})$ , and a pre-QSS correction factor, $\eta(T_{\rm t})$ . The value of $\eta(T_{\rm t})$ is evaluated by comparing the predictions of the proposed expression against existing master equation simulations of a 0-D isothermal and isochoric reactor for the case of $\rm H_{2}$ dissociation with the third-bodies $\rm H_{2}$ , H, and He. Despite its simple functional form, the proposed expression is able to reproduce the master equation results for the majority of the tested conditions. The best fit of $k_{\rm d,nr}(T_{\rm t})$ is then evaluated by conducting a detailed literature review. Data from a wide range of experimental and computational studies are considered for the third-bodies $\rm H_{2}$ , H, and inert gases, and fits that are valid from 200 to 20,000 K are proposed. From this review, the uncertainty of the proposed fits are estimated to be less than a factor of two.

^†^†preprint: AIP/123-QED

I Introduction

The dissociation of diatomic species is a relevant physical process for many fields of research, including astrophysics Draine, Roberge, and Dalgarno (1983); Hollenbach and McKee (1989); Chang and Martin (1991); Bell and Cowan (2018); Kristensen et al. (2023), plasma physics Matveyev and Silakov (1995); Silakov et al. (1996); Capitelli et al. (2002); Shakhatov and Lebedev (2018); Smith et al. (2024), detonations Taylor et al. (2013); Voelkel et al. (2016); Shi et al. (2016); Vargas et al. (2022); Baumgart, Yao, and Blanquart (2025), and hypersonic entry flows Park (2011, 2012); Palmer, Prabhu, and Cruden (2014); Erb, West, and Johnston (2019); Carroll and Brandis (2023); Scoggins, Hinkle, and Shellabarger (2024); Steer et al. (2024); Steuer et al. (2024); Rataczak, Boyd, and McMahon (2024). In each of these applications, large temperature gradients lead to internal energy excitation and dissociation of molecular species. Understanding and accurately modeling these non-equilibrium phenomena is complicated by the fact that the processes of internal energy relaxation and dissociation are inherently coupled. For example, high internal energy states are preferentially dissociated, which lead to non-Boltzmann state distributions Marrone and Treanor (1963); Macheret et al. (1994).

To obtain a detailed description of these non-equilibrium phenomena, ab-initio computational methods can be used. The most popular of these methods involves first constructing accurate potential energy surfaces (PESs) that are then used in quasi-classical trajectory (QCT) calculations to obtain state-specific cross-sections/ rate constants that explicitly describe the internal energy excitation and dissociation processes. An overview on the earliest uses of the QCT method is provided by Porter Porter (1974). These state-specific rates can then be used in master equation simulations, where the transitions/ reactions between internal energy states are explicitly computed. A related method is the direct molecular simulation (DMS) approach proposed by Koura Koura (1997), where trajectory calculations on a given PES are performed during a flow simulation.

Recently, several studies in the hypersonic flow literature have used either the master equation approach Kim, Kwon, and Park (2009); Kim and Boyd (2012); Panesi et al. (2013); Andrienko and Boyd (2015) or the DMS approach Norman, Valentini, and Schwartzentruber (2013); Valentini et al. (2014); Schwartzentruber, Grover, and Valentini (2018) to study non-equilibrium dissociation at high temperatures (> 10,000 K). These studies have highlighted several key physical features in a dissociating post-shock flow. Namely, the internal state distribution of the dissociating molecule initially begins at a Boltzmann distribution defined at the pre-shock temperature. Then, excited internal energy states become populated, eventually leading to a local balance between excitation and dissociation known as the quasi-steady-state (QSS) phenomena. Finally, the distribution evolves towards thermochemical equilibrium, defined by a Boltzmann distribution at the post-shock temperature.

While the physical insights gained from these ab-initio studies are valuable, they are unfortunately too computationally expensive to use in most practical computational fluid dynamics (CFD) calculations. The best strategy for incorporating the necessary physics from these detailed approaches into computationally tractable kinetic models for CFD calculations is still an open question. Some models, such as the Park two-temperature Park (1988, 1989a, 1989b), the modified Marrone-Treanor (MMT) Chaudhry and Candler (2019); Chaudhry et al. (2020) and the Singh-Schwartzentruber Singh and Schwartzentruber (2020a, b) models, rely on a multi-temperature framework where rates are described as a function of the translational temperature and one or more internal energy temperatures. Others Magin et al. (2012); Sahai et al. (2017); Venturi et al. (2020) rely on the more general multi-group/ coarse-graining approach described by Liu et al. Liu et al. (2015), where individual rovibrational states of the dissociating molecules are binned together into a number of groups. While the multi-temperature and coarse-graining approaches have been used successfully in the past to reproduce the dissociation dynamics from detailed ab-initio calculations, there are still a few limitations. Namely, they rely on multiple coupled input parameters that must be evaluated consistently, and they also require the solution of additional species and/ or internal energy conservation equations.

Separately, the task of verifying dissociation rate constants across multiple different sources from the literature can be difficult. For hypersonic entry flows in particular, the dissociation rates for $\rm N_{2}$ / $\rm O_{2}$ and $\rm CO_{2}$ chemistry have been studied and reviewed in the past due to their relevance for Earth and Mars entry flows Park (1993); Park et al. (1994); Park, Jaffe, and Partridge (2001); Sarma (2000); Pietanza et al. (2021). Comparatively less work has been done to review and validate relevant data for $\rm H_{2}$ / He systems. $\rm H_{2}$ / He chemistry is of particular interest in the planetary entry community, as the recently published 2023-2032 decadal survey for planetary exploration identified probe missions to the giant plants Saturn and Uranus (whose atmospheres are composed primarily of $\rm H_{2}$ / He) as priorities National Academies of Sciences, Engineering, and Medicine (2022). As highlighted in the recent aerothermal heating uncertainty studies by Palmer et al. Palmer, Prabhu, and Cruden (2014) and Rataczak et al. Rataczak, Boyd, and McMahon (2024), uncertainties in $\rm H_{2}$ dissociation rate constants can have a significant impact on predicted convective and radiative heating rates. This presents a problem, as the accurate prediction of heating rates is crucial for the design of mass efficient thermal protection systems for entry probes. Many computational studies in the ice and gas giant entry flow literature Palmer, Prabhu, and Cruden (2014); Higdon et al. (2018); Liu et al. (2021); Hansson et al. (2021); Carroll and Brandis (2023); Carroll et al. (2023); Coelho and da Silva (2023); Rataczak, Boyd, and McMahon (2024) have used the dissociation rate constants proposed by Leibowitz Leibowitz (1973), despite the fact that the Leibowitz rates are only based on a single set of shock tube experiments (namely those of Jacobs et al. Jacobs, Giedt, and Cohen (1967)). Reviews of experimentally-reported $\rm H_{2}$ dissociation rate constants have been performed by Cohen and Westburg Cohen and Westberg (1983) and Baulch et al. Baulch (1972); Baulch et al. (2005). However, these reviews were primarily concerned with temperature ranges relevant for combustion applications (1,000 - 3,000 K), and did not consider the results of any computational studies. As will be shown later, computational studies are essential for constraining the gaps in the experimental data sets, especially at the high temperatures most relevant for hypersonic entry flow applications.

In light of these issues, the objective of the present work is twofold. First, to develop and validate a computationally inexpensive dissociation model that still captures the necessary physics from detailed ab-initio approaches. Second, to perform an extensive review of the corresponding model inputs for the case of $\rm H_{2}$ dissociation with the third-bodies $\rm H_{2}$ , H, and He/ inert gases. The paper is organized as follows. First, section II discusses the theory of non-equilibrium dissociation and the mathematical formulation of the master equations. From the master equation formulation, a simplified non-equilibrium rate constant expression (that is only a function of the translational temperature and the fraction of dissociation) is proposed. Next, in section III, the proposed expression is validated for the case of $\rm H_{2}$ dissociation against existing master equation results. Then, in section IV, a review of the available rate constant data for $\rm H_{2}$ dissociation from both experimental and computational sources is presented. Following this review, new rate constant fits that are valid across the temperature range from 200 - 20,000 K are proposed.

II Theory of Non-Equilibrium Dissociation

Non-equilibrium effects on dissociation are first described in terms of state-specific rates and the corresponding master equations in section II.1. Then, three important limits/ regimes of dissociation rates are discussed: the thermal equilibrium limit (II.2), the quasi-steady-state (QSS) regime (II.3), and finally the pre-QSS regime (II.4). Throughout this section, only the case of dissociation in a single-component mixture/ an infinitely dilute bath of a single third-body is considered. The case of dissociation in a multi-component mixture is discussed separately in Appendix A.

II.1 State-Specific Rates and Master Equations

For an arbitrary homonuclear diatomic molecule $\rm A_{2}$ , the dissociation/ recombination reaction through collisions with a third-body, M, can be written as follows

\rm A_{2}+M\leftrightarrow 2A+M.

(1)

The chemical source terms for this reaction are given by

\frac{dn_{\rm A_{2}}}{dt}=-n_{\rm A_{2}}n_{\rm M}k_{\rm d}+n_{\rm A}^{2}n_{\rm M}k_{\rm r}

(2)

and

\frac{dn_{\rm A}}{dt}=-2\frac{dn_{\rm A_{2}}}{dt}=2n_{\rm A_{2}}n_{\rm M}k_{\rm d}-2n_{\rm A}^{2}n_{\rm M}k_{\rm r},

(3)

where $t$ is time, $n_{\rm i}$ is the number density for species i, and $k_{\rm d}$ and $k_{\rm r}$ are the dissociation and recombination rate constants, respectively.

While reaction (1) and Eq. (2) and (3) are adequate for describing the dissociation/ recombination of $\rm A_{2}$ in an aggregate sense, they do not consider the detailed internal energy distributions of any of the reactants or products. Therefore, to capture the evolution of the rovibrational distribution of $\rm A_{2}$ , the same dissociation process of reaction (1) is written in terms of the reactions of the individual internal energy states of $\rm A_{2}$ ,

\rm A_{2}{\it(J,\nu)}+M\leftrightarrow 2A+M,

(4)

along with a separate reaction to describe the relaxation of the internal energy modes of $\rm A_{2}$ ,

\rm A_{2}{\it(J,\nu)}+M\leftrightarrow A_{2}{\it(J^{\prime},\nu^{\prime})}+M.

(5)

Here, $J$ and $\nu$ represent the rotational and vibrational quantum numbers, respectively. Throughout this work, only the electronic ground states of $\rm A_{2}$ and A are considered. The corresponding chemical source terms/ master equations for these state-specific reactions are given by

	$\displaystyle\frac{dn_{\rm A_{2}\it(J,\nu)}}{dt}=$		$\displaystyle n_{\rm M}\sum_{\nu^{\prime}=0}^{\nu_{\rm max}}\sum_{J^{\prime}=0}^{J_{\rm max}(\nu^{\prime})}n_{\rm A_{2}\it(J^{\prime},\nu^{\prime})}k(J^{\prime},\nu^{\prime}\rightarrow J,\nu)+n_{\rm M}n_{\rm A}^{2}k(c\rightarrow J,\nu)$		(6)
			$\displaystyle-n_{\rm M}\sum_{\nu^{\prime}=0}^{\nu_{\rm max}}\sum_{J^{\prime}=0}^{J_{\rm max}(\nu^{\prime})}n_{\rm A_{2}\it(J,\nu)}k(J,\nu\rightarrow J^{\prime},\nu^{\prime})-n_{\rm M}n_{\rm A_{2}\it(J,\nu)}k(J,\nu\rightarrow c)$		(6)

and

\frac{dn_{\rm A}}{dt}=2n_{\rm M}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}n_{\rm A_{2}\it(J,\nu)}k(J,\nu\rightarrow c)-2n_{\rm M}n_{\rm A}^{2}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}k(c\rightarrow J,\nu).

(7)

In these equations, the $(J,\nu)$ state-specific relaxation rate constants are denoted as $k(J,\nu\rightarrow J^{\prime},\nu^{\prime})$ and $k(J^{\prime},\nu^{\prime}\rightarrow J,\nu)$ , and the state-specific dissociation and recombination rate constants are denoted as $k(J,\nu\rightarrow c)$ and $k(c\rightarrow J,\nu)$ respectively, where $c$ indicates the dissociated state of $\rm A_{2}$ .

The forward and backward rate constants for the state-specific processes are related via micro-reversibility as

\frac{k(J,\nu\rightarrow c)}{k(c\rightarrow J,\nu)}=\frac{n_{\rm A,eq}^{2}}{n_{\rm A_{2}\it(J,\nu),\rm eq}}=\frac{Q_{\rm t,A}^{2}Q_{\rm A}^{2}}{Q_{\rm t,A_{2}}Q_{\rm A_{2}\it(J,\nu)}}\exp\left(-\frac{\theta_{\rm d,A_{2}}}{T_{\rm t}}\right)

(8)

and

\frac{k(J,\nu\rightarrow J^{\prime},\nu^{\prime})}{k(J^{\prime},\nu^{\prime}\rightarrow J,\nu)}=\frac{n_{\rm A_{2}\it(J^{\prime},\nu^{\prime}),\rm eq}}{n_{\rm A_{2}\it(J,\nu),\rm eq}}=\frac{Q_{\rm A_{2}\it(J^{\prime},\nu^{\prime})}}{Q_{\rm A_{2}\it(J,\nu)}}.

(9)

Here, the subscript eq refers to the locally defined equilibrium state evaluated at the translational temperature, $T_{\rm t}$ , and $\theta_{\rm d,A_{2}}$ is the characteristic dissociation temperature of $\rm A_{2}$ measured from the ground state. $Q_{\rm t,i}=(2\pi m_{\rm i}k_{\rm B}T_{\rm t}/h^{2})^{3/2}$ is the translational partition function with species’ mass $m_{i}$ and the Boltzmann and Planck constants, $k_{\rm B}$ and $h$ , respectively, while $Q_{\rm i}$ is the internal partition function for species i evaluated at $T_{\rm t}$ . The state-specific internal partition function, $Q_{\rm A_{2}\it(J,\nu)}$ , is related to the bulk $\rm A_{2}$ internal partition function, $Q_{\rm A_{2}}$ , as

Q_{\rm A_{2}}=\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}Q_{\rm A_{2}\it(J,\nu)}.

(10)

For the case of H and $\rm H_{2}$ (considering only electronic ground states), $Q_{\rm H}=2$ and $Q_{{\rm H_{2}}(J,\nu)}=g_{\rm s}(2J+1)\exp(-\theta_{\rm H_{2}\it(J,\nu)}/T_{\rm t})$ , where $\theta_{{\rm H_{2}}(J,\nu)}$ is the characteristic rovibrational temperature of ${\rm H_{2}}(J,\nu)$ , and the statistical weight due to nuclear spin-splitting, $g_{\rm s}$ , is 1/4 for all even $J$ (para-hydrogen) and 3/4 for all odd $J$ (ortho-hydrogen) Colonna, D’Angola, and Capitelli (2012); Popovas and Jørgensen (2016).

The state-specific rate constants, $k(J,\nu\rightarrow c)$ , $k(c\rightarrow J,\nu)$ , $k(J,\nu\rightarrow J^{\prime},\nu^{\prime})$ , and $k(J^{\prime},\nu^{\prime}\rightarrow J,\nu)$ , are all only functions of the translational temperature, $T_{\rm t}$ . This is a direct consequence of the fact that the reaction of a given $\rm A_{2}\it(J,\nu)$ state is only a function of its impact with the third-body, M. Therefore, if the internal energy of M does not change during a collision, then the state-specific rate constants can only be a function of $T_{\rm t}$ . In the event that M is another molecule that could simultaneously undergo a relaxation or reaction, the state-specific rate constants can still be written as just functions of $T_{\rm t}$ , provided that the rate constants are written to index the internal energy states of both $\rm A_{2}\it(J,\nu)$ and M $(J,\nu)$ Kim and Boyd (2012). For the sake of simplicity, the less general case of M with a fixed internal energy will be considered for the analysis presented here.

Comparing the atomic master equation of Eq. (7) to Eq. (3), it is evident that the aggregate rate constants can be expressed in terms of the state-specific rate constants as

k_{\rm d}=\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\frac{n_{\rm A_{2}\it(J,\nu)}}{n_{\rm A_{2}}}k(J,\nu\rightarrow c)

(11)

and

k_{\rm r}=\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}k(c\rightarrow J,\nu),

(12)

where

n_{\rm A_{2}}=\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}n_{\rm A_{2}\it(J,\nu)}.

(13)

Using the micro-reversibility relation of Eq. (8), Eq. (12) can equivalently be written as

k_{\rm r}=K_{\rm eq}^{-1}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\frac{Q_{\rm A_{2}\it(J,\nu)}}{Q_{\rm A_{2}}}k(J,\nu\rightarrow c).

(14)

Here, we have used the fact that $n_{\rm A_{2},eq}$ and $n_{\rm A_{2}{\it(J,\nu)},eq}$ are related via their internal partition functions as

\frac{n_{\rm A_{2}\it(J,\nu)\rm,eq}}{n_{\rm A_{2},eq}}=\frac{Q_{\rm A_{2}\it(J,\nu)}}{Q_{\rm A_{2}}},

(15)

and that $n_{\rm A_{2},eq}$ and $n_{\rm A,eq}$ are related via the equilibrium constant, $K_{\rm eq}$ , as

K_{\rm eq}=\frac{n_{\rm A,eq}^{2}}{n_{\rm A_{2},eq}}=\frac{Q_{\rm t,A}^{2}Q_{\rm A}^{2}}{Q_{\rm t,A_{2}}Q_{\rm A_{2}}}\exp\left(-\frac{\theta_{\rm d,A_{2}}}{T_{\rm t}}\right).

(16)

Because Eq. (11) and (14) come from the master equations of Eq. (6) and (7), these results are general, and hold for any degree of thermochemical non-equilibrium.

In a compression/ shocked flow environment, the rovibrational distribution of $\rm A_{2}$ initially begins as a Boltzmann distribution at the pre-shock temperature, which then evolves towards a Boltzmann distribution at the post-shock equilibrium temperature. In the general case, this evolution can occur through a series of non-Boltzmann distributions, in contrast to the assumption of a typical multi-temperature model, which assumes that it must occur through a series of Boltzmann distributions. Because the distribution of $\rm A_{2}\it(J,\nu)$ can continually change in this post-shock environment, Eq. (11) gives that the effective aggregate $k_{\rm d}$ can continually change as well. Therefore, unlike the state-specific rate constants, $k_{\rm d}$ is a function of the instantaneous distribution of $\rm A_{2}\it(J,\nu)$ .

In contrast, Eq. (12) gives that $k_{\rm r}$ must always be a function of the state-specific rates alone, and therefore, can only ever be a function of $T_{\rm t}$ . This will be of particular importance when interpreting the results of the QSS analysis presented in section II.3.

II.2 Thermal Rates

In the limit of thermal equilibrium, the rovibrational distribution of $\rm A_{2}$ can be described simply by Eq. (15). Therefore, substituting Eq. (15) into Eq. (11) gives the thermal dissociation rate constant,

k_{\rm d,th}=\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\frac{Q_{\rm A_{2}\it(J,\nu)}}{Q_{\rm A_{2}}}k(J,\nu\rightarrow c).

(17)

This thermal rate is referred to as the “one-way” rate by Park Park (1989b) and Kim and Boyd Kim and Boyd (2012, 2013). Taking the ratio of Eq. (17) over the general expression for $k_{\rm r}$ given by Eq. (14), the expected result that macroscopic detailed balance must be satisfied at thermal equilibrium is obtained,

\frac{k_{\rm d,th}}{k_{\rm r}}=K_{\rm eq}.

(18)

Since Eq. (14) is a general relation for $k_{\rm r}$ , this result implies that regardless of the degree of non-equilibrium, $k_{\rm r}$ must always satisfy Eq. (18).

II.3 QSS Rates

II.3.1 Formulation and Solution for $\bm{A_{2}(J,\nu)}$

A second regime of importance for the distribution of $\rm A_{2}\it(J,\nu)$ and therefore $k_{\rm d}$ is the quasi-steady-state or QSS condition. The QSS condition arises when the rate of change of $n_{\rm A_{2}{\it(J,\nu)}}$ is significantly smaller than the total production term in the right-hand side of Eq. (6),

\left|\frac{dn_{\rm A_{2}\it(J,\nu)}}{dt}\right|\ll n_{\rm M}\sum_{\nu^{\prime}=0}^{\nu_{\rm max}}\sum_{J^{\prime}=0}^{J_{\rm max}(\nu^{\prime})}n_{\rm A_{2}\it(J^{\prime},\nu^{\prime})}k(J^{\prime},\nu^{\prime}\rightarrow J,\nu)+n_{\rm M}n_{\rm A}^{2}k(c\rightarrow J,\nu),

(19)

or equivalently the total consumption term,

\left|\frac{dn_{\rm A_{2}\it(J,\nu)}}{dt}\right|\ll n_{\rm M}\sum_{\nu^{\prime}=0}^{\nu_{\rm max}}\sum_{J^{\prime}=0}^{J_{\rm max}(\nu^{\prime})}n_{\rm A_{2}\it(J,\nu)}k(J,\nu\rightarrow J^{\prime},\nu^{\prime})+n_{\rm M}n_{\rm A_{2}\it(J,\nu)}k(J,\nu\rightarrow c).

(20)

When this condition is satisfied, the total consumption term is approximately equal to the total production term such that

	$\displaystyle\sum_{\nu^{\prime}=0}^{\nu_{\rm max}}\sum_{J^{\prime}=0}^{J_{\rm max}(\nu^{\prime})}n_{\rm A_{2}\it(J,\nu)}k(J,\nu\rightarrow J^{\prime},\nu^{\prime})+n_{\rm A_{2}\it(J,\nu)}k(J,\nu\rightarrow c)$
	$\displaystyle\approx\sum_{\nu^{\prime}=0}^{\nu_{\rm max}}\sum_{J^{\prime}=0}^{J_{\rm max}(\nu^{\prime})}n_{\rm A_{2}\it(J^{\prime},\nu^{\prime})}k(J^{\prime},\nu^{\prime}\rightarrow J,\nu)+n_{\rm A}^{2}k(c\rightarrow J,\nu).$		(21)

To obtain the corresponding dissociation rate constant at QSS, a similar setup to the one proposed by Park Park (1989b) is used here. In particular, the local equilibrium normalized number densities,

\phi_{\rm A_{2}\it(J,\nu)}\equiv\frac{n_{\rm A_{2}\it(J,\nu)}}{n_{\rm A_{2}\it(J,\nu),\rm eq}},

(22)

\phi_{\rm A}\equiv\frac{n_{\rm A}}{n_{\rm A,eq}},

(23)

and

\chi\equiv\frac{n_{\rm A_{2}}}{n_{\rm A_{2},eq}}

(24)

are defined. For convenience, the equilibrium normalized fractional density for ${\rm A_{2}}(J,\nu)$ is also defined as

\psi_{\rm A_{2}\it(J,\nu)}\equiv\frac{\phi_{\rm A_{2}\it(J,\nu)}}{\chi}.

(25)

Then, Eq. (21) is rewritten using these definitions and Eq. (8), (9), and (15) to get

	$\displaystyle\phi_{\rm A_{2}\it(J,\nu)}\left[\sum_{\nu^{\prime}=0}^{\nu_{\rm max}}\sum_{J^{\prime}=0}^{J_{\rm max}(\nu^{\prime})}k(J,\nu\rightarrow J^{\prime},\nu^{\prime})+k(J,\nu\rightarrow c)\right]-\sum_{\nu^{\prime}=0}^{\nu_{\rm max}}\sum_{J^{\prime}=0}^{J_{\rm max}(\nu^{\prime})}\phi_{\rm A_{2}\it(J^{\prime},\nu^{\prime})}k(J,\nu\rightarrow J^{\prime},\nu^{\prime})$
	$\displaystyle=\phi_{\rm A}^{2}k(J,\nu\rightarrow c).$		(26)

For $N$ total rovibrational states of $\rm A_{2}$ , Eq. (26) gives $N$ constraint equations. However, Eq. (13) gives an additional constraint for the total number density of $\rm A_{2}$ that must be satisfied as well. Written in terms of $\phi_{\rm A_{2}\it(J,\nu)}$ , this additional constraint is expressed as

\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\phi_{\rm A_{2}\it(J,\nu)}n_{\rm A_{2}\it(J,\nu)\rm,eq}=\chi n_{\rm A_{2},eq},

(27)

or equivalently

\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\phi_{\rm A_{2}\it(J,\nu)}\frac{Q_{\rm A_{2}\it(J,\nu)}}{Q_{\rm A_{2}}}=\chi.

(28)

This gives an over-constrained system with $N+1$ equations for $N$ unknowns. To resolve this over-constraint, the QSS condition given by Eq. (26) for the ground state, $\rm A_{2}$ ( $J$ =0, $\nu$ =0), is dropped. This follows from the understanding that the QSS condition arises as a competition between the production process of excitation and the consumption process of dissociation, as the magnitudes of the excitation and dissociation terms in Eq. (26) are initially much larger than the magnitudes of the de-excitation and recombination terms in compression/ shocked flows. As a consequence, the QSS condition is least likely to be satisfied by the ground state, as unlike all other excited rovibrational states of $\rm A_{2}\it(J,\nu)$ , there is no excitation term that can repopulate the ground state and compete with the dissociation process.

Equations (26) and (28) can then be written in the matrix form

\bm{M}\vec{\phi}_{\rm A_{2}}=\vec{C}\phi_{\rm A}^{2}+\vec{D}\chi.

(29)

The entries in the first row of matrix $\bm{M}$ are given by $Q_{\rm A_{2}\it(J,\nu)}/Q_{\rm A_{2}}$ , and the first entries of the column vectors $\vec{C}$ and $\vec{D}$ are 0 and 1, respectively. For all other rows, the entries of $\bm{M}$ are given by the coefficients in the left-hand side of Eq. (26), the entries of $\vec{C}$ are given by the coefficients in the right-hand side of Eq. (26), and $\vec{D}=0$ . The solution to Eq. (29) can be written as the sum of two terms,

\vec{\psi}_{\rm A_{2},nr}\equiv\bm{M}^{-1}\vec{D}

(30)

and

\vec{\psi}_{\rm A_{2},r}\equiv\bm{M}^{-1}\vec{C},

(31)

such that

\vec{\phi}_{\rm A_{2}}=\vec{\psi}_{\rm A_{2},nr}\chi+\vec{\psi}_{\rm A_{2},r}\phi_{\rm A}^{2}.

(32)

$\vec{\psi}_{\rm A_{2},nr}$ represents the QSS solution in the absence of recombination/ in the non-recombining (nr) limit, and $\vec{\psi}_{\rm A_{2},r}$ represents the contribution of recombination (r) to the QSS solution. In Park Park (1989b) as well as the subsequent master equation studies by Kim et al. Kim, Kwon, and Park (2009, 2010) and Kim and Boyd Kim and Boyd (2012, 2013), the solution to Eq. (29) was instead described as the sum of a homogeneous component, $\vec{\psi}_{\rm A_{2},h}\equiv{\bm{M}}^{-1}\vec{D}\chi$ , and a particular component, $\vec{\psi}_{\rm A_{2},p}\equiv{\bm{M}}^{-1}\vec{C}$ . The formulation of the present work differs in that $\chi$ is pulled out of the definition of the homogeneous solution, such that $\vec{\psi}_{\rm A_{2},nr}$ and $\vec{\psi}_{\rm A_{2},r}$ both represent particular solutions to Eq. (29) that are only functions of ${\bm{M}}$ , $\vec{D}$ , and $\vec{C}$ .

Two key observations are made here regarding the QSS formulation described by Eq. (30)-(32). Firstly, $\bm{M}$ , $\vec{C}$ , and $\vec{D}$ are only functions of the state-specific rates, which in turn are only functions of $T_{\rm t}$ . As a result, $\vec{\psi}_{\rm A_{2},nr}$ and $\vec{\psi}_{\rm A_{2},r}$ are also only functions of $T_{\rm t}$ . Therefore, even though the QSS distribution may correspond to internal rovibrational temperatures, $T_{\rm r}$ and $T_{\rm v}$ , that are not equal to $T_{\rm t}$ , the QSS solution is none-the-less functionalized by $T_{\rm t}$ alone. Secondly, since the QSS assumption of approximately equal production and consumption terms as given by Eq. (21) is true at equilibrium, Eq. (32) must also be valid at equilibrium. By construction, $\phi_{\rm A}$ and $\chi$ are equal to 1 and $\vec{\phi}_{\rm A_{2}}$ is the unity vector at equilibrium. Hence, Eq. (32) implies that

1=\vec{\psi}_{\rm A_{2},nr}+\vec{\psi}_{\rm A_{2},r}.

(33)

II.3.2 Solution for $\bm{k_{d}}$

To compute the dissociation rate constant at QSS, the general expression for $k_{\rm d}$ given by Eq. (11) is first rewritten in terms of $\psi_{{\rm A_{2}}(J,\nu)}$ using Eq. (15), (22), (24), and (25) as

k_{\rm d}=\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\psi_{\rm A_{2}\it(J,\nu)}\frac{Q_{\rm A_{2}{\it(J,\nu)}}}{Q_{\rm A_{2}}}k(J,\nu\rightarrow c).

(34)

Separately, Eq. (25), (32), and (33) are combined to give

\vec{\psi}_{\rm A_{2}}=\vec{\psi}_{\rm A_{2},nr}\left(1-\frac{\phi_{\rm A}^{2}}{\chi}\right)+\frac{\phi_{\rm A}^{2}}{\chi}.

(35)

Equation (35) is then substituted into Eq. (34) to give the final expression for $k_{\rm d}$ at QSS,

k_{\rm d}=\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\left(\psi_{\rm A_{2}{\it(J,\nu)},nr}\left(1-\frac{\phi_{\rm A}^{2}}{\chi}\right)+\frac{\phi_{\rm A}^{2}}{\chi}\right)\frac{Q_{\rm A_{2}\it(J,\nu)}}{Q_{\rm A_{2}}}k(J,\nu\rightarrow c).

(36)

Since $\vec{\psi}_{\rm A_{2},nr}$ represents the QSS solution in the absence of recombination, the corresponding dissociation rate constant in the non-recombining limit, $k_{\rm d,nr}$ , can be computed simply from Eq. (34) as

k_{\rm d,nr}=\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\psi_{\rm A_{2}\it(J,\nu),\rm nr}\frac{Q_{\rm A_{2}\it(J,\nu)}}{Q_{\rm A_{2}}}k(J,\nu\rightarrow c).

(37)

Using Eq. (37) in addition to the expression for the thermal rate constant, Eq. (17), Eq. (36) can be written more compactly as

k_{\rm d}=k_{\rm d,nr}\left(1-\frac{\phi_{\rm A}^{2}}{\chi}\right)+k_{\rm d,th}\frac{\phi_{\rm A}^{2}}{\chi}.

(38)

Here, the coefficient

\frac{\phi_{\rm A}^{2}}{\chi}=\frac{n_{\rm A}^{2}}{n_{\rm A_{2}}}\frac{n_{\rm A_{2},eq}}{n_{\rm A,eq}^{2}}=\frac{n_{\rm A}^{2}}{n_{\rm A_{2}}}\frac{1}{K_{\rm eq}}

(39)

is effectively a progress variable that varies from 0 in the limit of no dissociation to 1 in the limit of equilibrium. As such, Eq. (38) gives a description of $k_{\rm d}$ that varies between the two limiting rates of $k_{\rm d,nr}$ and $k_{\rm d,th}$ . Since $k_{\rm d,nr}$ and $k_{\rm d,th}$ are only functions of $T_{\rm t}$ , the description of $k_{\rm d}$ given by Eq. (38) is only a function of $T_{\rm t}$ and $\phi_{\rm A}^{2}/\chi$ . Despite this lack of dependence on any other internal temperatures (e.g., $T_{\rm r}$ or $T_{\rm v}$ ), the $k_{\rm d}$ given by Eq. (38) only satisfies macroscopic detailed balance with $k_{\rm r}$ exactly in the limit where $\phi_{\rm A}^{2}/\chi=1$ and hence $k_{\rm d}=k_{\rm d,th}$ . This implies that at this limit of the QSS condition, thermal and chemical equilibria are reached simultaneously.

Finally, it is worth noting that the description of $k_{\rm d}$ given by Eq. (38) is valid everywhere where the QSS assumption of Eq. (21) is valid, namely everywhere outside of the pre-QSS region. Some master equation studies in the literature, such as the study by Venturi et al. Venturi et al. (2020) for the dissociation of $\rm O_{2}$ with M = O, distinguish between a “QSS” and “post-QSS” region. The discussion above suggests that it would be more appropriate to classify these regions as the non-recombining (nr) and recombining (r) regions, where $k_{\rm d}\approx k_{\rm d,nr}$ or $k_{\rm d,nr}<k_{\rm d}\leq k_{\rm d,th}$ , respectively.

II.3.3 Chemical Source Term

While Eq. (38) provides an expression for $k_{\rm d}$ that is valid for all regimes where the QSS assumption of Eq. (21) is valid, it implies that separate fits for both $k_{\rm d,nr}$ and $k_{\rm d,th}$ in terms of $T_{\rm t}$ are required. However, if we are only concerned with reproducing the total chemical source term for each species, i.e., $dn_{\rm A_{2}}/dt$ and $dn_{\rm A}/dt$ , as is typically the case in CFD calculations, a further simplification can be made that eliminates the dependence on $k_{\rm d,th}$ . In particular, if Eq. (38) is substituted into the general aggregate source term equation, Eq. (3), $dn_{\rm A}/dt$ can be written as

\frac{dn_{\rm A}}{dt}=2n_{\rm A_{2}}n_{\rm M}\left(k_{\rm d,nr}\left(1-\frac{\phi_{\rm A}^{2}}{\chi}\right)+k_{\rm d,th}\frac{\phi_{\rm A}^{2}}{\chi}\right)-2n_{\rm A}^{2}n_{\rm M}k_{\rm r}.

(40)

If Eq. (39) is then used to substitute for $\phi_{\rm A}^{2}/\chi$ , Eq. (40) can be written as

\begin{split}\frac{dn_{\rm A}}{dt}&=2n_{\rm A_{2}}n_{\rm M}k_{\rm d,nr}-2n_{\rm A}^{2}n_{\rm M}\frac{k_{\rm d,nr}}{K_{\rm eq}}.\end{split}

(41)

Therefore, with no additional simplifying assumptions, the chemical source term from the full formulation of Eq. (40) can be computed equivalently as solely a function of $k_{\rm d,nr}$ . What’s more, the dependence on $k_{\rm d,nr}$ in Eq. (41) is exactly the same as it would be if one had “incorrectly” assumed that $k_{\rm d}=k_{\rm d,nr}$ and applied macroscopic detailed balance to compute the consumption term. However, the above analysis clearly shows that the consumption term with $k_{\rm d,nr}/K_{\rm eq}$ does not correspond to the physical recombination rate constant, $k_{\rm r}$ , but instead comes from the negative term of the dissociation rate constant given by Eq. (38).

While not explicitly stated or justified, the alternate formulation given by Eq. (41) is inherently what Park Park (1989b) and Kim Kim, Kwon, and Park (2009, 2010); Kim and Boyd (2012, 2013) have used. In their formulation, the recombination rate constant was simply defined as

k_{\rm r,Park}=\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}(1-\psi_{\rm A_{2}\it(J,\nu),\rm r})\frac{n_{\rm A_{2},eq}}{n_{\rm A,eq}^{2}}\frac{Q_{\rm A_{2}\it(J,\nu)}}{Q_{\rm A_{2}}}k(J,\nu\rightarrow c).

(42)

Using Eq. (16), (33), and (37), one can show that the above expression is equivalent to

k_{\rm r,Park}=\frac{k_{\rm d,nr}}{K_{\rm eq}}.

(43)

Therefore, $k_{\rm r,Park}$ is defined such that it satisfies macroscopic detailed balance with $k_{\rm d,nr}$ by construction.

There are two key implications of these results. First, the same $k_{\rm d,nr}$ expression should be used to compute both the production and consumption terms of Eq. (41). Second, since $k_{\rm d,nr}$ is only a function of $T_{\rm t}$ in the QSS regime, one significant reason to use a fit of $k_{\rm d,nr}$ that depends on other variables is to capture dissociation in regions where the QSS assumption is not valid, i.e., in the pre-QSS regime.

II.3.4 Transition Between $\bm{k}_{\bm{d},nr}$ and $\bm{k}_{\bm{d},th}$

Leveraging the above analysis, it is possible to determine explicitly the degree of dissociation that is reached when $k_{\rm d}$ exceeds the $k_{\rm d,nr}$ limit. Mathematically, this limit is defined as the degree of dissociation reached when $k_{\rm d}$ exceeds some $\delta$ fraction of $k_{\rm d,nr}$ , i.e., $k_{\rm d}=k_{\rm d,nr}(1+\delta)$ . Substituting this limit for $k_{\rm d}$ into Eq. (38) gives the limiting value for $\phi_{\rm A}^{2}/\chi$ as

\frac{\phi_{\rm A}^{2}}{\chi}\biggr{|}_{\rm nr}=\frac{\delta k_{\rm d,nr}}{k_{\rm d,th}-k_{\rm d,nr}}.

(44)

Next, the degree of dissociation, $\alpha$ , is defined as

\alpha\equiv\frac{n_{\rm A}}{\tilde{n}_{\rm A}},

(45)

where $\tilde{n}_{\rm A}\equiv 2n_{\rm A_{2}}+n_{\rm A}$ is the total number density of A atoms. From this, the number densities of $\rm A_{2}$ and A can be written as

n_{\rm A_{2}}=\frac{1-\alpha}{2}\tilde{n}_{\rm A},\quad n_{\rm A}=\alpha\tilde{n}_{\rm A}\,.

(46)

By construction, $\alpha$ is related to $\phi_{\rm A}$ as

\frac{\alpha}{\alpha_{\rm eq}}=\frac{n_{\rm A}}{\tilde{n}_{\rm A}}\frac{\tilde{n}_{\rm A}}{n_{\rm A,eq}}=\phi_{\rm A},

(47)

such that $\phi_{\rm A}$ is a direct measure of the fraction of dissociation relative to the locally defined equilibrium at $T_{\rm t}$ .

Using Eq. (46) along with Eq. (39), $\phi_{\rm A}^{2}/\chi$ takes the form

\frac{\phi_{\rm A}^{2}}{\chi}=\frac{2\alpha^{2}\tilde{n}_{\rm A}}{(1-\alpha)K_{\rm eq}},

(48)

with the corresponding solution for $\alpha$ given by

\alpha=\frac{1}{2}\left[-\frac{K_{\rm eq}}{2\tilde{n}_{\rm A}}\frac{\phi_{\rm A}^{2}}{\chi}+\sqrt{\frac{K_{\rm eq}}{2\tilde{n}_{\rm A}}\frac{\phi_{\rm A}^{2}}{\chi}\left(\frac{K_{\rm eq}}{2\tilde{n}_{\rm A}}\frac{\phi_{\rm A}^{2}}{\chi}+4\right)}\right].

(49)

To account for the fact that the equilibrium degree of dissociation can vary as a function of $T_{\rm t}$ and $\tilde{n}_{\rm A}$ , $\alpha$ is normalized by $\alpha_{\rm eq}$ (computed using Eq. (49) with $\phi_{\rm A}^{2}/\chi$ = 1) and the value of $\phi_{A}^{2}/\chi$ in the non-recombining limit from Eq. (44) is substituted to give

\phi_{\rm A}\big{|}_{\rm nr}=\frac{\alpha}{\alpha_{\rm eq}}\biggr{|}_{\rm nr}=\frac{\phi_{\rm A}^{2}}{\chi}\biggr{|}_{\rm nr}\frac{\sqrt{1+8\tilde{n}_{\rm A}/K_{\rm eq}/(\phi_{\rm A}^{2}/\chi)|_{\rm nr}}-1}{\sqrt{1+8\tilde{n}_{\rm A}/K_{\rm eq}}-1}.

(50)

Practically, to use this analysis, one first selects a threshold value of $\delta$ to define the $k_{\rm d,nr}$ limit. Then, the corresponding ratio $(\phi_{A}^{2}/\chi)|_{\rm nr}$ is evaluated using Eq. (44). Finally, Eq. (50) is used to give the fraction of dissociation that occurs in the non-recombining limit. Importantly, these results are only a function of $T_{\rm t}$ and $\tilde{n}_{\rm A}$ , and hence it is not necessary to solve the full system of master equations to know how these variables will impact the transition of $k_{\rm d}$ from the non-recombining to the thermal limit.

II.3.5 Solution for $\bm{e}_{\bm{r}v}$

Since Eq. (35) gives a complete description of the ${\rm A_{2}}(J,\nu)$ distribution in the QSS regime, the corresponding evolution of the average rovibrational energy of $\rm A_{2}$ can be computed explicitly as well. In general, the average rovibrational energy of $\rm A_{2}$ is given by

e_{\rm rv}=k_{\rm B}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\frac{n_{\rm A_{2}\it(J,\nu)}}{n_{\rm A_{2}}}\theta_{\rm A_{2}\it(J,\nu)}=k_{\rm B}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\psi_{\rm A_{2}{\it(J,\nu)}}\frac{Q_{\rm A_{2}\it(J,\nu)}}{Q_{\rm A_{2}}}\theta_{\rm A_{2}\it(J,\nu)},

(51)

where $\theta_{{\rm A_{2}}(J,\nu)}$ is the characteristic rovibrational temperature of ${\rm A_{2}}(J,\nu)$ . Using Eq. (35), $e_{\rm rv}$ in the QSS regime can then be expressed as

e_{\rm rv}=e_{\rm rv,nr}\left(1-\frac{\phi_{\rm A}^{2}}{\chi}\right)+e_{\rm rv,th}\frac{\phi_{\rm A}^{2}}{\chi},

(52)

where $e_{\rm rv}$ in the non-recombining and thermal limits has been defined as

e_{\rm rv,nr}\equiv k_{\rm B}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\psi_{\rm A_{2}{\it(J,\nu)},nr}\frac{Q_{\rm A_{2}\it(J,\nu)}}{Q_{\rm A_{2}}}\theta_{\rm A_{2}\it(J,\nu)}

(53)

and

e_{\rm rv,th}\equiv k_{\rm B}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\frac{Q_{\rm A_{2}\it(J,\nu)}}{Q_{\rm A_{2}}}\theta_{\rm A_{2}\it(J,\nu)},

(54)

respectively. As with $k_{\rm d,nr}$ and $k_{\rm d,th}$ , Eq. (53) and (54) give that $e_{\rm rv,nr}$ and $e_{\rm rv,th}$ are only functions of $T_{\rm t}$ , and hence $e_{\rm rv}$ in the QSS regime is functionalized by $T_{\rm t}$ and $\phi_{\rm A}^{2}/\chi$ alone.

Unlike $k_{\rm d,nr}$ , $e_{\rm rv,nr}$ is not a commonly reported quantity in the literature. Fortunately, it is still possible to obtain estimates of $e_{\rm rv,nr}$ if energy-equivalent rotational and vibrational temperatures, $T_{\rm r}$ and $T_{\rm v}$ , are reported. In a vibration-preferred framework where $\theta_{\rm r,A_{2}\it(J,\nu)}\equiv\theta_{\rm A_{2}\it(J,\nu)}-\theta_{\rm A_{2}\it(J={\rm 0},\nu)}$ and $\theta_{\rm v,A_{2}\it(J,\nu)}\equiv\theta_{\rm A_{2}\it(J={\rm 0},\nu)}$ , $T_{\rm r}$ and $T_{\rm v}$ are defined such that the average rotational and vibrational energies, $e_{\rm r}$ and $e_{\rm v}$ , are given by

e_{\rm r}\equiv k_{\rm B}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\frac{n_{\rm A_{2}\it(J,\nu)}}{n_{\rm A_{2}}}\theta_{\rm r,A_{2}\it(J,\nu)}\equiv k_{\rm B}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\frac{Q_{\rm A_{2}\it(J,\nu)}(T_{\rm r},T_{\rm v})}{Q_{\rm A_{2}}(T_{\rm r},T_{\rm v})}\theta_{\rm r,A_{2}\it(J,\nu)}

(55)

and

e_{\rm v}\equiv k_{\rm B}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\frac{n_{\rm A_{2}\it(J,\nu)}}{n_{\rm A_{2}}}\theta_{\rm v,A_{2}\it(J,\nu)}\equiv k_{\rm B}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\frac{Q_{\rm A_{2}\it(J,\nu)}(T_{\rm r},T_{\rm v})}{Q_{\rm A_{2}}(T_{\rm r},T_{\rm v})}\theta_{\rm v,A_{2}\it(J,\nu)},

(56)

where

e_{\rm rv}=e_{\rm r}+e_{\rm v}

(57)

by construction. Therefore, as long as estimates of $T_{\rm r,nr}(T_{\rm t})$ and $T_{\rm v,nr}(T_{\rm t})$ can be obtained, $e_{\rm rv,nr}$ can be computed equivalently as

e_{\rm rv,nr}=k_{\rm B}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\frac{Q_{\rm A_{2}\it(J,\nu)}(T_{\rm r,nr},T_{\rm v,nr})}{Q_{\rm A_{2}}(T_{\rm r,nr},T_{\rm v,nr})}\theta_{\rm rv,A_{2}\it(J,\nu)}.

(58)

II.4 Pre-QSS Rates

II.4.1 Formulation

The theory of the previous sections provides a complete description of $k_{\rm d}$ and $e_{\rm rv}$ under the QSS assumption up to and including thermal equilibrium. However, this description is of course only valid once the QSS condition of Eq. (21) is met. The region before QSS is established is referred to as the “pre-QSS” region. If it is assumed that recombination is negligible in the pre-QSS region as it occurs before the non-recombining QSS limit, the master equations for ${\rm A_{2}}(J,\nu)$ from Eq. (6) can be simplified to

	$\displaystyle\frac{dn_{\rm A_{2}\it(J,\nu)}}{dt}=$		$\displaystyle n_{\rm M}\sum_{\nu^{\prime}=0}^{\nu_{\rm max}}\sum_{J^{\prime}=0}^{J_{\rm max}(\nu^{\prime})}n_{\rm A_{2}\it(J^{\prime},\nu^{\prime})}k(J^{\prime},\nu^{\prime}\rightarrow J,\nu)-n_{\rm M}\sum_{\nu^{\prime}=0}^{\nu_{\rm max}}\sum_{J^{\prime}=0}^{J_{\rm max}(\nu^{\prime})}n_{\rm A_{2}\it(J,\nu)}k(J,\nu\rightarrow J^{\prime},\nu^{\prime})$		(59)
			$\displaystyle-n_{\rm M}n_{\rm A_{2}\it(J,\nu)}k(J,\nu\rightarrow c).$		(59)

The inherent time-dependence and coupled nature of Eq. (59) implies that it is only possible to obtain an exact solution in the pre-QSS by solving Eq. (59) for all states of ${\rm A_{2}}(J,\nu)$ simultaneously. However, under several simplifying assumptions, useful insights into the pre-QSS behavior can still be obtained.

The isothermal and isochoric assumptions are first made such that all equilibrium number densities can be treated as constants. Under this assumption, Eq. (59) can be rewritten in terms of $\phi_{\rm A_{2}(J,\nu)}$ as

	$\displaystyle\frac{d\phi_{\rm A_{2}\it(J,\nu)}}{dt}=$		$\displaystyle n_{\rm M}\sum_{\nu^{\prime}=0}^{\nu_{\rm max}}\sum_{J^{\prime}=0}^{J_{\rm max}(\nu^{\prime})}\phi_{\rm A_{2}\it(J^{\prime},\nu^{\prime})}k(J,\nu\rightarrow J^{\prime},\nu^{\prime})$		(60)
			$\displaystyle-n_{\rm M}\phi_{\rm A_{2}\it(J,\nu)}\left[\sum_{\nu^{\prime}=0}^{\nu_{\rm max}}\sum_{J^{\prime}=0}^{J_{\rm max}(\nu^{\prime})}k(J,\nu\rightarrow J^{\prime},\nu^{\prime})+k(J,\nu\rightarrow c)\right],$		(60)

or equivalently in matrix form as

\frac{d\vec{\phi}_{\rm A_{2}}}{dt}=-n_{\rm M}\bm{\tilde{M}}\vec{\phi}_{\rm A_{2}},

(61)

where the entries of $\bm{\tilde{M}}$ correspond to the state-specific rate constant terms in the right-hand side of Eq. (60). This definition of $\bm{\tilde{M}}$ is nearly identical to the definition of $\bm{M}$ used previously in the QSS formulation of Eq. (29), except that the first row of $\bm{\tilde{M}}$ does not correspond to the constraint equation of Eq. (28) as it does for $\bm{M}$ . Next, the rate of change of $\chi$ is expressed in terms of $\vec{\phi}_{\rm A_{2}}$ by taking the time derivative of Eq. (28) and applying Eq. (61) to get

\frac{d\chi}{dt}=\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\frac{d\phi_{\rm A_{2}\it(J,\nu)}}{dt}\frac{Q_{\rm A_{2}\it(J,\nu)}}{Q_{\rm A_{2}}}=-{n_{\rm M}}^{t}\vec{Q}_{\rm A_{2}}\bm{\tilde{M}}\vec{\phi}_{\rm A_{2}}.

(62)

Here, ${}^{t}\vec{Q}_{\rm A_{2}}$ corresponds to the transpose of the vector with entries given by $Q_{\rm A_{2}\it(J,\nu)}/Q_{\rm A_{2}}$ , and Eq. (28) gives that ${}^{t}\vec{Q}_{\rm A_{2}}\vec{\psi}_{\rm A_{2}}=1$ . These results can be combined to express the time rate of change of $\vec{\psi}_{\rm A_{2}}$ as

\frac{d\vec{\psi}_{\rm A_{2}}}{dt}=\chi^{-1}\frac{d\vec{\phi}_{\rm A_{2}}}{dt}+\vec{\phi}_{\rm A_{2}}\frac{d\chi^{-1}}{dt}=-n_{\rm M}\left[\tilde{M}\vec{\psi}_{\rm A_{2}}-\left({}^{t}\vec{Q}_{\rm A_{2}}\bm{\tilde{M}}\vec{\psi}_{\rm A_{2}}\right)\vec{\psi}_{\rm A_{2}}\right].

(63)

It can be shown that all eigenvalues of $\bm{\tilde{M}}$ are negative, and hence as $t\rightarrow\infty$ , the solution to Eq. (63) converges to a steady state given by

\tilde{M}\vec{\psi}_{\rm A_{2},\infty}=\left({}^{t}\vec{Q}_{\rm A_{2}}\bm{\tilde{M}}\vec{\psi}_{\rm A_{2},\infty}\right)\vec{\psi}_{\rm A_{2},\infty}.

(64)

In other words, $\vec{\psi}_{\rm A_{2},\infty}$ is an eigenvector of $\bm{\tilde{M}}$ with the corresponding eigenvalue

\lambda_{0}=\ ^{t}\vec{Q}_{\rm A_{2}}\bm{\tilde{M}}\vec{\psi}_{\rm A_{2},\infty}.

(65)

A linear stability analysis of Eq. (61) would show that $\lambda_{0}$ is the smallest eigenvalue and hence corresponds to the slowest decaying mode of $\bm{\tilde{M}}$ .

II.4.2 Approximate Solution for $\bm{A}_{2}(J,\nu)$ and $\bm{k}_{\bm{d}}$

Until this point in the analysis of the pre-QSS, only the isothermal and isochoric assumptions have been made. To compute a solution for the evolution of $A_{2}(J,\nu)$ without numerically integrating the master equations, three additional approximations are made here. Firstly, the nonlinear term of Eq. (63) is approximated using Eq. (65) as

\left({}^{t}\vec{Q}_{\rm A_{2}}\bm{\tilde{M}}\vec{\psi}_{\rm A_{2}}\right)\vec{\psi}_{\rm A_{2}}\approx\lambda_{0}\vec{\psi}_{\rm A_{2}}.

(66)

This is equivalent to assuming that all but the slowest modes of $\bm{\tilde{M}}$ have decayed. Secondly, it is assumed that $\vec{\psi}_{\rm A_{2},\infty}\approx\vec{\psi}_{\rm A_{2},nr}$ . $\vec{\psi}_{\rm A_{2},\infty}$ represents the steady solution to the pre-QSS master equations (Eq. (60)), whereas $\vec{\psi}_{\rm A_{2},nr}$ is the solution to the QSS master equations (Eq. (26)) in the absence of recombination. Due to the fact that $\bm{M}$ and $\bm{\tilde{M}}$ are not exactly equal, $\vec{\psi}_{\rm A_{2},\infty}$ and $\vec{\psi}_{\rm A_{2},nr}$ are not formally equivalent. For practical applications however, they are the same. Thirdly, the third-body M is treated as inert such that $n_{\rm M}$ is constant.

By multiplying Eq. (64) by $n_{\rm M}$ and adding it to Eq. (63), then applying Eq. (66), the following expression is obtained

\frac{d\vec{\psi}_{\rm A_{2}}}{dt}=-n_{\rm M}\left(\bm{\tilde{M}}-\lambda_{0}\bm{I}\right)\left(\vec{\psi}_{\rm A_{2}}-\vec{\psi}_{\rm A_{2},\infty}\right).

(67)

Using the fact that $d\vec{\psi}_{\rm A_{2},nr}/dt=0$ and applying the approximation that $\vec{\psi}_{\rm A_{2},\infty}\approx\vec{\psi}_{\rm A_{2},nr}$ then gives

\frac{d}{dt}\left(\vec{\psi}_{\rm A_{2}}-\vec{\psi}_{\rm A_{2},nr}\right)=-n_{\rm M}\left(\bm{\tilde{M}}-\lambda_{0}\bm{I}\right)\left(\vec{\psi}_{\rm A_{2}}-\vec{\psi}_{\rm A_{2},nr}\right),

(68)

where the time-dependent solution is given by

\vec{\psi}_{\rm A_{2}}=\vec{\psi}_{\rm A_{2},nr}-\exp\left[-n_{\rm M}(\bm{\tilde{M}}-\lambda_{0}\bm{I})t\right]\left(\vec{\psi}_{\rm A_{2},nr}-\vec{\psi}_{\rm A_{2},0}\right).

(69)

Under the discussed assumptions, the pre-QSS internal state distribution of $\rm A_{2}$ is only a function of $n_{\rm M}$ , $T_{\rm t}$ , and $t$ . In particular, the pre-QSS solution has a simple exponential dependence on both $n_{\rm M}$ and $t$ . Mathematically, the convergence rate of this solution is controlled by the smallest eigenvalue of $(\bm{\tilde{M}}-\lambda_{0}\bm{I})$ Haug, Truhlar, and Blais (1987). Hence, this matrix can be approximated by $(\lambda_{1}-\lambda_{0}){\bm{I}}$ , where $\lambda_{1}$ is the second smallest eigenvalue of $\bm{\tilde{M}}$ . As will be shown in section III.1, the resulting functional form from this approximation is sufficient for capturing the majority of the pre-QSS behavior for $\rm H_{2}$ dissociation. Then, Eq. (69) is substituted into the general expression for $k_{\rm d}$ given by Eq. (34) to yield

	$\displaystyle k_{\rm d,pre-QSS}$		$\displaystyle=k_{\rm d,nr}-\exp(-n_{\rm M}(\lambda_{1}-\lambda_{0})t)(k_{\rm d,nr}-k_{\rm d,0})$		(70)
			$\displaystyle=k_{\rm d,nr}[1-(1-\epsilon)\exp(-n_{\rm M}(\lambda_{1}-\lambda_{0})t)].$		(70)

Here, $k_{\rm d,0}$ corresponds to the dissociation rate constant evaluated at the post-shock $T_{\rm t}$ but weighted by the pre-shock Boltzmann distribution, and $\epsilon\equiv k_{\rm d,0}/k_{\rm d,nr}$ . For any practical pre-shock temperatures (e.g., 100 - 300 K), $k_{\rm d,0}$ will be magnitudes smaller than $k_{\rm d,nr}$ , such that $\epsilon\ll 1$ .

While Eq. (70) gives a compact expression for $k_{\rm d,pre-QSS}$ , its dependence on $t$ makes it inconvenient for practical use. To remedy this, an alternate formulation is proposed where $t$ is replaced with the fraction of dissociation, $\phi_{\rm A}$ . This is done by considering the chemical source term for A in the pre-QSS,

\frac{dn_{\rm A}}{dt}=2n_{\rm A_{2}}n_{\rm M}k_{\rm d,pre-QSS}.

(71)

Equation (71) is then rewritten using Eq. (46) and Eq. (70) as

\frac{d\alpha}{dt}=(1-\alpha)n_{\rm M}k_{\rm d,nr}[1-(1-\epsilon)\exp(-n_{\rm M}(\lambda_{1}-\lambda_{0})t)],

(72)

where $\alpha$ is the degree of dissociation as defined previously in Eq. (45). Next, using a first-order Taylor expansion to approximate $\exp(-n_{\rm M}(\lambda_{1}-\lambda_{0})t)\approx 1-n_{\rm M}(\lambda_{1}-\lambda_{0})t$ , the solution to Eq. (72) (in the limit of a negligible $\epsilon$ ) is given by

t=\frac{1}{n_{\rm M}}\sqrt{-\frac{2\ln(1-\alpha)}{k_{\rm d,nr}(\lambda_{1}-\lambda_{0})}}.

(73)

Substituting this result back into Eq. (70) gives

k_{\rm d,pre-QSS}=k_{\rm d,nr}\left[1-(1-\epsilon)\exp\left(-\sqrt{-\frac{2(\lambda_{1}-\lambda_{0})\ln(1-\alpha)}{k_{\rm d,nr}}}\right)\right].

(74)

Finally, because the equilibrium degree of dissociation, $\alpha_{\rm eq}$ , is close to 1 (i.e., complete dissociation) for the majority of the conditions where pre-QSS effects would be relevant, $\alpha/\alpha_{\rm eq}=\phi_{\rm A}$ is approximated as $\alpha\approx\phi_{\rm A}$ . Performing this substitution in Eq. (74) yields the final expression for the pre-QSS dissociation rate constant,

k_{\rm d,pre-QSS}=k_{\rm d,nr}\left[1-(1-\epsilon)\exp\left(-\sqrt{-\frac{\ln(1-\phi_{\rm A})}{\eta}}\right)\right],

(75)

where $\eta(T_{\rm t})\equiv k_{\rm d,nr}/(2(\lambda_{1}-\lambda_{0}))$ . For most uses, the $\epsilon$ in Eq. (75) can be neglected. However, when Eq. (75) is used in a numerical integration (as in a CFD calculation), a non-zero value of $\epsilon$ is required such that $k_{\rm d,pre-QSS}$ is not 0 at $\phi_{\rm A}=0$ . Practically, for the case of $\rm H_{2}$ dissociation discussed later in section III.3, it was found that there is no sensitivity to the value of $\epsilon$ below the threshold of $\epsilon\leq 10^{-3}$ .

A key implication of these results is that under the discussed assumptions, $k_{\rm d,pre-QSS}$ can be expressed simply as $k_{\rm d,nr}$ modified by a multiplicative non-equilibrium factor that is only a function of $T_{\rm t}$ and $\phi_{\rm A}$ . Importantly, the dependence of $k_{\rm d,pre-QSS}$ on $n_{\rm M}$ as seen in Eq. (70) drops out when $t$ is substituted by $\phi_{\rm A}$ . This means that as long as an accurate fit of $\eta(T_{\rm t})$ can be obtained, the fraction of dissociation that occurs in the pre-QSS can be estimated easily as just a function of $T_{\rm t}$ . Analogously to the QSS analysis presented in section II.3.4, the fraction of dissociation that occurs in the pre-QSS region can be defined as the values of $\phi_{\rm A}=\alpha/\alpha_{\rm eq}$ for which $k_{\rm d,pre-QSS}$ is within some $\delta$ fraction of $k_{\rm d,nr}$ , i.e., $k_{\rm d,pre-QSS}=k_{\rm d,nr}(1-\delta)$ . Substituting this constraint into Eq. (75) with $\epsilon=0$ gives the desired result,

\phi_{\rm A}\big{|}_{\rm pre-QSS}=\frac{\alpha}{\alpha_{\rm eq}}\biggr{|}_{\rm pre-QSS}=1-\exp(-\eta(\ln(\delta))^{2}).

(76)

II.4.3 Chemical Source Term

There are two key implications of the analysis presented above for the calculation of the chemical source term. First, since Eq. (75) gives that $k_{\rm d,pre-QSS}$ is simply $k_{\rm d,nr}$ multiplied by a correction factor, it is easy to incorporate the expression for $k_{\rm d,pre-QSS}$ into the previous QSS rate constant expression. Namely, by replacing $k_{\rm d,nr}$ in Eq. (38) with $k_{\rm d,pre-QSS}$ , a complete description of $k_{\rm d}$ that is valid from the pre-QSS to the thermal limit can be obtained. This result is given by

k_{\rm d}=k_{\rm d,nr}\left[1-(1-\epsilon)\exp\left(-\sqrt{-\frac{\ln(1-\phi_{\rm A})}{\eta}}\right)\right]\left(1-\frac{\phi_{\rm A}^{2}}{\chi}\right)+k_{\rm d,th}\frac{\phi_{\rm A}^{2}}{\chi}.

(77)

As expected, in the limit of no dissociation when $\phi_{\rm A}=\chi=0$ , this expression gives that $k_{\rm d}=\epsilon k_{\rm d,nr}=k_{\rm d,0}$ , and in the limit of equilibrium when $\phi_{\rm A}=\chi=1$ , this expression gives that $k_{\rm d}=k_{\rm d,th}$ .

Second, by combining this description of $k_{\rm d}$ with the formula for $k_{\rm r}$ from Eq. (18), the general source terms for $\rm A_{2}$ and A can be computed directly from Eq. (2) and (3). Equivalently, by taking advantage of the same source term equivalence as discussed for Eq. (40) and (41), the chemical source term can be computed practically as

\frac{dn_{\rm A}}{dt}=2n_{\rm A_{2}}n_{\rm M}k_{\rm d,pre-QSS}-2n_{\rm A}^{2}n_{\rm M}\frac{k_{\rm d,pre-QSS}}{K_{\rm eq}}.

(78)

This gives a complete description of the non-equilibrium chemical source term that is only a function of the translational temperature, $T_{\rm t}$ , and the fraction of dissociation, $\phi_{\rm A}$ , via the expression for $k_{\rm d,pre-QSS}$ from Eq. (75).

II.4.4 Approximate Solution for $\bm{e}_{\bm{r}v}$

Analogously to the QSS case, the evolution of the average rovibrational energy, $e_{\rm rv}$ , can also be computed from the approximate solution of the $\rm A_{2}\it(J,\nu)$ distribution in the pre-QSS regime. Using the same substitutions as Eq. (75) for $\bm{\tilde{M}}$ and $t$ , the rovibrational distribution of $\rm A_{2}$ can be computed from Eq. (69) as

\vec{\psi}_{\rm A_{2}}=\vec{\psi}_{\rm A_{2},nr}-\exp\left(-\sqrt{-\frac{\ln(1-\phi_{\rm A})}{\eta}}\right)\left(\vec{\psi}_{\rm A_{2},nr}-\vec{\psi}_{\rm A_{2},0}\right).

(79)

Substituting this result into Eq. (51) and using the definition of $e_{\rm rv,nr}$ from Eq. (53) then gives

e_{\rm rv,pre-QSS}=e_{\rm rv,nr}\left[1-\left(1-\frac{e_{\rm rv,th,0}}{e_{\rm rv,nr}}\right)\exp\left(-\sqrt{-\frac{\ln(1-\phi_{\rm A})}{\eta}}\right)\right],

(80)

where $e_{\rm rv,th,0}$ corresponds to $e_{\rm rv,th}$ evaluated via Eq. (54) at the pre-shock temperature. Similarly to the discussion for $k_{\rm d,pre-QSS}$ , this result gives that $e_{\rm rv,pre-QSS}$ is equal to $e_{\rm rv,nr}$ multiplied by a pre-QSS correction factor. Therefore, by applying this correction factor in the previous QSS description of $e_{\rm rv}$ from Eq. (52), the corresponding complete description of $e_{\rm rv}$ from the pre-QSS to the thermal limit is given by

e_{\rm rv}=e_{\rm rv,nr}\left[1-\left(1-\frac{e_{\rm rv,th,0}}{e_{\rm rv,nr}}\right)\exp\left(-\sqrt{-\frac{\ln(1-\phi_{\rm A})}{\eta}}\right)\right]\left(1-\frac{\phi_{\rm A}^{2}}{\chi}\right)+e_{\rm rv,th}\frac{\phi_{\rm A}^{2}}{\chi}.

(81)

III Application to the Dissociation of $\rm\bf H_{2}$

In this section, the QSS and pre-QSS theories described in section II are applied to the case of $\rm H_{2}$ dissociation. In section III.1, the rate expressions are first validated against master equation results from the literature. Then, in section III.2, the fraction of dissociation that occurs in the QSS and pre-QSS regimes are estimated as a function of temperature and number density. Finally in section III.3, the QSS and pre-QSS source term expressions are used to compute number density and average rovibrational energy profiles which are compared directly to the master equation results.

III.1 Extracted Rate Constants

To validate the analysis presented in the previous sections, the predictions of the QSS and pre-QSS theories are compared to the results of the master equation calculations for the dissociation of $\rm H_{2}$ by Kim and Boyd Kim and Boyd (2012) (M = $\rm H_{2}$ ) and Kim Kim (2015) (M = H and He). In these computational studies, state-specific rate constants were evaluated using the QCT method. For the case of M = $\rm H_{2}$ by Kim and Boyd Kim and Boyd (2012), the complete set of state-specific rate constants for both target and projectile $\rm H_{2}$ molecules were evaluated using the QCT method along with response surfaces designed by the ordinary Kriging model. For each third-body, the full set of master equations was used to simulate a 0-D isothermal and isochoric reactor in a heating environment where $T_{\rm t}$ and $\tilde{n}_{\rm H}$ were held constant. $\rm H_{2}$ molecules were allowed to dissociate, and the evolution to equilibrium of the initial rovibrational distribution of $\rm H_{2}$ (characterized by $T_{\rm r,0}=T_{\rm v,0}<T_{\rm t}$ ) was observed. The resulting number density and internal temperature profiles from these calculations are reproduced in Fig. 1 and 2, respectively¹¹1Number densities, rovibrational temperatures, and rate constants for these studies were only reported as plots. Therefore, these values have all been obtained via digitizations.. For the 10,000 K and 16,000 K cases with M = H and He, rovibrational temperatures are obtained from the earlier study by Kim et al. Kim, Kwon, and Park (2009), as temperatures were reported over a wider range of simulation times. For all plotted cases, rovibrational temperatures were initialized at $T_{\rm r,0}=T_{\rm v,0}$ = 1,000 K, and the number density of $\rm H_{2}$ was initialized at $1.0\times 10^{18}$ $\rm cm^{-3}$ . For the cases with M = H and He, the number density of the third-body was additionally set to $5.0\times 10^{17}$ $\rm cm^{-3}$ . For the higher temperature cases with all third-bodies, the results of Fig. 2 highlight distinct plateaus in the $T_{\rm r}$ and $T_{\rm v}$ profiles between the initial value at 1,000 K and the final equilibrium values at $T_{\rm t}$ . These plateaus correspond to the non-recombining QSS limit, and the regions before and after the plateaus correspond to the pre-QSS and recombining regions, respectively.

For several of these master equation cases, the corresponding aggregate $k_{\rm d}$ is calculated and shown as solid black lines in Fig. 3 ²²2For the case of M= $\rm H_{2}$ , a modified equilibrium constant, $K_{\rm eq,mod}=2K_{\rm eq}$ , was required to reproduce the equilibrium number densities plotted by Kim and Boyd Kim and Boyd (2012). Additionally, for the case of M=H, $k_{\rm d}$ values were extracted from the number density profiles assuming $n_{\rm M}$ was fixed and equal to $n_{\rm H,0}$ .. These values are computed by substituting $k_{\rm d,th}$ into Eq. (18) to compute $k_{\rm r}$ , after which $k_{\rm r}$ and the corresponding number densities from Fig. 1 and their derivatives are substituted into Eq. (3) to compute $k_{\rm d}$ . For the case of M = $\rm H_{2}$ , $k_{\rm d,th}$ values are taken directly from Kim and Boyd Kim and Boyd (2012). However, in the cases of M = H and He, $k_{\rm d,th}$ values were not reported by Kim Kim (2015). Therefore, $k_{\rm d,th}$ values are instead taken from the recent QCT study by Vargas et al. Vargas, Monge-Palacios, and Lacoste (2024). The master equation-extracted $k_{\rm d}$ values show similar trends for all third-bodies. Namely, the rate constants start at a low value, increase to reach a plateau, and then equilibrate to the thermal limit. As with the $T_{\rm r}$ and $T_{\rm v}$ plots from Fig. 2, the initial increase corresponds to the pre-QSS period, the plateau corresponds to the non-recombining QSS limit, and the final increase corresponds to the recombining region. For the lower temperature cases, the gradual evolution of $k_{\rm d}$ from the non-recombining QSS to the thermal limit is apparent; at higher temperatures this transition occurs rapidly near $\phi_{\rm H}=1$ . Consistent with the QSS theory discussed previously, these results show that $k_{\rm d}$ only reaches $k_{\rm d,th}$ when $\phi_{\rm H}=1$ . In other words, thermal and chemical equilibria are reached simultaneously. At the other end near $\phi_{\rm H}=0$ , the pre-QSS region is negligible for the lower temperature cases. As $T_{\rm t}$ is increased, the fraction of dissociation that occurs in the pre-QSS region increases.

To compare to the master equation results, predicted $k_{\rm d}$ values from the QSS as well as the pre-QSS theories are computed. These are shown in Fig. 3 as dash-dotted red and dashed blue lines, respectively. For the QSS predicted values, $k_{\rm d}$ is computed using Eq. (38); for the QSS values with the pre-QSS correction, Eq. (77) is used instead with $\epsilon=0$ . The $\eta$ values used in the pre-QSS expressions are computed for each third-body and temperature case using a least squares fit to the master equation results. As expected, the $k_{\rm d}$ values from the QSS theory without the pre-QSS correction are able to reproduce the results of the master equations everywhere outside of the pre-QSS region. Namely, the $k_{\rm d,nr}$ values reflect the plateaus in the master equation results, and the transition from the non-recombining QSS to the thermal limit is captured correctly for all third-bodies and temperatures by the QSS theory. Additionally, with the simple pre-QSS correction of Eq. (75), the $k_{\rm d}$ values are also able to reproduce the results of the master equation calculations in the pre-QSS regions.

III.2 Fraction of Dissociation in the Pre-QSS and QSS Regions

To estimate the fraction of dissociation that occurs in the pre-QSS for any temperature, fits of $\eta(T_{\rm t})$ for each third-body are first computed. The definition of $\eta(T_{\rm t})\equiv k_{\rm d,nr}/(2(\lambda_{1}-\lambda_{0}))$ gives that $\eta$ is effectively the ratio of two rate constants, as $\lambda_{0}$ and $\lambda_{1}$ both correspond to eigenvalues of the rate matrix $\bm{\tilde{M}}$ . Therefore, a modified Arrhenius form is used for these fits. The parameters for the Arrhenius fits are determined by fitting $\eta(T_{\rm t})$ as a function of temperature (in K) across each of the master equation cases discussed in the previous section. For M = $\rm H_{2}$ , H, and He, these fits are given by $\eta_{\rm M=H_{2}}(T_{\rm t})=6.489\times 10^{1}T_{\rm t}^{-0.50}\exp(-2.221\times 10^{4}/T_{\rm t})$ , $\eta_{\rm M=H}(T_{\rm t})=1.897\times 10^{3}T_{\rm t}^{-0.66}\exp(-4.236\times 10^{4}/T_{\rm t})$ , and $\eta_{\rm M=He}(T_{\rm t})=1.668\times 10^{6}T_{\rm t}^{-1.41}\exp(-4.035\times 10^{4}/T_{\rm t})$ , respectively. Figure 4 compares these fits (in solid lines) to the master equation-extracted $\eta$ values (in open symbols). The Arrhenius fits are in excellent agreement with the master equation-extracted $\eta$ values. In fact, using these Arrhenius fits in Eq. (77) instead of the individual master equation-extracted $\eta$ values produced effectively identical predictions of the pre-QSS dissociation rate constants for the master equation cases discussed in the previous section.

Using these $\eta(T_{\rm t})$ fits, the fraction of dissociation that occurs in the pre-QSS, non-recombining QSS, and recombining regions as a function of $T_{\rm t}$ and $\tilde{n}_{\rm H}$ is estimated using Eq. (44), (50), and (76). In these equations, the exact value of $\delta$ and hence the cutoffs between the pre-QSS, non-recombining QSS, and recombining regions is arbitrary. As will be shown in section IV, the uncertainties in the rate constants from the literature are approximately bounded by a factor of two. In light of these uncertainties, a $\delta$ value of 0.25 is used for the present analysis. For the M = $\rm H_{2}$ case, continuous fits of $k_{\rm d,nr}$ and $k_{\rm d,th}$ are obtained by fitting the reported values from Kim and Boyd Kim and Boyd (2012) to a modified Arrhenius form. For the M = H and He cases, fits of $k_{\rm d,nr}$ and $k_{\rm d,th}$ are taken from Kim et al. Kim, Kwon, and Park (2009) and Vargas et al. Vargas, Monge-Palacios, and Lacoste (2024), respectively. The results of these calculations are shown in Fig. 5. The selected temperature and number density ranges correspond to the range of conditions relevant for post-shock hypersonic entry flows as well as the shock tube experiments discussed later in section IV.1. The $\tilde{n}_{\rm H}=2.0\times 10^{18}$ $\rm cm^{-3}$ case in particular corresponds to the master equation calculations discussed in the previous section. These plots highlight that for most temperatures and number densities, the majority of dissociation for all three third-bodies occurs near the non-recombining QSS limit where $k_{\rm d}\approx k_{\rm d,nr}$ . As the temperature decreases and/ or as the number density increases, the fraction of dissociation that occurs in the recombining region increases. Luckily, as discussed in section II.3.3, the source terms for the entire QSS region can be captured by a one-temperature fit of $k_{\rm d,nr}(T_{\rm t})$ , as the dependence on number densities for the transition of $k_{\rm d}$ between the non-recombining and thermal limits is captured implicitly by Eq. (41). Following Eq. (76), the fraction of dissociation that occurs in the pre-QSS region is only a function of temperature. This is consistent with the empirical observations from the $\rm O_{2}$ + O master equation calculations of Venturi et al. Venturi et al. (2020). While the fraction of dissociation that occurs in the pre-QSS is negligible at low temperatures, it increases monotonically with temperature, and in the case of M = H, becomes comparable to the non-recombining contribution at 20,000 K.

III.3 Predicted Number Densities and Rovibrational Energies

The QSS source term expressions with and without the pre-QSS correction (Eq. (78) and (41), respectively) are numerically integrated, and the resulting number density profiles are compared to those of the master equation calculations in Fig. 6. For the pre-QSS expression, $\epsilon$ is set to $10^{-3}$ and the $\eta(T_{\rm t})$ fits obtained in the previous section are used. For the lowest temperature cases (4,000 K and 6,000 K), both sets of the predicted profiles reproduce the master equation results accurately for all three third-bodies. This serves as a verification of the source term equivalence discussed in section II.3.3, as Fig. 3 shows that the aggregate $k_{\rm d}$ varies considerably between the $k_{\rm d,nr}$ and $k_{\rm d,th}$ limits for these low temperature cases, and yet integrating Eq. (41) and (78) (which are only functions of $k_{\rm d,nr}$ and $k_{\rm d,pre-QSS}$ , respectively) still yields the correct number density profiles. For the higher temperature cases ( $T_{\rm t}\geq$ 10,000 K), a larger difference appears between the predicted profiles with and without the pre-QSS correction. In particular, the QSS expression without the pre-QSS correction results in dissociation that is too fast when compared to the master equation calculations. However, with the pre-QSS correction, the predicted profiles are in much better agreement with the master equation results, especially for the M = H and He cases where there is excellent agreement for the higher temperature cases.

Figure 7 shows the corresponding average rovibrational energy profiles for the same simulation cases. For the master equation calculations, the reported $T_{\rm r}$ and $T_{\rm v}$ values from Fig. 2 are used with Eq. (55)-(57) to compute $e_{\rm rv}$ values. For the predicted values, $T_{\rm r,nr}$ and $T_{\rm v,nr}$ are first estimated from Fig. 2, then are used in Eq. (58) to compute $e_{\rm rv,nr}$ . Then, values of $e_{\rm rv}$ are computed using either Eq. (81) or Eq. (52) with the predicted number densities from Fig. 6 for the cases with and without the pre-QSS correction, respectively. Overall, the predicted profiles with the pre-QSS correction reproduce the trends seen in the master equation calculations, without any refitting of $\eta(T_{\rm t})$ . In the recombining regions for all third-bodies, there is excellent agreement in the transition of $e_{\rm rv}$ between the $e_{\rm rv,nr}$ and $e_{\rm rv,th}$ limits.

The results of Fig. 6 and 7 show that the relatively simple pre-QSS corrected expressions of Eq. (75) and (81) are sufficient for capturing the majority of the non-equilibrium behavior of the full master equation calculations for $\rm H_{2}$ dissociation. Importantly, the rate constant expression only requires fits for $k_{\rm d,nr}(T_{\rm t})$ and $\eta(T_{\rm t})$ , and is ultimately only a function of $T_{\rm t}$ and the fraction of dissociation, $\phi_{\rm H}=n_{\rm H}/n_{\rm H,eq}$ . In CFD calculations, $T_{\rm t}$ and $n_{\rm H}$ are already computed via the energy and species’ continuity equations. Therefore, the non-equilibrium dissociation behavior (including QSS and pre-QSS effects) may be captured with just a one-temperature formulation that only considers bulk species. This approach offers benefits over other approaches from the literature.

The most widely used class of non-equilibrium dissociation models in the literature is known broadly as “multi-temperature” models. Examples of such models include the two-temperature formulation by Park Park (1988, 1989a, 1989b), as well as the more recent QCT-informed modified Marrone-Treanor (MMT) Chaudhry and Candler (2019); Chaudhry et al. (2020) and Singh-Schwartzentruber models Singh and Schwartzentruber (2020a, b). Multi-temperature models require the solution of additional energy equations (one per independently tracked energy mode/ temperature) to describe the coupling between internal energy relaxation and dissociation. As a consequence, multi-temperature models (at a minimum) require inputs for relaxation time-scales and energy feedback factors, as well as multi-temperature fits of rate constants. These three model inputs are inherently coupled and cannot be estimated independently. In particular, they must be evaluated in a consistent manner to recover QSS dissociation rates accurately. By contrast, this coupling is already captured by design in the one-temperature formulation of the present work, as the expressions for $k_{\rm d}$ (Eq. (75)) and $e_{\rm rv}$ (Eq. (81)) are derived from the same non-equilibrium internal energy distribution solution.

A more general class of non-equilibrium dissociation models from the literature are based on the multi-group/ coarse-graining approach outlined by Liu et al. Liu et al. (2015). In this approach, individual rovibrational states are binned into a number of “groups”. Typically, this is done by assuming that the states within a given group follow a Boltzmann distribution about a locally defined internal group temperature. Then, the master equations are solved for the grouped states instead of the entire set of rovibrational states. In the limit where each “group” consists of a single rovibrational state, the full master equation/ state-to-state approach is recovered. Conversely, in the limit where all rovibrational states are binned together into one group, the multi-temperature approach is recovered. While this approach has been used successfully in the past to reproduce the dissociation dynamics of full master equation calculations with a fraction of the species Magin et al. (2012); Munafò, Panesi, and Magin (2014); Sahai et al. (2017, 2019); Venturi et al. (2020), there are a few limitations. Firstly, similar to the multi-temperature approach, the coarse-graining approach requires the solution of additional continuity and internal energy conservation equations for each group of states. Even if relatively few groups are required, the associated computational cost is still larger than the one-temperature approach. Secondly, finding the optimal grouping strategy to capture the correct physics while minimizing computational cost is not a trivial task Sahai et al. (2017). This is not a problem for the present formulation, as only bulk species need to be tracked explicitly.

IV Review of Rate Constants for $\rm\bf H_{2}$ Dissociation

While the master equation studies of Kim and Boyd Kim and Boyd (2012) and Kim Kim (2015)/ Kim et al. Kim, Kwon, and Park (2009) discussed in the previous section provide a good validation of the QSS and pre-QSS theories and framework, they only represent one set of calculations for determining $\rm H_{2}$ dissociation rate constants. Even detailed master equation studies such as those by Kim and Boyd Kim and Boyd (2012) and Kim Kim (2015)/ Kim et al. Kim, Kwon, and Park (2009) can have associated uncertainties due to potential inaccuracies in the underlying PESs or QCT calculations. As a consequence, accurate estimates of rate constants should rely on compilations of data from as many sources as possible. To implement the rate constant models proposed in section II.4, only fits of $k_{\rm d,nr}(T_{\rm t})$ and $\eta(T_{\rm t})$ for each third-body are required. Since $\eta(T_{\rm t})$ can only be computed by fitting directly to master equation results, the expressions presented previously in section III.2 already represent the best available estimates. However, estimates of $k_{\rm d,nr}(T_{\rm t})$ can be obtained from a much wider variety of sources. Therefore, in this section, rate constants for the dissociation of $\rm H_{2}$ from both experimental and computational studies are reviewed for each relevant third-body. Then, from this review, new fits of $k_{\rm d,nr}(T_{\rm t})$ that are valid from 200 - 20,000 K are proposed.

The available data sets in the literature can be split into three primary categories: high temperature (2,000 - 8,000 K) rate constants extracted from shock tube experiments (reviewed in section IV.1), low temperature ( $\lesssim$ 350 K) rate constants extracted from discharge-flow tube experiments (reviewed in section IV.2), and computational studies (reviewed in section IV.3). A comparison of the rate constants from these sources along with new fits of the rate constants for each third-body are presented in section IV.4.

IV.1 High Temperature/ Shock Tube Experiments

Table 1 summarizes the relevant shock tube studies in the literature for $\rm H_{2}$ dissociation. While there were differences in the measurement techniques and measured quantities used to determine the degree of dissociation in these studies, the overall method used to compute rate constants was the same. Namely, measurements were taken behind incident shocks of varying strengths and compositions, for mixtures composed of $\rm H_{2}$ and a heavy noble gas (i.e., Ar, Xe, or Kr). Then, corresponding inviscid 1-D shock calculations were performed with varying rate constants assuming the linear mixture rule (see Appendix A) until the calculations could reproduce the measured results. For most studies, this was done by assuming that the dissociation rate constants followed the modified Arrhenius form given by

k_{\rm d}=A_{\rm d}T^{n_{\rm d}}\exp\left(-\frac{\theta_{\rm d,H_{2}}}{T}\right),

(82)

or alternatively

k_{\rm r}=A_{\rm r}T^{n_{\rm r}}

(83)

for recombination rate constants. Here, the scalars $A_{\rm d,r}$ and $n_{\rm d,r}$ are fitting parameters specific to each third-body. Because these experiments were conducted behind strong incident shocks, the $\rm H_{2}$ dissociation process was dominant over recombination. Nevertheless, rate constants were often still reported in the recombination direction instead of dissociation by assuming $k_{\rm r}=k_{\rm d}/K_{\rm eq}$ .

Table 1: Experimental Sources for High Temperature Rate Constants

Source	Measurement Method	M	$T_{\rm t}$ [K]	$P$ [atm]
Gardiner and Kistiakowsky, 1961 Gardiner and Kistiakowsky (1961)	X-ray densitometry	$\rm H_{2}$ , H, Xe	3,000 - 4,500	3.6e-1 - 6.2e-1
Sutton, 1962 Sutton (1962)	light interferometry	$\rm H_{2}$ ¹¹1Rate constants were obtained from the fits by Baulch et al. Baulch et al. (2005)., H²²2Rate constants were obtained via digitizations of the published plots., Ar¹¹1Rate constants were obtained from the fits by Baulch et al. Baulch et al. (2005).	2,800 - 4,500	2.6e-2 - 1.3e-1
Rink, 1962 Rink (1962)	X-ray densitometry	$\rm H_{2}$ , H, Ar, Xe, Kr	2,800 - 5,000	1.2e-2 - 6.7e-2
Patch, 1962 Patch (1962)	UV absorption ( $\rm H_{2}$ )	$\rm H_{2}$ , H, Ar	2,950 - 5,330	2.5e-3 - 1.1e-1
Jacobs, Geidt, and Cohen, 1967 Jacobs, Giedt, and Cohen (1967)	IR emission (trace HCl)	$\rm H_{2}$ , H, Ar	2,900 - 4,700	1.4e0 - 1.8e0
Myerson and Watt, 1968 Myerson and Watt (1968)	VUV absorption (H)	$\rm H_{2}$ , Ar	2,290 - 3,790	2.0e-1 - 1.0e0
Hurle, Jones, and Rosenfeld, 1969 Hurle, Jones, and Rosenfeld (1969)	Na or Cr line-reversal	$\rm H_{2}$ , H²²2Rate constants were obtained via digitizations of the published plots., Ar	2,500 - 7,000	6.6e-3 - 2.6e-2
Breshears and Bird, 1973 Breshears and Bird (1972)	laser schlieren	$\rm H_{2}$ , H, Ar, Xe	3,500 - 8,000	3.9e-3 - 6.6e-2

For all of the shock tube studies, one major source of uncertainty arises from the fact that the rate constants were fitted simultaneously for all three third-bodies ( $\rm H_{2}$ , H, and the chosen noble gas). Therefore, an underprediction of a rate constant for one third-body would result in an overprediction of the rate constant for a different third-body. Additionally, due to the limited temperature ranges, Gardiner and Kistiakowsky Gardiner and Kistiakowsky (1961), Rink Rink (1962), Patch Patch (1962), and Jacobs et al. Jacobs, Giedt, and Cohen (1967) reported that the temperature exponential term in the Arrhenius fits could not be determined reliably. As a consequence, these four studies all assumed a value of $n_{\rm r}=-1$ for their reported $k_{\rm r}$ fits with all third-bodies. Therefore, while the magnitude of the reported rate constants from these four shock tube studies may be meaningful, the temperature dependencies may not be as accurate across larger temperature ranges. This will be discussed further in the context of other low temperature and computational data sets in section IV.4.

Based on the QSS and pre-QSS theories presented earlier in section II.3 and II.4, a relevant question is what impact non-equilibrium effects may have had on the reported rate constants from these shock tube studies. Based on Fig. 5, it is clear that pre-QSS effects were likely negligible even at the highest temperatures (7,000 - 8,000 K) investigated in the shock tube studies. However, the estimated fraction of dissociation that occurred in the non-recombining versus recombining regions varies significantly based on the number density conditions. To illustrate this point, Fig. 8 shows contours of the estimated fraction of dissociation that occurred in the non-recombining limit, computed using Eq. (44) and (50) for the case of M = $\rm H_{2}$ . These limits are computed with $k_{\rm d,nr}$ and $k_{\rm d,th}$ taken from Kim and Boyd Kim and Boyd (2012), and with $\delta$ = 0.25. Overlaid are the frozen $T_{\rm t}$ and $\tilde{n}_{\rm H}$ values from the shot conditions in the studies by Gardiner and Kistiakowsky Gardiner and Kistiakowsky (1961), Rink Rink (1962), and Patch Patch (1962), computed using the CEA code Gordon and McBride (1996). For the sake of clarity, only the shots from these three studies are plotted, however, the shot conditions from any of the other studies in Table 1 could be plotted similarly. From Fig. 8, it is clear that the fraction of dissociation that occurred in the non-recombining limit varied drastically for different shot conditions. For example, for some of the highest temperature cases by Patch, 95% of the dissociation is estimated to have occurred in the non-recombining limit; however, for the lowest temperature cases by Gardiner and Kistiakowski, this is just 15%.

Luckily, in all of the considered shock tube studies, it was assumed that dissociation and recombination rate constants were related via macroscopic detailed balance. As discussed in section II.3.3, and in particular as given by Eq. (41), this implies that the reported rate constant values from the shock tube studies correspond to $k_{\rm d,nr}$ (and $k_{\rm d,nr}/K_{\rm eq}$ where recombination rate constants were reported instead), regardless of the temperature and number density conditions in which the experiments were conducted.

IV.2 Low Temperature/ Discharge-Flow Tube Experiments

Table 2 summarizes the discharge-flow tube studies considered in the present work for the determination of rate constants at low temperatures ( $T_{\rm t}\lesssim$ 350 K). These experiments relied on a conventional dynamical flow tube setup, wherein atomic hydrogen was first generated, then continually pumped along a tube with measurements taken as H atom concentrations decayed with distance. The only exception to this is the study by Lynch et al. Lynch, Schwab, and Michael (1976), which instead relied on a static method with H atom decay measured instead as a function of time. It should be noted that there are a number of earlier discharge-flow tube studies in the literature (prior to 1960) that show a wide range of possible recombination rates near $T_{\rm t}\approx$ 300 K with approximately a magnitude of disagreement for the third-bodies M = $\rm H_{2}$ and Ar. However, the previous reviews by Cohen and Westburg Cohen and Westberg (1983) and Baulch et al. Baulch et al. (2005) suggest that the rate constants reported by these earlier investigations are in doubt, and therefore they are not considered in the present work.

Table 2: Experimental Sources for High Temperature Rate Constants

Source	Measurement Method	M	$T_{\rm t}$ [K]	$P$ [atm]
Larkin, 1968 Larkin (1968)	calorimetric probe	$\rm H_{2}$	291	8.4e-3
		Ar	213 - 349	6.0e-3 - 7.1e-3
Bennett and Blackmore, 1968 Bennett and Blackmore (1968)	ESR spectroscopy¹¹1Electron Spin Resonance spectroscopy.	$\rm H_{2}$ , H²²2Only an upper bound value was provided.	300	1.3e-3 - 1.3e-2
Bennett and Blackmore, 1970 Bennett and Blackmore (1970)	ESR spectroscopy¹¹1Electron Spin Resonance spectroscopy.	$\rm H_{2}$	300	1.8e-3 - 7.7e-3
Bennett and Blackmore, 1971 Bennett and Blackmore (1971)	ESR spectroscopy¹¹1Electron Spin Resonance spectroscopy.	$\rm H_{2}$ , He, Ar	298	6.5e-2 - 3.3e-1
Trainor, Ham, and Kaufman, 1973 Trainor, Ham, and Kaufman (1973)	calorimetric probe	$\rm H_{2}$ , He, Ar	77 - 298	2.6e-3 - 2.0e-2
Lynch, Schwab, and Michael, 1976 Lynch, Schwab, and Michael (1976)	H Lyman- $\alpha$ absorption	$\rm H_{2}$ , He, Ne, Ar, Kr	298	6.5e-1 - 2.0e0
Mitchell and LeRoy, 1977 Mitchell and LeRoy (1977)	capillary flow-meter	He³³3Both lower and upper bound values were reported.	296.5	7.9e-3 - 2.4e-2

Rate constants are available from the listed sources for $\rm H_{2}$ and the noble gases He, Ar, and Kr as third-bodies. Unfortunately, due to low temperature limitations, rates with H as a third-body could not reliably be determined using discharge-flow tube systems. Bennett and Blackmore Bennett and Blackmore (1968) proposed a room temperature recombination rate for M = H, but only as an upper bound. However, this value has been called into question by Baulch et al. Baulch (1972).

As was done implicitly in the previous reviews by Cohen and Westburg Cohen and Westberg (1983) and Baulch et al. Baulch et al. (2005), it is assumed that at the low temperatures encountered in the discharge-flow tube studies ( $T_{\rm t}\lesssim$ 350 K), there is no distribution dependence for the reported rate constants. As a consequence, the recombination rate constants reported in these studies are interpreted directly as estimates of $k_{\rm d,nr}/K_{\rm eq}$ . Additionally, it is assumed that recombination is a third-order reaction. The validity of this assumption is discussed in Appendix B.

IV.3 Computational Studies

Table 3 summarizes the computational studies considered in the present work. Each of these studies followed a similar methodology to that already discussed previously in section III.1 for the works of Kim and Boyd Kim and Boyd (2012) and Kim Kim (2015)/ Kim et al. Kim, Kwon, and Park (2009). Namely, $(J,\nu)$ state-specific relaxation and reaction rates were first computed using QCT calculations on detailed ab-initio PESs. Once converged state-specific rates were obtained, aggregate rate constants were computed. The specific method used to compute these aggregate rate constants differs slightly between the different studies.

Table 3: Computational Sources for Rate Constants

Source	M	$T_{\rm t}$ [K]	PES
Schwenke, 1990 Schwenke (1990)	$\rm H_{2}$	1,000 - 5,000	Schwenke, 1988 Schwenke (1988)
	H	3,000 - 5,000	Schwenke, 1988 Schwenke (1988)
Martin, Schwarz, and Mandy, 1996 Martin, Schwarz, and Mandy (1996)	H¹¹1The sum of dissociation via both collision-induced and tunneling contributions is considered in the present work.	450 - 45,000	Liu, Siegbahn, Truhlar, and Horowitz, 1978 Liu (1973); Siegbahn and Liu (1978); Truhlar and Horowitz (1978, 1979)
Furudate, Fujita, and Abe, 2006 Furudate, Fujita, and Abe (2006)	$\rm H_{2}$ ²²2Rate constants were obtained via digitizations of the published plots.³³3 $k_{\rm d,nr}$ was computed using two methods. Only the first (equivalent to Eq. (37)) is considered in the present work.	5,000 - 50,000	Schwenke, 1988 Schwenke (1988)
Kim, Kwon, and Park, 2009 Kim, Kwon, and Park (2009)	H	2,000 - 16,000	Boothroyd, Keogh, Martin, and Peterson, 1996 Boothroyd et al. (1996)
	He	2,000 - 16,000	Boothroyd, Martin, and Peterson, 2003 Boothroyd, Martin, and Peterson (2003)
Kim, Kwon, and Park, 2010 Kim, Kwon, and Park (2010)	$\rm H_{2}$	4,000 - 30,000	Schwenke, 1988 Schwenke (1988)
Kim and Boyd, 2012 Kim and Boyd (2012)	$\rm H_{2}$ ²²2Rate constants were obtained via digitizations of the published plots.	1,000 - 30,000	Schwenke, 1988 Schwenke (1988)

For the studies by Furudate et al. Furudate, Fujita, and Abe (2006), Kim and Boyd Kim and Boyd (2012), and both studies by Kim et al. Kim, Kwon, and Park (2009, 2010), aggregate dissociation rate constants were computed by using the same QSS formulation outlined in section II.3. As a consequence, the reported rate constants correspond directly to values of $k_{\rm d,nr}$ . For the study by Furudate et al. Furudate, Fujita, and Abe (2006), two sets of rate constants were reported: one based directly on QCT calculations, and one where the same underlying state-specific rate constants were modified by an arbitrary scaling factor to match the state-specific rate constants reported by Sharma Sharma (1994). These were labeled as the “uncorrected” and “corrected” rate constants, respectively. Both sets of rate constants are presented in section IV.4.

For the study by Martin et al. Martin, Schwarz, and Mandy (1996), recombination terms were not included in their formulation of the master equations, so their “steady-state” dissociation rate constants also correspond to $k_{\rm d,nr}$ . The authors reported two limiting dissociation rate constants: a “low density” rate constant for $n_{\rm H}\lesssim 1$ $\rm cm^{-3}$ and a “high density” rate constant for $n_{\rm H}\gtrsim 10^{4}$ $\rm cm^{-3}$ , with a smooth transition in between. This dependence on $n_{\rm H}$ appeared because radiative transitions were included in their calculations. For $n_{\rm H}\lesssim 10^{4}$ $\rm cm^{-3}$ , radiative transitions were non-negligible and served to depopulate the excited ${\rm H_{2}}(J,\nu)$ state distribution, thereby lowering $k_{\rm d,nr}$ . However, the values of $n_{\rm H}$ that are relevant for hypersonic entry calculations are generally several magnitudes greater than $n_{\rm H}\approx 10^{4}$ $\rm cm^{-3}$ . Therefore, only the “high density” fit of $k_{\rm d,nr}$ is considered in the present work.

Finally, for the study by Schwenke Schwenke (1990), a different approach was taken in which isothermal simulations were first performed for a $\rm H_{2}$ / H mixture using the master equations. Then, using several different approaches, corresponding phenomenological rate constants were extracted and compared. Of the tested methods, the most meaningful was determined to be the one in which $\rm H_{2}$ and H number densities from the master equation calculations were substituted into the aggregate source term equations, Eq. (2) and (3), assuming the linear mixture rule (see Appendix A) and assuming that $k_{\rm d}/k_{\rm r}=K_{\rm eq}$ . This method was described as the “small master equation” approach. As shown in section II.3, this approach is equivalent to fitting $k_{\rm d,nr}$ values using Eq. (41). Because the maximum temperature simulated by Schwenke Schwenke (1990) was 5,000 K, pre-QSS effects (as estimated by Fig. 5) were negligible. Hence, the rate constants reported by Schwenke can also be interpreted directly as values of $k_{\rm d,nr}$ .

Of the four studies that computed rate constants for M = $\rm H_{2}$ , only Kim and Boyd Kim and Boyd (2012) considered the complete set of state-to-state rate coefficients for both target and projectile $\rm H_{2}$ molecules. The other three studies assumed that the internal state distribution of the projectile $\rm H_{2}$ followed a Boltzmann distribution defined by either by the translational temperature (Schwenke Schwenke (1990) and Furudate et al. Furudate, Fujita, and Abe (2006)) or a separate non-equilibrium internal temperature (Kim et al. Kim, Kwon, and Park (2010)). As discussed by Kim and Boyd Kim and Boyd (2012), the reported dissociation rate constants from Kim et al. Kim, Kwon, and Park (2010) were still similar to those of Kim and Boyd Kim and Boyd (2012).

Outside of the master equation studies discussed above, there are a number of earlier studies in the literature that also relied on master equation calculations to compute $\rm H_{2}$ dissociation rate constants. However, these earlier studies generally made several more simplifying assumptions. These assumptions include the local equilibrium/ Boltzmann treatment of the rotational mode Shui, Appleton, and Keck (1971); Shui and Appleton (1971); Shui (1973); Roberge and Dalgarno (1982); Esposito, Gorse, and Capitelli (1999), the calculation of state-specific rate constants using simple scaling laws instead of trajectory calculations Dove and Jones (1972); Lepp and Shull (1983), or the explicit treatment of only para- $\rm H_{2}$ (as opposed to both ortho and para- $\rm H_{2}$ ) Blais and Truhlar (1979); Haug, Truhlar, and Blais (1987); Dove and Raynor (1979); Dove et al. (1987). Because the validity of these assumptions is questionable, these earlier studies are not considered for the review in the present work.

IV.4 Comparison and Fits of Rate Constant Data

Figure 9 shows the values of $k_{\rm d,nr}/K_{\rm eq}$ from all of the sources discussed in the previous sections. The rate constant data is plotted as $k_{\rm d,nr}/K_{\rm eq}$ instead of $k_{\rm d,nr}$ , as this significantly reduces the total change in magnitude of the plotted values between the low and high temperature regimes, making comparisons between different data sets easier. However, these plots should not be mistaken as plots of $k_{\rm r}$ , as $k_{\rm d,nr}/K_{\rm eq}\neq k_{\rm r}$ (see section II.3). Additionally, $k_{\rm d,nr}/K_{\rm eq}$ values are only plotted as functions of $T_{\rm t}$ , as both $k_{\rm d,nr}$ and $K_{\rm eq}$ are functions of $T_{\rm t}$ alone. For the cases with a noble gas as a third-body, i.e., M = He, Ne, Ar, Xe, and Kr, the rate constants have all been plotted together. While there have been some studies in the literature Shui and Appleton (1971); Whitlock, Muckerman, and Roberts (1974) that have indicated that the rate constants for $\rm H_{2}$ dissociation are not identical between different noble gas third-bodies (especially between He and the heavier noble gases), the scatter of the rate constant data is seemingly larger than any of these differences. Therefore, it is assumed that any differences in dissociation rate constants for varying noble gas third-bodies are negligible.

There are two key observations that can be made about the rate constant data shown in Fig. 9. Firstly, while there is considerable scatter in both the discharge-flow tube and shock tube data, the computational data from the master equation studies are generally in good agreement with each other and roughly lie in the middle of the scatter of the experimental data. The only exception to this is the computational data of Furudate et al. Furudate, Fujita, and Abe (2006) (M = $\rm H_{2}$ ), where there is disagreement between even the “corrected” and “uncorrected” versions of their rate constants. These two rate constants roughly bracket the rate constants from the other computational studies.

A second observation is that the temperature dependencies of most of the shock tube data are similar to that of the computational data between 2,000 and 5,000 K. Recall that Gardiner and Kistiakowsky Gardiner and Kistiakowsky (1961), Rink Rink (1962), Patch Patch (1962), and Jacobs et al. Jacobs, Giedt, and Cohen (1967) assumed a value of $n_{\rm r}=-1$ for their reported $k_{\rm r}$ fits, as discussed previously in section IV.1. The other shock tube studies did not however, and report rate constants with similar magnitudes and trends. There are a couple of notable exceptions. For M = $\rm H_{2}$ , Breshears and Bird Breshears and Bird (1972) report a rate constant that increases with temperature, which is inconsistent with all of the other studies. For M = H, Hurle et al. Hurle, Jones, and Rosenfeld (1969) argued that the extreme temperature dependence in their data (and that of Sutton Sutton (1962)), along with the low-temperature, “upper bound” estimate by Bennett and Blackmore Bennett and Blackmore (1968), suggested a non-monotonic trend in the M = H rate constant with a maximum value between approximately 2,000 and 3,000 K. This argument is not supported by the computational results of Martin et al. Martin, Schwarz, and Mandy (1996) and Kim et al. Kim, Kwon, and Park (2009), which predict monotonically decreasing trends with temperature. Shui Shui (1973) suggested that the high-temperature results of Hurle et al. Hurle, Jones, and Rosenfeld (1969) may be in error due to the reliance of their measurements on the assumption that the vibrational mode of $\rm H_{2}$ was in thermal equilibrium with the translational mode. From section III.1, it is clear that for the higher temperatures in the experiments by Hurle et al. Hurle, Jones, and Rosenfeld (1969) (approaching 7,000 K), an assumption of thermal equilibrium would be invalid in the non-recombining QSS limit where $T_{\rm v}$ < $T_{\rm t}$ . Ultimately, this suggests that due to the lack of reliable low-temperature data for M = H, the temperature dependence of the M = H fit must be determined almost entirely from the computational studies of Martin et al. Martin, Schwarz, and Mandy (1996) and Kim et al. Kim, Kwon, and Park (2009).

Based on the results of the review, new rate constant fits are proposed and plotted along with the rate constants of Leibowitz Leibowitz (1973) in Fig. 9. The Leibowitz rate constants are the most widely used in the ice and gas giant entry flow literature Palmer, Prabhu, and Cruden (2014); Higdon et al. (2018); Liu et al. (2021); Hansson et al. (2021); Carroll and Brandis (2023); Carroll et al. (2023); Coelho and da Silva (2023); Rataczak, Boyd, and McMahon (2024), and are based on an extrapolation of the rate constants of Jacobs et al. Jacobs, Giedt, and Cohen (1967). A more detailed discussion of the Leibowitz rate constants is provided in Appendix C. The rate constant fits were chosen to lie in the middle of the scatter of both the low-temperature and high-temperature experimental data where available, while still following the larger temperature dependencies predicted by the computational studies. The proposed rate constant fits are summarized in Table 4. These fits are for $k_{\rm d,nr}/K_{\rm eq}$ instead of $k_{\rm d,nr}$ , as this was found to provide a better fit to the rate constant data. The M = H fit below 450 K and the M = {noble gas} fit above 16,000 K should be used with some caution, as the fits in these regions are extrapolations that go beyond the temperature ranges of the reviewed data.

Table 4: Fits of

k_{\rm d,nr}

K_{\rm eq}

M	$k_{\rm d,nr}/K_{\rm eq}$ $\rm[cm^{6}/mol^{2}/s]$
$\rm H_{2}$	$1.779\times 10^{17}T_{\rm t}^{-0.69}$
H	$5.461\times 10^{17}T_{\rm t}^{-0.61}$
He, Ne, Ar, Xe, Kr	$8.168\times 10^{17}T_{\rm t}^{-0.95}$

Figure 10 shows the same rate constant data as Fig. 9, but now normalized by the fits of the present work. The majority of the reviewed data falls within a factor of two of the proposed fits, and in general, the fits reproduce the reviewed data more accurately than the rate constants of Leibowitz Leibowitz (1973). This is especially true in the lower temperature ranges, where the Leibowitz rate constants severely overpredict the reviewed data. This overprediction is not surprising, as the shock tube study of Jacobs et al. Jacobs, Giedt, and Cohen (1967) (upon which the Leibowitz rate constants are based) only measured rates to a minimum temperature of 2,900 K. The aerothermal heating uncertainty study for Saturn and Uranus entry probes by Palmer et al. Palmer, Prabhu, and Cruden (2014) assumed an uncertainty factor of $\pm$ 1 order of magnitude on the rate constants for $\rm H_{2}$ dissociation with all third-bodies. The more recent uncertainty quantification study for ice giant aerocapture applications by Rataczak et al. Rataczak, Boyd, and McMahon (2024) instead assumed uncertainty bounds of 0.25 and 3 times the baseline rate constants. Both studies used the rate constants of Leibowitz Leibowitz (1973) as a baseline. The results of Fig. 9 highlight that, with the fits of the present work, an uncertainty factor of two is more appropriate than the wider uncertainty bounds used in either of these previous uncertainty quantification studies.

Figure 10 also shows that none of the rate constants from any one of the shock tube studies consistently over or underpredicts the fits for all three third-bodies. For example, in the case of Patch Patch (1962), the reported rate constants overpredict the fits for the M = $\rm H_{2}$ and H cases, but underpredict the fit for the M = {noble gas} case. The only exception are the results of Myerson and Watt Myerson and Watt (1968), but only because they did not report a rate constant for the M = H case. This result suggests that while the total dissociation rates measured by the shock tube studies may be accurate, the reported third-body efficiencies may not be, possibly due to the simultaneous third-body fitting procedure that was used in all of the shock tube studies.

V Conclusion

In this work, a detailed description of the non-equilibrium dissociation of diatomic molecules was presented, then applied to the case of $\rm H_{2}$ dissociation. The master equation formulation was used to analyze general non-equilibrium rate constant expressions in the thermal equilibrium limit, the QSS regime, and the pre-QSS regime. In the QSS regime (up to and including thermal equilibrium), it was found that the chemical source term could be described entirely as a function of the steady dissociation rate constant in the absence of recombination, $k_{\rm d,nr}$ , which in turn is only a function of the translational temperature, $T_{\rm t}$ . To account for pre-QSS effects, a simple modification of the QSS rate constant expression was proposed through the use of an additional rate constant term, $\eta(T_{\rm t})$ . The complete expression, which is ultimately only a function of $T_{\rm t}$ and the fraction of dissociation, $\phi_{\rm H}$ , was able to reproduce the results of detailed master equation simulations of a 0-D isothermal and isochoric reactor for the case of $\rm H_{2}$ dissociation with the third-bodies $\rm H_{2}$ , H, and He. Finally, an extensive review of $k_{\rm d,nr}(T_{\rm t})$ for $\rm H_{2}$ dissociation was performed by considering both experimental and computational data sources. From this review, fits of $k_{\rm d,nr}/K_{\rm eq}$ that could describe the majority of the reviewed rate constants from the literature within a factor of two were proposed for each relevant third-body.

Practically, if pre-QSS effects can be neglected, the use of the QSS source term expression of Eq. (41) along with the fits of $k_{\rm d,nr}/K_{\rm eq}$ from Table 4 are recommended. By implementing the fits of $k_{\rm d,nr}/K_{\rm eq}$ as “recombination” rate constants, the dissociation rate constants, $k_{\rm d,nr}$ , can be computed by applying macroscopic detailed balance, where $K_{\rm eq}$ can be computed readily using the NASA9 polynomial fits of McBride et al. McBride, Zehe, and Gordon (2002). If pre-QSS effects are non-negligible, the same approach can be applied, but instead through the use of the pre-QSS corrected source term (Eq. (78)) and rate constant (Eq. (75) with $\epsilon=10^{-3}$ and the fits of $\eta(T_{\rm t})$ from section III.2) expressions.

There are two caveats that should be considered when using the rate constant expressions proposed in the present work. First, for the simulation of mixtures, the analysis presented is only strictly valid when the dissociation of the diatomic molecule of interest, e.g., $\rm H_{2}$ , occurs in an infinitely dilute bath composed solely of one third-body. However, in most practical CFD simulations and shock tube experiments, dissociation does not occur in an infinitely dilute bath, but instead occurs in a multi-component mixture. If the linear mixture rule is assumed such that the chemical source terms can be computed as just a linear sum over each of the third-bodies in the mixture, then, the same expressions proposed in this work can be used directly. This is discussed in more detail in Appendix A.

Separately, it is not clear if the proposed expressions are sufficient for the accurate simulation of recombination-dominated flows at high temperatures. As highlighted in the recent study by Macdonald Macdonald (2024), recombining flows are characterized by internal state distributions with an overpopulation of excited states (relative to a Boltzmann distribution). This is the opposite of the QSS behavior in dissociation-dominated flows, which is characterized by an underpopulation of excited states Singh and Schwartzentruber (2017, 2020c). As implied by Eq. (11), this difference in rovibrational distributions will lead to different behaviors of the aggregate dissociation rate constant, $k_{\rm d}$ . This means that even if the recombination rate constant is always given by $k_{\rm r}=k_{\rm d,th}/K_{\rm eq}$ , the chemical source term will not be accurately predicted without considering the impact of this difference in dissociation rates. The further analysis of this topic will be the subject of future work.

Acknowledgements.

This work is supported by the NASA entry systems modeling project under grant number 80NSSC21K1751.

Author Declarations

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Alex T. Carroll: Conceptualization (equal), Methodology (equal), Data curation (lead), Formal analysis (lead), Investigation (lead), Software (lead), Validation (lead), Visualization (lead), Writing - original draft (lead), Writing - review and editing (equal). Jacob Wolmer: Data curation (supporting), Investigation (supporting), Visualization (supporting). Guillaume Blanquart: Conceptualization (equal), Methodology (equal), Formal analysis (supporting), Writing - review and editing (equal), Funding acquisition (lead), Supervision (lead). Aaron M. Brandis: Writing - review and editing (supporting), Supervision (supporting). Brett A. Cruden: Writing - review and editing (supporting), Supervision (supporting).

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Appendix A Dissociation in Multi-Component Mixtures

Throughout the present work, only dissociation that occurs in a single-component mixture, i.e., dissociation in an infinitely dilute bath composed solely of one third-body, has been considered. However, for most practical simulations and experiments of interest, it is necessary to consider dissociation that occurs in a multi-component mixture. Considering again the general diatom $\rm A_{2}$ , for a mixture, the source term for A can be written as

\frac{dn_{\rm A}}{dt}=\sum_{\rm M}\left[2n_{\rm M}n_{\rm A_{2}}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}\frac{n_{\rm A_{2}\it(J,\nu)}}{n_{\rm A_{2}}}k_{\rm M}(J,\nu\rightarrow c)-2n_{\rm M}n_{\rm A}^{2}\sum_{\nu=0}^{\nu_{\rm max}}\sum_{J=0}^{J_{\rm max}(\nu)}k_{\rm M}(c\rightarrow J,\nu)\right],

(84)

where $k_{\rm M}(J,\nu\rightarrow c)$ and $k_{\rm M}(c\rightarrow J,\nu)$ are the state-specific dissociation and recombination rate constants with the third-body M. Because the recombination term in Eq. (84) is not a function of the internal state distribution of $\rm A_{2}$ , it is the sum of the recombination terms of all third-bodies considered independently. This is not true for the dissociation term in Eq. (84), as $n_{{\rm A_{2}}(J,\nu)}$ is now determined through collisions with all third-bodies in the mixture. Stated differently, the mixture-dependent $n_{{\rm A_{2}}(J,\nu)}$ in Eq. (84) is not necessarily the same as the distributions obtained in the single-component cases, denoted here as $n_{{\rm A_{2}}(J,\nu),{\rm M}}$ .

The most common method for relating $n_{{\rm A_{2}}(J,\nu)}$ to $n_{{\rm A_{2}}(J,\nu),{\rm M}}$ is through the use of the linear mixture rule. In the linear mixture rule, it is assumed that the source term for A can be computed as

\frac{dn_{\rm A}}{dt}=\sum_{\rm M}\frac{dn_{\rm A}}{dt}\biggr{|}_{\rm M}=\sum_{\rm M}\left[2n_{\rm M}n_{\rm A_{2}}k_{\rm d,M}-2n_{\rm M}n_{\rm A}^{2}k_{\rm r,M}\right].

(85)

Here, $dn_{\rm A}/dt|_{\rm M}$ , $k_{\rm d,M}$ , and $k_{\rm r,M}$ are the source term and aggregate dissociation and recombination rate constants for the single-component cases. To get Eq. (85) from the general expression of Eq. (84), it is sufficient to assume that $n_{{\rm A_{2}}(J,\nu)}$ is related to $n_{{\rm A_{2}}(J,\nu),{\rm M}}$ by

n_{{\rm A_{2}}(J,\nu)}=\frac{\sum_{\rm M}n_{{\rm A_{2}}(J,\nu){\rm,M}}n_{\rm M}k_{\rm M}(J,\nu\rightarrow c)}{\sum_{\rm M}n_{\rm M}k_{\rm M}(J,\nu\rightarrow c)}.

(86)

In other words, $n_{{\rm A_{2}}(J,\nu)}$ is a linear combination of $n_{{\rm A_{2}}(J,\nu),{\rm M}}$ weighted by the state-specific branching ratios.

As discussed by Boyd Boyd (1977), Dove et al. Dove, Halperin, and Raynor (1984), and Carruthers and Teitelbaum Carruthers and Teitelbaum (1988), the range of validity of the linear mixture rule is an open question. Unfortunately, to the authors’ knowledge, there are currently no detailed/ master equation simulation data available for $\rm H_{2}$ / H/ He mixtures. Therefore, the recommended practice remains to use Eq. (85) for the simulation of mixtures.

Appendix B Recombination as a Second or Third-Order Process

Throughout this work, it has been implicitly assumed that recombination is a third-order process. This assumption has been made in the master equations and the subsequent QSS formulation presented in sections II and III, as well as in the interpretation of rate constant data from the literature for $\rm H_{2}$ dissociation throughout section IV. In particular, it has been assumed that recombination occurs as the direct inverse process of dissociation, and hence that recombination rates can be characterized entirely by applying the micro-reversibility relation of Eq. (8).

Under certain pressure and temperature conditions however, recombination can occur primarily as a second-order process. This can be seen by considering recombination that occurs through a series of two-body collisions, as opposed to a “direct” three-body collision. As reviewed by Mirahmadi and Pérez-Ríos Mirahmadi and Pérez-Ríos (2022), this “indirect” approach treats recombination in two steps. First, two bodies collide to form an intermediate complex,

\rm A+A\leftrightarrow A_{2}^{*}.

(87)

Then, this intermediate complex is stabilized into a bound state/ recombined molecule via a collision with a third-body,

\rm A_{2}^{*}+M\leftrightarrow A_{2}+M.

(88)

Here, $\rm A_{2}$ is a generic diatom as before, and $\rm A_{2}^{*}$ is an intermediate complex/ quasibound state. These quasibound states are bound classically, but have a finite lifetime quantum mechanically due to tunneling. If the backwards direction of reaction (87) is dominant, recombination will occur primarily as a third-order process. Conversely, if the forwards direction of reaction (88) is dominant, recombination will occur primarily as a second-order process.

In detailed chemical models, the pressure dependence of a recombination reaction rate is often modeled with the Lindemann form Lindemann et al. (1922) or the more general Troe form Gilbert, Luther, and Troe (1983). In the Lindemann expression, low-pressure and high-pressure rate constants, $k_{\rm 0}$ and $k_{\infty}$ respectively, are first defined. Then, the rate constant at any intermediate pressure is computed as

k=\frac{k_{0}n_{\rm M}}{P_{\rm r}+1},

(89)

where

P_{\rm r}\equiv\frac{k_{0}n_{\rm M}}{k_{\infty}}

(90)

is the reduced pressure. For an infinitely dilute mixture with the third-body, M, $n_{\rm M}=P/(k_{\rm B}T_{\rm t})$ . For the recombination of a diatom, the high-pressure and low-pressure limits correspond to the limits dominated by two-body and three-body collisions, respectively. For $\rm H_{2}$ , $k_{\rm 0}$ can be estimated by using the rate constant fits from section IV.4. $k_{\infty}$ on the other hand can be estimated by using the upper-bound value of the bimolecular collision rate for homonuclear collisions,

k_{\infty}\approx\frac{1}{2}\sigma_{\rm H_{2}}\sqrt{\frac{8k_{\rm B}T_{\rm t}}{\pi\mu}},

(91)

where $\sigma_{\rm H_{2}}=\pi d_{\rm H}^{2}$ is the hard-sphere collision cross section with $d_{\rm H}\approx 3\rm\mathring{A}$ Jasper and Miller (2014), and $\mu=m_{\rm H}/2=8.37\times 10^{-28}$ kg is the reduced mass. Of all of the discharge-flow tube studies considered in section IV.2, the experiments by Lynch et al. Lynch, Schwab, and Michael (1976) at $T_{\rm t}=$ 298 K and $P=2$ atm exhibit the largest values of $P_{\rm r}$ and hence the furthest departures from the low-pressure limit. However, even at this condition, $P_{\rm r}$ as estimated by Eq. (90) is still approximately $10^{-3}$ . Therefore, from Eq. (89), $k\approx k_{0}n_{\rm M}$ , and the choice to interpret the rate constants from the discharge-flow tube studies as third-order is appropriate.

Outside of the discharge-flow tube studies, there have been a number of computational studies Pack, Snow, and Smith (1972); Whitlock, Muckerman, and Roberts (1972, 1974); Orel (1987); Schwenke (1988); Esposito and Capitelli (2009); Forrey (2013) that have directly computed the third-order recombination rate constant of $\rm H_{2}$ by using the resonance complex theory of Roberts et al. Roberts, Bernstein, and Curtiss (1969). In this theory, it is assumed that recombination occurs exclusively through the indirect pathway. Then, a small number of select quasibound states are considered to be in equilibrium with the continuum (i.e., reaction (87)), and the downward flux of these quasibound states (i.e., reaction (88)) are assumed to be predominantly responsible for the recombination rate. Unfortunately, the accuracy of the phenomenological PESs used in several of these studies is questionable, and as pointed out by Schwenke Schwenke (1988, 1990), the reported rate constants from these studies are not always consistent with the high-temperature experimental and computational rate constants reviewed in sections IV.1 and IV.3, respectively. For these reasons, the studies based on the resonance complex theory are not considered for the review in the present work.

Appendix C Rate Constants of Leibowitz

Three different sets of $\rm H_{2}$ dissociation rate constants have been proposed by Leibowitz, namely the models in Leibowitz Leibowitz (1973), Leibowitz et al. Leibowitz, Menard, and Stickford (1973), and Leibowitz and Kuo Leibowitz and Kuo (1976). The reported rate constants are summarized below in Table 5. In all three versions, the rate constants for the third-bodies M = ${\rm H_{2}}$ and H are evaluated from the M = He case using a constant multiplicative factor.

Table 5: Rate Constants of Leibowitz in

\rm cm^{3}/mol/s

	M = He	M = $\rm H_{2}$	M = H
Leibowitz, 1973 Leibowitz (1973)	$4.17\times 10^{18}T_{\rm t}^{-1}\exp(-51952/T_{\rm t})$ ¹¹1The term inside the exponential is mistakenly reported as a positive value in Leibowitz Leibowitz (1973).	$2.5k_{\rm d,M=He}$	$20.0k_{\rm d,M=He}$
Leibowitz et al. 1973 Leibowitz, Menard, and Stickford (1973)	$1.1\times 10^{18}T_{\rm t}^{-1}[1-\exp(-1.5\times 10^{8}/T_{\rm t}^{2})]\exp(-52340/T_{\rm t})$ ²²2The leading $T_{\rm t}^{-1}$ term is mistakenly reported as just $T_{\rm t}$ in both Leibowitz et al. Leibowitz, Menard, and Stickford (1973) and Leibowitz and Kuo Leibowitz and Kuo (1976).	$2.5k_{\rm d,M=He}$	$14.0k_{\rm d,M=He}$
Leibowitz and Kuo, 1976 Leibowitz and Kuo (1976)	$4.3\times 10^{18}T_{\rm t}^{-1}[1-\exp(-1.5\times 10^{8}/T_{\rm t}^{2})]\exp(-52340/T_{\rm t})$ ²²2The leading $T_{\rm t}^{-1}$ term is mistakenly reported as just $T_{\rm t}$ in both Leibowitz et al. Leibowitz, Menard, and Stickford (1973) and Leibowitz and Kuo Leibowitz and Kuo (1976).³³3The term $[1-\exp(-1.5\times 10^{8}/T_{\rm t}^{2})]$ is mistakenly reported as $[1-\exp(1.5\times 10^{8}/T_{\rm t}^{2})]$ in Leibowitz and Kuo Leibowitz and Kuo (1976).	$2.5k_{\rm d,M=He}$	$14.0k_{\rm d,M=He}$

It is important to note that the three different sets of rate constants are not based on the same source of rate constant data. For the first set, Leibowitz Leibowitz (1973) cites the shock tube study of Jacobs et al. Jacobs, Giedt, and Cohen (1967). For the second set, Leibowitz et al. Leibowitz, Menard, and Stickford (1973) cites the shock tube study of Hurle et al. Hurle, Jones, and Rosenfeld (1969), adding that the rate constants were selected to be “consistent with shock tube measurements below 7,000 K, reach a maximum at 20,000 K, and then decay with increasing temperature”. Finally, for the third set by Leibowitz and Kuo Leibowitz and Kuo (1976), the same temperature dependence and third-body efficiencies are used as in the second set by Leibowitz et al. Leibowitz, Menard, and Stickford (1973), but the rate constants are all increased by a factor of four. Leibowitz and Kuo Leibowitz and Kuo (1976) do not provide an explanation for this change, and simply cite the two previous sets of rate constants by Leibowitz Leibowitz (1973) and Leibowitz and Kuo Leibowitz, Menard, and Stickford (1973) as the source for these rate constants.

Figure 11 shows the rate constants (plotted as $k_{\rm d,nr}/K_{\rm eq}$ ) for all three formulations compared to the fits from the shock tube studies of Jacobs et al. Jacobs, Giedt, and Cohen (1967) and Hurle et al. Hurle, Jones, and Rosenfeld (1969). While the rate constants of Leibowitz Leibowitz (1973) and Leibowitz and Kuo Leibowitz and Kuo (1976) are similar to each other and to the rate constants of Jacobs et al. Jacobs, Giedt, and Cohen (1967) for the temperature range between approximately 3,000 and 5,000 K, the rate constants of Leibowitz et al. Leibowitz, Menard, and Stickford (1973) are consistently lower for all three third-bodies. Additionally, despite citing the shock tube study by Hurle et al. Hurle, Jones, and Rosenfeld (1969), the rate constants of Leibowitz et al. Leibowitz and Kuo (1976) do not follow the rate constants proposed by Hurle et al. Hurle, Jones, and Rosenfeld (1969).

Although Fig. 11 highlights that there are apparent differences between the three sets of rate constants proposed by Leibowitz Leibowitz (1973); Leibowitz, Menard, and Stickford (1973); Leibowitz and Kuo (1976), many of the studies in the ice and gas giant entry flow literature cite the three different formulations interchangeably. To the authors’ knowledge, all of the studies that have implemented and used the rate constants attributed to Leibowitz Palmer, Prabhu, and Cruden (2014); Higdon et al. (2018); Liu et al. (2021); Hansson et al. (2021); Carroll and Brandis (2023); Carroll et al. (2023); Coelho and da Silva (2023); Rataczak, Boyd, and McMahon (2024) have used the rate constants from the first formulation of Leibowitz Leibowitz (1973), as they follow the modified Arrhenius form and thus can be readily implemented in most CFD codes. For this reason, comparisons are only made to this first set of rate constants from Leibowitz Leibowitz (1973) in section IV.4.

References

Draine, Roberge, and Dalgarno (1983) B. T. Draine, W. G. Roberge, and A. Dalgarno, “Magnetohydrodynamic shock waves in molecular clouds,” The Astrophysical Journal 264, 485 (1983).
Hollenbach and McKee (1989) D. Hollenbach and C. F. McKee, “Molecule formation and infrared emission in fast interstellar shocks. III - results for J shocks in molecular clouds,” The Astrophysical Journal 342, 306 (1989).
Chang and Martin (1991) C. A. Chang and P. G. Martin, “Partially dissociative jump shocks in molecular hydrogen,” The Astrophysical Journal 378, 202 (1991).
Bell and Cowan (2018) T. J. Bell and N. B. Cowan, “Increased heat transport in ultra-hot Jupiter atmospheres through $\rm H_{2}$ dissociation and recombination,” The Astrophysical Journal Letters 857, L20 (2018).
Kristensen et al. (2023) L. E. Kristensen, B. Godard, P. Guillard, A. Gusdorf, and G. Pineau des Forêts, “Shock excitation of $\rm H_{2}$ in the James Webb space telescope era,” Astronomy & Astrophysics 675, A86 (2023).
Matveyev and Silakov (1995) A. A. Matveyev and V. P. Silakov, “Kinetic processes in a highly-ionized non-equilibrium hydrogen plasma,” Plasma Sources Science and Technology 4, 606–617 (1995).
Silakov et al. (1996) V. P. Silakov, A. A. Matveyev, A. V. Chebotarev, and D. K. Otorbaev, “Non-equilibrium properties of a flowing hydrogen cascaded arc plasma: Kinetic modelling,” Journal of Physics D: Applied Physics 29, 2111–2118 (1996).
Capitelli et al. (2002) M. Capitelli, R. Celiberto, F. Esposito, A. Laricchiuta, K. Hassouni, and S. Longo, “Elementary processes and kinetics of $\rm H_{2}$ plasmas for different technological applications,” Plasma Sources Science and Technology 11, A7–A25 (2002).
Shakhatov and Lebedev (2018) V. A. Shakhatov and Y. A. Lebedev, “Kinetics of populations of singlet and triplet states in non-equilibrium hydrogen plasma,” Journal of Physics D: Applied Physics 51, 213001 (2018).
Smith et al. (2024) G. J. Smith, P. Diomede, A. R. Gibson, S. J. Doyle, V. Guerra, M. J. Kushner, T. Gans, and J. P. Dedrick, “Low-pressure inductively coupled plasmas in hydrogen: Impact of gas heating on the spatial distribution of atomic hydrogen and vibrationally excited states,” Plasma Sources Science and Technology 33, 025002 (2024).
Taylor et al. (2013) B. Taylor, D. Kessler, V. Gamezo, and E. Oran, “Numerical simulations of hydrogen detonations with detailed chemical kinetics,” Proceedings of the Combustion Institute 34, 2009–2016 (2013).
Voelkel et al. (2016) S. Voelkel, D. Masselot, P. L. Varghese, and V. Raman, “Analysis of hydrogen-air detonation waves with vibrational nonequilibrium,” in AIP Conference Proceedings (2016).
Shi et al. (2016) L. Shi, H. Shen, P. Zhang, D. Zhang, and C. Wen, “Assessment of vibrational non-equilibrium effect on detonation cell size,” Combustion Science and Technology 189, 841–853 (2016).
Vargas et al. (2022) J. Vargas, R. Mével, M. Lino da Silva, and D. A. Lacoste, “Development of a steady detonation reactor with state-to-state thermochemical modeling,” Shock Waves 32, 679–689 (2022).
Baumgart, Yao, and Blanquart (2025) A. Baumgart, M. X. Yao, and G. Blanquart, “Tabulated chemistry approach for detonation simulations,” Combustion and Flame 272, 113878 (2025).
Park (2011) C. Park, “Nonequilibrium chemistry and radiation for Neptune entry,” Journal of Spacecraft and Rockets 48, 897–903 (2011).
Park (2012) C. Park, “Nonequilibrium ionization and radiation in hydrogen-helium mixtures,” Journal of Thermophysics and Heat Transfer 26, 231–243 (2012).
Palmer, Prabhu, and Cruden (2014) G. Palmer, D. Prabhu, and B. A. Cruden, “Aeroheating uncertainties in Uranus and Saturn entries by the Monte Carlo method,” Journal of Spacecraft and Rockets 51, 801–814 (2014).
Erb, West, and Johnston (2019) A. J. Erb, T. K. West, and C. O. Johnston, “Investigation of Galileo probe entry heating with coupled radiation and ablation,” in AIAA Aviation 2019 Forum (American Institute of Aeronautics and Astronautics, 2019).
Carroll and Brandis (2023) A. T. Carroll and A. M. Brandis, “Stagnation point convective heating correlations for entry into $\rm H_{2}$ /He atmospheres,” in AIAA Scitech 2023 Forum (American Institute of Aeronautics and Astronautics, 2023).
Scoggins, Hinkle, and Shellabarger (2024) J. B. Scoggins, A. D. Hinkle, and E. Shellabarger, “Aeroheating environment of aerocapture systems for Uranus orbiters,” in AIAA Scitech 2024 Forum (American Institute of Aeronautics and Astronautics, 2024).
Steer et al. (2024) J. Steer, P. Collen, A. Glenn, C. Hambidge, L. J. Doherty, M. McGilvray, T. Sopek, S. Loehle, and L. Walpot, “Commissioning of upgrades to T6 to study giant planet entry,” Journal of Spacecraft and Rockets , 1–18 (2024).
Steuer et al. (2024) D. C. Steuer, N. Donaldson, C. A. Dürnhofer, D. Leiser, S. Fasoulas, and S. Loehle, “Investigation on the usability of arc-jet generators for the experimental simulation of giant planet entry,” in AIAA Aviation FORUM AND ASCEND 2024 (American Institute of Aeronautics and Astronautics, 2024).
Rataczak, Boyd, and McMahon (2024) J. A. Rataczak, I. D. Boyd, and J. W. McMahon, “Parametric sensitivity analysis and uncertainty quantification of ice giant aerocapture aerothermodynamics,” in AIAA Aviation FORUM AND ASCEND 2024 (2024).
Marrone and Treanor (1963) P. V. Marrone and C. E. Treanor, “Chemical relaxation with preferential dissociation from excited vibrational levels,” The Physics of Fluids 6, 1215–1221 (1963).
Macheret et al. (1994) S. O. Macheret, A. A. Fridman, I. V. Adamovich, J. W. Rich, and C. E. Treanor, “Mechanisms of nonequilibrium dissociation of diatomic molecules,” in 6th Joint Thermophysics and Heat Transfer Conference (American Institute of Aeronautics and Astronautics, 1994).
Porter (1974) R. N. Porter, “Molecular trajectory calculations,” Annual Review of Physical Chemistry 25, 317–355 (1974).
Koura (1997) K. Koura, “Monte Carlo direct simulation of rotational relaxation of diatomic molecules using classical trajectory calculations: Nitrogen shock wave,” Physics of Fluids 9, 3543–3549 (1997).
Kim, Kwon, and Park (2009) J. G. Kim, O. J. Kwon, and C. Park, “Master equation study and nonequilibrium chemical reactions for H+ $\rm H_{2}$ and He+ $\rm H_{2}$ ,” Journal of Thermophysics and Heat Transfer 23, 443–453 (2009).
Kim and Boyd (2012) J. G. Kim and I. D. Boyd, “State-resolved thermochemical nonequilibrium analysis of hydrogen mixture flows,” Physics of Fluids 24 (2012), 10.1063/1.4747340.
Panesi et al. (2013) M. Panesi, R. L. Jaffe, D. W. Schwenke, and T. E. Magin, “Rovibrational internal energy transfer and dissociation of $\mathbf{\rm N_{2}(^{1}\Sigma_{g}^{+})-\rm N(^{4}S_{u})}$ system in hypersonic flows,” The Journal of Chemical Physics 138 (2013), 10.1063/1.4774412.
Andrienko and Boyd (2015) D. Andrienko and I. D. Boyd, “Investigation of oxygen vibrational relaxation by quasi-classical trajectory method,” Chemical Physics 459, 1–13 (2015).
Norman, Valentini, and Schwartzentruber (2013) P. Norman, P. Valentini, and T. Schwartzentruber, “GPU-accelerated classical trajectory calculation direct simulation Monte Carlo applied to shock waves,” Journal of Computational Physics 247, 153–167 (2013).
Valentini et al. (2014) P. Valentini, P. Norman, C. Zhang, and T. E. Schwartzentruber, “Rovibrational coupling in molecular nitrogen at high temperature: An atomic-level study,” Physics of Fluids 26 (2014), 10.1063/1.4875279.
Schwartzentruber, Grover, and Valentini (2018) T. E. Schwartzentruber, M. S. Grover, and P. Valentini, “Direct molecular simulation of nonequilibrium dilute gases,” Journal of Thermophysics and Heat Transfer 32, 892–903 (2018).
Park (1988) C. Park, “Assessment of a two-temperature kinetic model for dissociating and weakly ionizing nitrogen,” Journal of Thermophysics and Heat Transfer 2, 8–16 (1988).
Park (1989a) C. Park, “Assessment of two-temperature kinetic model for ionizing air,” Journal of Thermophysics and Heat Transfer 3, 233–244 (1989a).
Park (1989b) C. H. Park, Nonequilibrium Hypersonic Aerothermodynamics (Wiley, 1989).
Chaudhry and Candler (2019) R. S. Chaudhry and G. V. Candler, “Statistical analyses of quasiclassical trajectory data for air dissociation,” in AIAA Scitech 2019 Forum (American Institute of Aeronautics and Astronautics, 2019).
Chaudhry et al. (2020) R. S. Chaudhry, I. D. Boyd, E. Torres, T. E. Schwartzentruber, and G. V. Candler, “Implementation of a chemical kinetics model for hypersonic flows in air for high-performance CFD,” in AIAA Scitech 2020 Forum (American Institute of Aeronautics and Astronautics, 2020).
Singh and Schwartzentruber (2020a) N. Singh and T. Schwartzentruber, “Consistent kinetic–continuum dissociation model I. kinetic formulation,” The Journal of Chemical Physics 152 (2020a), 10.1063/1.5142752.
Singh and Schwartzentruber (2020b) N. Singh and T. Schwartzentruber, “Consistent kinetic-continuum dissociation model. II. continuum formulation and verification,” The Journal of Chemical Physics 152 (2020b), 10.1063/1.5142754.
Magin et al. (2012) T. E. Magin, M. Panesi, A. Bourdon, R. L. Jaffe, and D. W. Schwenke, “Coarse-grain model for internal energy excitation and dissociation of molecular nitrogen,” Chemical Physics 398, 90–95 (2012).
Sahai et al. (2017) A. Sahai, B. Lopez, C. O. Johnston, and M. Panesi, “Adaptive coarse graining method for energy transfer and dissociation kinetics of polyatomic species,” The Journal of Chemical Physics 147 (2017), 10.1063/1.4996654.
Venturi et al. (2020) S. Venturi, M. P. Sharma, B. Lopez, and M. Panesi, “Data-inspired and physics-driven model reduction for dissociation: Application to the $\rm O_{2}$ + O system,” The Journal of Physical Chemistry A 124, 8359–8372 (2020).
Liu et al. (2015) Y. Liu, M. Panesi, A. Sahai, and M. Vinokur, “General multi-group macroscopic modeling for thermo-chemical non-equilibrium gas mixtures,” The Journal of Chemical Physics 142 (2015), 10.1063/1.4915926.
Park (1993) C. Park, “Review of chemical-kinetic problems of future NASA missions. I - Earth entries,” Journal of Thermophysics and Heat Transfer 7, 385–398 (1993).
Park et al. (1994) C. Park, J. T. Howe, R. L. Jaffe, and G. V. Candler, “Review of chemical-kinetic problems of future NASA missions. II - Mars entries,” Journal of Thermophysics and Heat Transfer 8, 9–23 (1994).
Park, Jaffe, and Partridge (2001) C. Park, R. L. Jaffe, and H. Partridge, “Chemical-kinetic parameters of hyperbolic Earth entry,” Journal of Thermophysics and Heat Transfer 15, 76–90 (2001).
Sarma (2000) G. Sarma, “Physico–chemical modelling in hypersonic flow simulation,” Progress in Aerospace Sciences 36, 281–349 (2000).
Pietanza et al. (2021) L. D. Pietanza, O. Guaitella, V. Aquilanti, I. Armenise, A. Bogaerts, M. Capitelli, G. Colonna, V. Guerra, R. Engeln, E. Kustova, A. Lombardi, F. Palazzetti, and T. Silva, “Advances in non-equilibrium $\rm CO_{2}$ plasma kinetics: A theoretical and experimental review,” The European Physical Journal D 75 (2021), 10.1140/epjd/s10053-021-00226-0.
National Academies of Sciences, Engineering, and Medicine (2022) National Academies of Sciences, Engineering, and Medicine, Origins, Worlds, and Life: A Decadal Strategy for Planetary Science and Astrobiology 2023-2032 (The National Academies Press, Washington, DC, 2022).
Higdon et al. (2018) K. J. Higdon, B. A. Cruden, A. M. Brandis, D. S. Liechty, D. B. Goldstein, and P. L. Varghese, “Direct Simulation Monte Carlo shock simulation of Saturn entry probe conditions,” Journal of Thermophysics and Heat Transfer 32, 680–690 (2018).
Liu et al. (2021) Y. Liu, C. M. James, R. G. Morgan, P. Jacobs, R. Gollan, and T. McIntyre, “Electron number density measurements in a Saturn entry condition,” in AIAA Scitech 2021 Forum (American Institute of Aeronautics and Astronautics, 2021).
Hansson et al. (2021) K. Hansson, A. T. Carroll, S. J. Poovathingal, and I. D. Boyd, “Analysis of chemical kinetic parameters for hydrogen atmospheres,” in AIAA Scitech 2021 Forum (American Institute of Aeronautics and Astronautics, 2021).
Carroll et al. (2023) A. T. Carroll, G. Blanquart, A. M. Brandis, and B. A. Cruden, “Kinetic and transport modeling for entry flows in hydrogen-helium atmospheres,” in AIAA Aviation 2023 Forum (American Institute of Aeronautics and Astronautics, 2023).
Coelho and da Silva (2023) J. Coelho and M. L. da Silva, “Aerothermodynamic analysis of Neptune ballistic entry and aerocapture flows,” Advances in Space Research 71, 3408–3432 (2023).
Leibowitz (1973) L. P. Leibowitz, “Measurements of the structure of an ionizing shock wave in a hydrogen-helium mixture,” The Physics of Fluids 16, 59–68 (1973).
Jacobs, Giedt, and Cohen (1967) T. A. Jacobs, R. R. Giedt, and N. Cohen, “Kinetics of hydrogen halides in shock waves. II. a new measurement of the hydrogen dissociation rate,” The Journal of Chemical Physics 47, 54–57 (1967).
Cohen and Westberg (1983) N. Cohen and K. R. Westberg, “Chemical kinetic data sheets for high-temperature chemical reactions,” Journal of Physical and Chemical Reference Data 12, 531–590 (1983).
Baulch (1972) D. Baulch, Evaluated Kinetic Data for High Temperature Reactions (Butterworths, 1972).
Baulch et al. (2005) D. L. Baulch, C. T. Bowman, C. J. Cobos, R. A. Cox, T. Just, J. A. Kerr, M. J. Pilling, D. Stocker, J. Troe, W. Tsang, R. W. Walker, and J. Warnatz, “Evaluated kinetic data for combustion modeling: Supplement II,” Journal of Physical and Chemical Reference Data 34, 757–1397 (2005).
Colonna, D’Angola, and Capitelli (2012) G. Colonna, A. D’Angola, and M. Capitelli, “Statistical thermodynamic description of $\rm H_{2}$ molecules in normal ortho/para mixture,” International Journal of Hydrogen Energy 37, 9656–9668 (2012).
Popovas and Jørgensen (2016) A. Popovas and U. G. Jørgensen, “Partition functions - I. improved partition functions and thermodynamic quantities for normal, equilibrium, and ortho and para molecular hydrogen,” Astronomy & Astrophysics 595, A130 (2016).
Kim and Boyd (2013) J. G. Kim and I. D. Boyd, “State-resolved master equation analysis of thermochemical nonequilibrium of nitrogen,” Chemical Physics 415, 237–246 (2013).
Kim, Kwon, and Park (2010) J. G. Kim, O. J. Kwon, and C. Park, “Master equation study and nonequilibrium chemical reactions for hydrogen molecule,” Journal of Thermophysics and Heat Transfer 24, 281–290 (2010).
Haug, Truhlar, and Blais (1987) K. Haug, D. G. Truhlar, and N. C. Blais, “Monte Carlo trajectory and master equation simulation of the nonequilibrium dissociation rate coefficient for Ar+ $\rm H_{2}$ $\rightarrow$ Ar+2H at 4500 K,” The Journal of Chemical Physics 86, 2697–2716 (1987).
Kim (2015) J. G. Kim, “Rovibrational energy transitions and coupled chemical reaction modeling of H+ $\rm H_{2}$ and He+ $\rm H_{2}$ in DSMC,” International Journal of Aeronautical and Space Sciences 16, 347–359 (2015).
Note (1) Number densities, rovibrational temperatures, and rate constants for these studies were only reported as plots. Therefore, these values have all been obtained via digitizations.
Note (2) For the case of M= $\rm H_{2}$ , a modified equilibrium constant, $K_{\rm eq,mod}=2K_{\rm eq}$ , was required to reproduce the equilibrium number densities plotted by Kim and Boyd Kim and Boyd (2012). Additionally, for the case of M=H, $k_{\rm d}$ values were extracted from the number density profiles assuming $n_{\rm M}$ was fixed and equal to $n_{\rm H,0}$ .
Vargas, Monge-Palacios, and Lacoste (2024) J. Vargas, M. Monge-Palacios, and D. A. Lacoste, “State-specific dissociation and inelastic rate constants for collisions of $\rm H_{2}$ with H and He,” Journal of Thermophysics and Heat Transfer 38, 210–221 (2024).
Munafò, Panesi, and Magin (2014) A. Munafò, M. Panesi, and T. E. Magin, “Boltzmann rovibrational collisional coarse-grained model for internal energy excitation and dissociation in hypersonic flows,” Physical Review E 89 (2014), 10.1103/physreve.89.023001.
Sahai et al. (2019) A. Sahai, B. Lopez, C. O. Johnston, and M. Panesi, “Reduced-order modeling of non-equilibrium kinetics and radiation for $\rm CO_{2}$ axisymmetric wake flows,” in 31ST INTERNATIONAL SYMPOSIUM ON RAREFIED GAS DYNAMICS: RGD31, Vol. 2132 (AIP Publishing, 2019) p. 130002.
Gardiner and Kistiakowsky (1961) W. C. Gardiner and G. B. Kistiakowsky, “Thermal dissociation rate of hydrogen,” The Journal of Chemical Physics 35, 1765–1770 (1961).
Sutton (1962) E. A. Sutton, “Measurement of the dissociation rates of hydrogen and deuterium,” The Journal of Chemical Physics 36, 2923–2931 (1962).
Rink (1962) J. P. Rink, “Shock tube determination of dissociation rates of hydrogen,” The Journal of Chemical Physics 36, 262–265 (1962).
Patch (1962) R. W. Patch, “Shock-tube measurement of dissociation rates of hydrogen,” The Journal of Chemical Physics 36, 1919–1924 (1962).
Myerson and Watt (1968) A. L. Myerson and W. S. Watt, “Atom-formation rates behind shock waves in hydrogen and the effect of added oxygen,” The Journal of Chemical Physics 49, 425–433 (1968).
Hurle, Jones, and Rosenfeld (1969) I. R. Hurle, A. Jones, and J. L. J. Rosenfeld, “Shock-wave observations of rate constants for atomic hydrogen recombination from 2500 to 7000 K: Collisional stabilization by exchange of hydrogen atoms,” Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 310, 253–276 (1969).
Breshears and Bird (1972) W. Breshears and P. Bird, “Precise measurements of diatomic dissociation rates in shock waves,” Symposium (International) on Combustion (1972), 10.2172/4653639.
Gordon and McBride (1996) S. Gordon and B. J. McBride, “Computer program for calculation of complex chemical equilibrium compositions and applications,” Tech. Rep. (NASA Reference Publication 1311, 1996).
Lynch, Schwab, and Michael (1976) K. P. Lynch, T. C. Schwab, and J. V. Michael, “Lyman- $\alpha$ absorption photometry at high pressure and atom density kinetic results for H recombination,” International Journal of Chemical Kinetics 8, 651–671 (1976).
Larkin (1968) F. S. Larkin, “Homogeneous rate of recombination of hydrogen atoms,” Canadian Journal of Chemistry 46, 1005–1015 (1968).
Bennett and Blackmore (1968) J. E. Bennett and D. R. Blackmore, “The measurement of the rate of recombination of hydrogen atoms at room temperature by means of e. s. r. spectroscopy,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 305, 553–574 (1968).
Bennett and Blackmore (1970) J. E. Bennett and D. R. Blackmore, “Gas-phase recombination rates of hydrogen and deuterium atoms,” The Journal of Chemical Physics 53, 4400–4401 (1970).
Bennett and Blackmore (1971) J. Bennett and D. Blackmore, “Rates of gas-phase hydrogen-atom recombination at room temperature in the presence of added gases,” Symposium (International) on Combustion 13, 51–59 (1971).
Trainor, Ham, and Kaufman (1973) D. W. Trainor, D. O. Ham, and F. Kaufman, “Gas phase recombination of hydrogen and deuterium atoms,” The Journal of Chemical Physics 58, 4599–4609 (1973).
Mitchell and LeRoy (1977) D. N. Mitchell and D. J. LeRoy, “An experimental test of the orbiting resonance theory of hydrogen atom recombination at room temperature,” The Journal of Chemical Physics 67, 1042–1050 (1977).
Schwenke (1990) D. W. Schwenke, “A theoretical prediction of hydrogen molecule dissociation-recombination rates including an accurate treatment of internal state nonequilibrium effects,” The Journal of Chemical Physics 92, 7267–7282 (1990).
Schwenke (1988) D. W. Schwenke, “A new method for the direct calculation of resonance parameters with application to the quasibound states of the $\rm H_{2}(X^{1}\Sigma_{g}^{+})$ system,” Theoretica Chimica Acta 74, 381–402 (1988).
Martin, Schwarz, and Mandy (1996) P. G. Martin, D. H. Schwarz, and M. E. Mandy, “Master equation studies of the collisional excitation and dissociation of $\rm H_{2}$ molecules by H atoms,” The Astrophysical Journal 461, 265 (1996).
Liu (1973) B. Liu, “Ab initio potential energy surface for linear $\rm H_{3}$ ,” The Journal of Chemical Physics 58, 1925–1937 (1973).
Siegbahn and Liu (1978) P. Siegbahn and B. Liu, “An accurate three-dimensional potential energy surface for $\rm H_{3}$ ,” The Journal of Chemical Physics 68, 2457–2465 (1978).
Truhlar and Horowitz (1978) D. G. Truhlar and C. J. Horowitz, “Functional representation of Liu and Siegbahn’s accurate ab-initio potential energy calculations for H+ $\rm H_{2}$ ,” The Journal of Chemical Physics 68, 2466–2476 (1978).
Truhlar and Horowitz (1979) D. G. Truhlar and C. J. Horowitz, “Erratum: Functional representation of Liu and Siegbahn’s accurate ab-initio potential energy calculations for H+ $\rm H_{2}$ ,” The Journal of Chemical Physics 71, 1514–1514 (1979).
Furudate, Fujita, and Abe (2006) M. Furudate, K. Fujita, and T. Abe, “Coupled rotational-vibrational relaxation of molecular hydrogen at high temperatures,” Journal of Thermophysics and Heat Transfer 20, 457–464 (2006).
Boothroyd et al. (1996) A. I. Boothroyd, W. J. Keogh, P. G. Martin, and M. R. Peterson, “A refined $\rm H_{3}$ potential energy surface,” The Journal of Chemical Physics 104, 7139–7152 (1996).
Boothroyd, Martin, and Peterson (2003) A. I. Boothroyd, P. G. Martin, and M. R. Peterson, “Accurate analytic He– $\rm H_{2}$ potential energy surface from a greatly expanded set of ab initio energies,” The Journal of Chemical Physics 119, 3187–3207 (2003).
Sharma (1994) S. P. Sharma, “Rotational relaxation of molecular hydrogen at moderate temperatures,” Journal of Thermophysics and Heat Transfer 8, 35–39 (1994).
Shui, Appleton, and Keck (1971) V. H. Shui, J. P. Appleton, and J. C. Keck, “The three-body recombination and dissociation of diatomic molecules: A comparison between theory and experiment,” Symposium (International) on Combustion 13, 21–35 (1971).
Shui and Appleton (1971) V. H. Shui and J. P. Appleton, “Gas-phase recombination of hydrogen: A comparison between theory and experiment,” The Journal of Chemical Physics 55, 3126–3132 (1971).
Shui (1973) V. H. Shui, “Thermal dissociation and recombination of hydrogen according to the reactions $\rm H_{2}$ +H $\leftrightarrow$ H+H+H,” The Journal of Chemical Physics 58, 4868–4879 (1973).
Roberge and Dalgarno (1982) W. Roberge and A. Dalgarno, “Collision-induced dissociation of $\rm H_{2}$ and CO molecules,” The Astrophysical Journal 255, 176 (1982).
Esposito, Gorse, and Capitelli (1999) F. Esposito, C. Gorse, and M. Capitelli, “Quasi-classical dynamics calculations and state-selected rate coefficients for H+ $\rm H_{2}$ (v, j) $\rightarrow$ 3H processes: Application to the global dissociation rate under thermal conditions,” Chemical Physics Letters 303, 636–640 (1999).
Dove and Jones (1972) J. Dove and D. Jones, “Effect of rotation on the computed thermal dissociation rate of $\rm H_{2}$ ,” Chemical Physics Letters 17, 134–136 (1972).
Lepp and Shull (1983) S. Lepp and J. M. Shull, “The kinetic theory of $\rm H_{2}$ dissociation,” The Astrophysical Journal 270, 578 (1983).
Blais and Truhlar (1979) N. C. Blais and D. G. Truhlar, “Monte Carlo trajectory study of Ar+ $\rm H_{2}$ : Vibrational selectivity of dissociative collisions at 4500 K and the characteristics of dissociation under equilibrium conditions,” The Journal of Chemical Physics 70, 2962–2978 (1979).
Dove and Raynor (1979) J. E. Dove and S. Raynor, “An ab initio calculation of the rate of vibrational relaxation and thermal dissociation of hydrogen by helium at high temperatures,” The Journal of Physical Chemistry 83, 127–133 (1979).
Dove et al. (1987) J. E. Dove, A. C. M. Rusk, P. H. Cribb, and P. G. Martin, “Excitation and dissociation of molecular hydrogen in shock waves at interstellar densities,” The Astrophysical Journal 318, 379 (1987).
Whitlock, Muckerman, and Roberts (1974) P. A. Whitlock, J. T. Muckerman, and R. E. Roberts, “Classical mechanics of recombination via the resonance complex mechanism: H + H + M → $\rm H_{2}$ + M for M = H, $\rm H_{2}$ , He, and Ar,” The Journal of Chemical Physics 60, 3658–3673 (1974).
McBride, Zehe, and Gordon (2002) B. J. McBride, M. J. Zehe, and S. Gordon, “NASA Glenn coefficicients for calculating thermodynamic properties of individual species,” Tech. Rep. (NASA NTRS Report/ Patent Number: NASA/TP-2002-211556, 2002).
Macdonald (2024) R. L. Macdonald, “State-to-state study of non-equilibrium recombination of oxygen and nitrogen molecules,” The Journal of Chemical Physics 160 (2024), 10.1063/5.0195238.
Singh and Schwartzentruber (2017) N. Singh and T. Schwartzentruber, “Nonequilibrium internal energy distributions during dissociation,” Proceedings of the National Academy of Sciences 115, 47–52 (2017).
Singh and Schwartzentruber (2020c) N. Singh and T. Schwartzentruber, “Non-boltzmann vibrational energy distributions and coupling to dissociation rate,” The Journal of Chemical Physics 152 (2020c), 10.1063/1.5142732.
Boyd (1977) R. K. Boyd, “The linear mixture formula in non-equilibrium chemical kinetics,” Canadian Journal of Chemistry 55, 802–811 (1977).
Dove, Halperin, and Raynor (1984) J. E. Dove, S. Halperin, and S. Raynor, “Deviations from the linear mixture rule in nonequilibrium chemical kinetics,” The Journal of Chemical Physics 81, 799–811 (1984).
Carruthers and Teitelbaum (1988) C. Carruthers and H. Teitelbaum, “The linear mixture rule in chemical kinetics. ii. thermal dissociation of diatomic molecules,” Chemical Physics 127, 351–362 (1988).
Mirahmadi and Pérez-Ríos (2022) M. Mirahmadi and J. Pérez-Ríos, “Three-body recombination in physical chemistry,” International Reviews in Physical Chemistry 41, 233–267 (2022).
Lindemann et al. (1922) F. A. Lindemann, S. Arrhenius, I. Langmuir, N. R. Dhar, J. Perrin, and W. C. McC. Lewis, “Discussion on “the radiation theory of chemical action”,” Trans. Faraday Soc. 17, 598–606 (1922).
Gilbert, Luther, and Troe (1983) R. G. Gilbert, K. Luther, and J. Troe, “Theory of thermal unimolecular reactions in the fall-off range. II. weak collision rate constants,” Berichte der Bunsengesellschaft für physikalische Chemie 87, 169–177 (1983).
Jasper and Miller (2014) A. W. Jasper and J. A. Miller, “Lennard–Jones parameters for combustion and chemical kinetics modeling from full-dimensional intermolecular potentials,” Combustion and Flame 161, 101–110 (2014).
Pack, Snow, and Smith (1972) R. T. Pack, R. L. Snow, and W. D. Smith, “On the mechanism of low-temperature termolecular atomic recombination,” The Journal of Chemical Physics 56, 926–932 (1972).
Whitlock, Muckerman, and Roberts (1972) P. Whitlock, J. Muckerman, and R. Roberts, “Classical dynamics of three-body recombination via the resonance complex mechanism,” Chemical Physics Letters 16, 460–463 (1972).
Orel (1987) A. E. Orel, “Nascent vibrational/ rotational distribution produced by hydrogen atom recombination,” The Journal of Chemical Physics 87, 314–318 (1987).
Esposito and Capitelli (2009) F. Esposito and M. Capitelli, “Selective vibrational pumping of molecular hydrogen via gas phase atomic recombination,” The Journal of Physical Chemistry A 113, 15307–15314 (2009).
Forrey (2013) R. C. Forrey, “Sturmian theory of three-body recombination: Application to the formation of $\rm H_{2}$ in primordial gas,” Physical Review A 88 (2013), 10.1103/physreva.88.052709.
Roberts, Bernstein, and Curtiss (1969) R. E. Roberts, R. B. Bernstein, and C. F. Curtiss, “Resonance theory of termolecular recombination kinetics: H+H+M $\rightarrow$ $\rm H_{2}$ +M,” The Journal of Chemical Physics 50, 5163–5176 (1969).
Leibowitz, Menard, and Stickford (1973) L. P. Leibowitz, W. A. Menard, and G. H. Stickford, “Radiative relaxation behind strong shock waves in hydrogen-helium mixtures,” in International Shock Tube Symposium (Stanford University Press, 1973) pp. 306–317.
Leibowitz and Kuo (1976) L. P. Leibowitz and T.-J. Kuo, “Ionizational nonequilibrium heating during outer planetary entries,” AIAA Journal 14, 1324–1329 (1976).

Treatment of Thermal Non-Equilibrium Dissociation Rates: Application to H2\text{H}_{2}

Abstract

I Introduction

II Theory of Non-Equilibrium Dissociation

II.1 State-Specific Rates and Master Equations

II.2 Thermal Rates

II.3 QSS Rates

II.3.1 Formulation and Solution for 𝑨𝟐​(𝑱,𝝂)\bm{A_{2}(J,\nu)}

II.3.2 Solution for 𝒌𝒅\bm{k_{d}}

II.3.3 Chemical Source Term

II.3.4 Transition Between 𝒌𝒅,n​r\bm{k}_{\bm{d},nr} and 𝒌𝒅,t​h\bm{k}_{\bm{d},th}

II.3.5 Solution for 𝒆𝒓​v\bm{e}_{\bm{r}v}

II.4 Pre-QSS Rates

II.4.1 Formulation

II.4.2 Approximate Solution for 𝑨2​(J,ν)\bm{A}_{2}(J,\nu) and 𝒌𝒅\bm{k}_{\bm{d}}

II.4.3 Chemical Source Term

II.4.4 Approximate Solution for 𝒆𝒓​v\bm{e}_{\bm{r}v}

III Application to the Dissociation of 𝐇𝟐\rm\bf H_{2}

III.1 Extracted Rate Constants

III.2 Fraction of Dissociation in the Pre-QSS and QSS Regions

III.3 Predicted Number Densities and Rovibrational Energies

IV Review of Rate Constants for 𝐇𝟐\rm\bf H_{2} Dissociation

IV.1 High Temperature/ Shock Tube Experiments

IV.2 Low Temperature/ Discharge-Flow Tube Experiments

IV.3 Computational Studies

IV.4 Comparison and Fits of Rate Constant Data

V Conclusion

Acknowledgements.

Author Declarations

Conflict of Interest

Author Contributions

Data Availability Statement

Appendix A Dissociation in Multi-Component Mixtures

Appendix B Recombination as a Second or Third-Order Process

Appendix C Rate Constants of Leibowitz

References

Treatment of Thermal Non-Equilibrium Dissociation Rates: Application to $\text{H}_{2}$

II.3.1 Formulation and Solution for $\bm{A_{2}(J,\nu)}$

II.3.2 Solution for $\bm{k_{d}}$

II.3.4 Transition Between $\bm{k}_{\bm{d},nr}$ and $\bm{k}_{\bm{d},th}$

II.3.5 Solution for $\bm{e}_{\bm{r}v}$

II.4.2 Approximate Solution for $\bm{A}_{2}(J,\nu)$ and $\bm{k}_{\bm{d}}$

II.4.4 Approximate Solution for $\bm{e}_{\bm{r}v}$

III Application to the Dissociation of $\rm\bf H_{2}$

IV Review of Rate Constants for $\rm\bf H_{2}$ Dissociation