2. The Master Equation

A full treatment of the energy states of each molecule is unfeasible for molecules larger than diatomics, as there are simply too many states. To simplify things we apply the RRKM approximation, which leaves the state of a molecule as a function of two quantities: the total energy \(E\) and total angular momentum quantum number \(J\). Frequently we will find that even this is too difficult, and will only keep the total energy \(E\) as an independent variable.

2.1. Isomers, Reactants, and Products

Throughout this document we will utilize the following terminology:

• An isomer is a unimolecular configuration on the potential energy surface.

• A reactant channel is a bimolecular configuration that associates to form an isomer. Dissociation from the isomer back to reactants is allowed.

• A product channel is a bimolecular configuration that is formed by dissociation of an isomer. Reassociation of products to the isomer is not allowed.

The isomers are the configurations for which we must model the energy states. We designate \(p_i(E, J, t)\) as the population of isomer \(i\) having total energy \(E\) and total angular momentum quantum number \(J\) at time \(t\). At long times, statistical mechanics requires that the population of each isomer approach a Boltzmann distribution \(b_i(E, J)\):

\[\lim_{t \rightarrow \infty} p_i(E, J, t) \propto b_i(E, J)\]

We can simplify by eliminating the angular momentum quantum number to get

\[p_i(E, t) = \sum_J p_i(E, J, t)\]

Let us also denote the (time-dependent) total population of isomer \(i\) by \(x_i(t)\):

\[x_i(t) \equiv \sum_J \int_0^\infty p_i(E, J, t) \ dE\]

The two molecules of a reactant or product channel are free to move apart from one another and interact independently with other molecules in the system. Accordingly, we treat these channels as fully thermalized, leaving as the only variable the total concentrations \(y_{n\ce{A}}(t)\) and \(y_{n\ce{B}}(t)\) of the molecules \(\ce{A}_n\) and \(\ce{B}_n\) of reactant channel \(n\). (Since the product channels act as infinite sinks, their populations do not need to be considered explicitly.)

Finally, we will use \(N_\mathrm{isom}\), \(N_\mathrm{reac}\), and \(N_\mathrm{prod}\) as the numbers of isomers, reactant channels, and product channels, respectively, in the system.

2.2. Collision Models

Bimolecular collisions with an inert species \(\ce{M}\) are the primary means by which an isomer molecule changes its energy. A reasonable estimate – although generally a bit of an underestimate – of the total rate of collisions \(k_{\mathrm{coll},i}(T)\) for each isomer \(i\) comes from Lennard-Jones collision theory:

\[k_{\mathrm{coll},i}(T) = \sqrt{\frac{8 k_\mathrm{B} T}{\pi \mu_i}} \pi d_i^2 \Omega_i^{(2,2)\ast}\]

Above, \(\mu_i\) is the reduced mass, \(d_i\) is the collision diameter, and \(k_\mathrm{B}\) is the Boltzmann constant. The collision diameter is generally taken as \(d \approx \frac{1}{2} (\sigma_i + \sigma_{\ce{M}})\), the arithmetic average of the Lennard-Jones \(\sigma\) parameter for the isomer and the bath gas. The parameter \(\Omega_i^{(2,2)\ast}\) represents a configurational integral, which is well-approximated by the expression

\[\Omega_i^{(2,2)\ast} = 1.16145 \tilde{T}^{-0.14874} + 0.52487 e^{-0.7732 \tilde{T}} + 2.16178 e^{-2.437887 \tilde{T}}\]

where \(\tilde{T} \equiv k_\mathrm{B} T / \sqrt{\epsilon_i \epsilon_{\ce{M}}}\) is a reduced temperature and \(\epsilon_i\) is the Lennard-Jones \(\epsilon\) parameter. Note that we have used a geometric average for the \(\epsilon\) parameters of the isomer and the bath gas in this expression. Assuming the total gas concentration to be constant and that the gas is ideal, we obtain an expression for the collision frequency \(\omega_i(T,P)\), which makes explicit the pressure dependence:

\[\omega_i(T,P) = k_{\mathrm{coll},i}(T) \frac{P}{k_\mathrm{B} T}\]

Now that we have an estimate for the total rate of collisions, we need to develop a model of the effect that these collisions have on the state of the isomer distribution. To this end, we define \(P(E, J, E^\prime, J^\prime)\) as the probability of a collision resulting in a transfer of a molecule from state \((E^\prime, J^\prime)\) to state \((E, J)\). There are two mathematical constraints on \(P(E, J, E^\prime, J^\prime)\). The first of these is normalization:

\[\sum_{J^\prime} (2 J^\prime + 1) \int_0^\infty P(E, J, E^\prime, J^\prime) \ dE^\prime = 1\]

The second of these is detailed balance, required in order to obtain the Boltzmann distribution at long times:

\[P(E, J, E^\prime, J^\prime) b(E^\prime, J^\prime) = P(E^\prime, J^\prime, E, J) b(E, J)\]

\[P_i(E^\prime, J^\prime, E, J) = \frac{\rho_i(E^\prime, J^\prime)}{\rho_i(E, J)} \exp \left(- \frac{E^\prime - E}{k_\mathrm{B} T} \right) P_i(E, J, E^\prime, J^\prime) \hspace{40pt} E < E^\prime\]

Rather than define models directly for \(P(E, J, E^\prime, J^\prime)\), we usually eliminate the angular momentum contribution and instead define \(P(E, E^\prime)\). This can be related to \(P(E, J, E^\prime, J^\prime)\) via

\[P(E, J, E^\prime, J^\prime) = P(E, E^\prime) \phi(E, J) = P(E, E^\prime) (2J + 1) \frac{\rho(E, J)}{\rho(E)}\]

where \(\rho(E) \equiv \sum_J (2J+1) \rho(E, J)\).

There are a variety of models used for \(P(E, E^\prime)\). By far the most common is the single exponential down model

\[P(E, E^\prime) = C(E^\prime) \exp \left( -\frac{E^\prime - E}{\alpha} \right) \hspace{40pt} E < E^\prime\]

where \(C(E^\prime)\) is determined from the normalization constraint. Note that this function has been defined for the deactivating direction (\(E < E^\prime\)) only, as the activating direction (\(E > E^\prime\)) is then set from detailed balance. The parameter \(\alpha\) corresponds to the average energy transferred in a deactivating collision \(\left< \Delta E_\mathrm{d} \right>\), which itself is a weak function of temperature.

Other models for \(P(E, J, E^\prime, J^\prime)\) include the Gaussian down

\[P(E, E^\prime) = C(E^\prime) \exp \left[ - \frac{(E^\prime - E)^2}{\alpha^2} \right] \hspace{40pt} E < E^\prime\]

and the double exponential down

\[P(E, E^\prime) = C(E^\prime) \left[ (1 - f) \exp \left( -\frac{E^\prime - E}{\alpha_1} \right) + f \exp \left( -\frac{E^\prime - E}{\alpha_2} \right) \right] \hspace{40pt} E < E^\prime\]

The parameters for these simple models generally contain so much uncertainty that more complex functional forms are generally not used.

2.3. Reaction Models

Chemical reaction events cause a change in molecular configuration at constant energy. The rate coefficient for this process must be determined as a function of energy rather than the usual temperature. Such a quantity is called a microcanonical rate coefficient and written as \(k(E, J)\). In the master equation we will differentiate between microcanonical rate coefficients for isomerization, dissociation, and association by using different letters: \(k_{ij}(E, J)\) for isomerization, \(g_{nj}(E, J)\) for dissociation, and \(f_{im}(E, J)\) for association. (By convention, we use indices \(i\) and \(j\) to refer to unimolecular isomers, \(m\) and \(n\) to refer to bimolecular reactant and product channels, and, later, \(r\) and \(s\) to refer to energy grains.)

As with collision models, the values of the microcanonical rate coefficients are constrained by detailed balance so that the proper equilibrium is obtained. The detailed balance expressions have the form

\[k_{ij}(E, J) \rho_j(E, J) = k_{ji}(E, J) \rho_i(E, J)\]

for isomerization and

\[f_{in}(E, J) \rho_n(E, J) = g_{ni}(E, J) \rho_i(E, J)\]

for association/dissociation, where \(\rho_i(E, J)\) is the density of states of the appropriate unimolecular or bimolecular configuration.

An alternative formulation incorporates the macroscopic equilibrium coefficient \(K_\mathrm{eq}(T)\) and equilibrium distributions \(b_i(E, J, T)\) at each temperature:

\[k_{ij}(E, J) b_j(E, J, T) = K_\mathrm{eq}(T) k_{ji}(E, J) b_i(E, J, T)\]

for isomerization and

\[f_{in}(E, J) b_n(E, J, T) = K_\mathrm{eq}(T) g_{ni}(E, J) b_i(E, J, T)\]

for association/dissociation. Note that these two formulations are equivalent if the molecular degrees of freedom are consistent with the macroscopic thermodynamic parameters. There are multiple reasons to use the latter formulation:

• Only the density of states of the unimolecular isomers need be computed. This is a result of the assumption of thermalized bimolecular channels, which means that we only need to compute the product \(f_{in} b_n\), and not the individual values of \(f_{in}\) and \(b_n\).

• Only the reactive rovibrational modes need be included in the density of states. Missing modes will not affect the observed equilibrium because we are imposing the macroscopic equilibrium via \(K_\mathrm{eq}(T)\).

• Constants of proportionality in the density of states become unimportant, as they cancel when taking the ratio \(\rho(E,J)/Q(\beta)\). For example, if the external rotational constants are unknown then we will include an active K-rotor in the density of states; this property means that the rotational constant of this active K-rotor cancels and is therefore arbitrary.

There are two common ways of determining values for \(k(E, J)\): the inverse Laplace transform method and RRKM theory. The latter requires detailed information about the transition state, while the former only requires the high-pressure limit rate coefficient \(k_\infty(T)\).

2.3.1. Inverse Laplace Transform

The microcanonical rate coefficient \(k(E)\) is related to the canonical high-pressure limit rate coefficient \(k_\infty(T)\) via a Boltzmann averaging

\[k_\infty(T) = \frac{\sum_J \int_0^\infty k(E) \rho(E, J) e^{-\beta E} \ dE}{\sum_J \int_0^\infty \rho(E, J) e^{-\beta E} \ dE}\]

where \(\rho(E, J)\) is the rovibrational density of states for the reactants and \(\beta \equiv (k_\mathrm{B} T)^{-1}\). Neglecting the angular momentum dependence, the above can be written in terms of Laplace transforms as

\[k_\infty(T) = \frac{\mathcal{L} \left[ k(E) \rho(E) \right]} {\mathcal{L} \left[ \rho(E) \right]} = \frac{\mathcal{L} \left[ k(E) \rho(E) \right]} {Q(\beta)}\]

where \(Q(\beta)\) is the rovibrational partition function for the reactants. The above implies that \(E\) and \(\beta\) are the transform variables. We can take an inverse Laplace transform in order to solve for \(k(E)\):

\[k(E) = \frac{\mathcal{L}^{-1} \left[ k_\infty(\beta) Q(\beta) \right] }{\rho(E)}\]

Hidden in the above manipulation is the assumption that \(k_\infty(\beta)\) is valid over a temperature range from zero to positive infinity.

The most common form of \(k_\infty(T)\) is the modified Arrhenius expression

\[k(T) = A T^n \exp \left( - \frac{E_\mathrm{a}}{k_\mathrm{B} T} \right)\]

where \(A\), \(n\), and \(E_\mathrm{a}\) are the Arrhenius preexpoential, temperature exponent, and activation energy, respectively. For \(n = 0\) and \(E_\mathrm{a} > 0\) the inverse Laplace transform can be easily evaluated to give

\[k(E) = A \frac{\rho(E - E_\mathrm{a})}{\rho(E)} \hspace{40pt} E > E_\mathrm{a}\]

We can also determine an expression when \(n > 0\) and \(E_\mathrm{a} > 0\) using a convolution integral:

\[k(E) = A \frac{\phi(E - E_\mathrm{a})}{\rho(E)} \hspace{40pt} E > E_\mathrm{a}\]

\[\phi(E) = \mathcal{L}^{-1} \left[ T^n Q(\beta) \right] = \frac{1}{k_\mathrm{B}^n \Gamma(n)} \int_0^E (E - x)^{n-1} \rho(x) \ dx\]

Finally, for cases where \(n < 0\) and/or \(E_\mathrm{a} < 0\) we obtain a rough estimate by lumping these contributions into the preexponential at the temperature we are working at. By redoing this at each temperature being considered we minimize the error introduced, at the expense of not being able to identify a single \(k(E)\).

2.3.2. RRKM Theory

RRKM theory – named for Rice, Ramsperger, Kassel, and Marcus – is a microcanonical transition state theory. Like canonical transition state theory, detailed information about the transition state and reactants are required, e.g. from a quantum chemistry calculation. If such information is available, then the microcanonical rate coefficient can be evaluated via the equation

\[k(E, J) = \frac{N^\ddagger(E, J)}{h \rho(E, J)}\]

where \(N^\ddagger(E, J)\) is the sum of states of the transition state, \(\rho(E, J)\) is the density of states of the reactant, and \(h\) is the Planck constant. Both the transition state and the reactants have been referenced to the same zero of energy. The sum of states is related to the density of states via

\[N(E, J) = \int_0^E \rho(x, J) \ dx\]

The angular momentum quantum number dependence can be removed via

\[k(E) = \sum_J (2J+1) k(E, J)\]

2.4. The Full Master Equation

The governing equation for the population distributions \(p_i(E, J, t)\) of each isomer \(i\) and the reactant concentrations \(y_{n\mathrm{A}}(t)\) and \(y_{n\mathrm{B}}(t)\) combines the collision and reaction models to give a linear integro-differential equation:

\[ \begin{align}\begin{aligned}\begin{split}\frac{d}{dt} p_i(E, J, t) &= \omega_i(T, P) \sum_{J^\prime} \int_0^\infty P_i(E, J, E^\prime, J^\prime) p_i(E^\prime, J^\prime, t) \ dE^\prime - \omega_i(T, P) p_i(E, J, t) \\ & \mbox{} + \sum_{j \ne i}^{N_\mathrm{isom}} k_{ij}(E, J) p_j(E, J, t) - \sum_{j \ne i}^{N_\mathrm{isom}} k_{ji}(E, J) p_i(E, J, t) \\ & \mbox{} + \sum_{n=1}^{N_\mathrm{reac}} y_{n\mathrm{A}}(t) y_{n\mathrm{B}}(t) f_{in}(E, J) b_n(E, J, t) - \sum_{n=1}^{N_\mathrm{reac} + N_\mathrm{prod}} g_{ni}(E, J) p_i(E, J, t) \\\end{split}\\\begin{split}\frac{d}{dt} y_{n\mathrm{A}}(t) = \frac{d}{dt} y_{n\mathrm{B}}(t) &= \sum_{i=1}^{N_\mathrm{isom}} \int_0^\infty g_{ni}(E, J) p_i(E, J, t) \ dE \\ & \mbox{} - \sum_{i=1}^{N_\mathrm{isom}} y_{n\mathrm{A}}(t) y_{n\mathrm{B}}(t) \int_0^\infty f_{in}(E, J) b_n(E, J, t) \ dE\end{split}\end{aligned}\end{align} \]

A summary of the variables is given below:

Variable	Meaning
\(p_i(E, J, t)\)	Population distribution of isomer \(i\)
\(y_{n\mathrm{A}}(t)\)	Total population of species \(\ce{A}_n\) in reactant channel \(n\)
\(\omega_i(T, P)\)	Collision frequency of isomer \(i\)
\(P_i(E, J, E^\prime, J^\prime)\)	Collisional transfer probability from \((E^\prime, J^\prime)\) to \((E, J)\) for isomer \(i\)
\(k_{ij}(E, J)\)	Microcanonical rate coefficient for isomerization from isomer \(j\) to isomer \(i\)
\(f_{im}(E, J)\)	Microcanonical rate coefficient for association from reactant channel \(m\) to isomer \(i\)
\(g_{nj}(E, J)\)	Microcanonical rate coefficient for dissociation from isomer \(j\) to reactant or product channel \(n\)
\(b_n(E, J, t)\)	Boltzmann distribution for reactant channel \(n\)
\(N_\mathrm{isom}\)	Total number of isomers
\(N_\mathrm{reac}\)	Total number of reactant channels
\(N_\mathrm{prod}\)	Total number of product channels

The above is called the two-dimensional master equation because it contains two dimensions: total energy \(E\) and total angular momentum quantum number \(J\). In the first equation (for isomers), the first pair of terms correspond to collision, the second pair to isomerization, and the final pair to association/dissociation. Similarly, in the second equation above (for reactant channels), the pair of terms refer to dissociation/association.

We can also simplify the above to the one-dimensional form, which only has \(E\) as a dimension:

\[ \begin{align}\begin{aligned}\begin{split}\frac{d}{dt} p_i(E, t) &= \omega_i(T, P) \int_0^\infty P_i(E, E^\prime) p_i(E^\prime, t) \ dE^\prime - \omega_i(T, P) p_i(E, t) \\ & \mbox{} + \sum_{j \ne i}^{N_\mathrm{isom}} k_{ij}(E) p_j(E, t) - \sum_{j \ne i}^{N_\mathrm{isom}} k_{ji}(E) p_i(E, t) \\ & \mbox{} + \sum_{n=1}^{N_\mathrm{reac}} y_{n\mathrm{A}}(t) y_{n\mathrm{B}}(t) f_{in}(E) b_n(E, t) - \sum_{n=1}^{N_\mathrm{reac} + N_\mathrm{prod}} g_{ni}(E) p_i(E, t) \\\end{split}\\\begin{split}\frac{d}{dt} y_{n\mathrm{A}}(t) = \frac{d}{dt} y_{n\mathrm{B}}(t) &= \sum_{i=1}^{N_\mathrm{isom}} \int_0^\infty g_{ni}(E) p_i(E, t) \ dE \\ & \mbox{} - \sum_{i=1}^{N_\mathrm{isom}} y_{n\mathrm{A}}(t) y_{n\mathrm{B}}(t) \int_0^\infty f_{in}(E) b_n(E, t) \ dE\end{split}\end{aligned}\end{align} \]

The equations as given are nonlinear, both due to the presence of the bimolecular reactants and because both \(\omega_i\) and \(P_i(E, E^\prime)\) depend on the composition, which is changing with time. The rate coefficients can be derived from considering the pseudo-first-order situation where \(y_{n\mathrm{A}}(t) \ll y_{n\mathrm{B}}(t)\), and all \(y(t)\) are negligible compared to the bath gas \(\ce{M}\). From these assumptions the changes in \(\omega_i\), \(P_i(E, E^\prime)\), and all \(y_{n\mathrm{B}}\) can be neglected, which yields a linear equation system.

2.5. The Energy-Grained Master Equation

Except for the simplest of unimolecular reaction networks, both the one-dimensional and two-dimensional master equation must be solved numerically. To do this we must discretize and truncate the energy domain into a finite number of discrete bins called grains. This converts the linear integro-differential equation into a system of first-order ordinary differential equations:

\[\begin{split}\renewcommand{\vector}[1]{\boldsymbol{\mathbf{#1}}} \renewcommand{\matrix}[1]{\boldsymbol{\mathbf{#1}}} \frac{d}{dt} \begin{bmatrix} \vector{p}_1 \\ \vector{p}_2 \\ \vdots \\ y_{1\mathrm{A}} \\ y_{2\mathrm{A}} \\ \vdots \end{bmatrix} = \begin{bmatrix} \matrix{M}_1 & \matrix{K}_{12} & \ldots & \matrix{F}_{11} \vector{b}_1 y_{1\mathrm{B}} & \matrix{F}_{12} \vector{b}_2 y_{2\mathrm{B}} & \ldots \\ \matrix{K}_{21} & \matrix{M}_2 & \ldots & \matrix{F}_{21} \vector{b}_1 y_{1\mathrm{B}} & \matrix{F}_{22} \vector{b}_2 y_{2\mathrm{B}} & \ldots \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \\ (\vector{g}_{11})^T & (\vector{g}_{12})^T & \ldots & h_1 & 0 & \ldots \\ (\vector{g}_{21})^T & (\vector{g}_{22})^T & \ldots & 0 & h_2 & \ldots \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \end{bmatrix} \begin{bmatrix} \vector{p}_1 \\ \vector{p}_2 \\ \vdots \\ y_{1\mathrm{A}} \\ y_{2\mathrm{A}} \\ \vdots \end{bmatrix}\end{split}\]

The diagonal matrices \(\matrix{K}_{ij}\) and \(\matrix{F}_{in}\) and the vector \(\vector{g}_{ni}\) contain the microcanonical rate coefficients for isomerization, association, and dissociation, respectively:

\[\begin{split}(\matrix{K}_{ij})_{rs} &= \begin{cases} \frac{1}{\Delta E_r} \int_{E_r - \Delta E_r/2}^{E_r + \Delta E_r/2} k_{ij}(E) \, dE & r = s \\ 0 & r \ne s \end{cases} \\ (\matrix{F}_{in})_{rs} &= \begin{cases} \frac{1}{\Delta E_r} \int_{E_r - \Delta E_r/2}^{E_r + \Delta E_r/2} f_{in}(E) \, dE & r = s \\ 0 & r \ne s \end{cases} \\ (\matrix{g}_{ni})_r &= \frac{1}{\Delta E_r} \int_{E_r - \Delta E_r/2}^{E_r + \Delta E_r/2} g_{ni}(E) \, dE\end{split}\]

The matrices \(\matrix{M}_i\) represent the collisional transfer probabilities minus the rates of reactive loss to other isomers and to reactants and products:

\[\begin{split}(\matrix{M}_i)_{rs} = \begin{cases} \omega_i \left[ P_i(E_r, E_r) - 1 \right] - \sum_{j \ne i}^{N_\mathrm{isom}} k_{ij}(E_r) - \sum_{n=1}^{N_\mathrm{reac} + N_\mathrm{prod}} g_{ni}(E_r) & r = s \\ \omega_i P_i(E_r, E_s) & r \ne s \end{cases}\end{split}\]

The scalars \(h_n\) are simply the total rate coefficient for loss of reactant channel \(n\) due to chemical reactions:

\[h_n = - \sum_{i=1}^{N_\mathrm{isom}} \sum_{r=1}^{N_\mathrm{grains}} y_{n\mathrm{B}} f_{in}(E_r) b_n(E_r)\]

2.6. Further Reading

The interested reader is referred to any of a variety of other sources for alternative presentations, of which an illustrative sampling is given here [Gilbert1990] [Baer1996] [Holbrook1996] [Forst2003] [Pilling2003].

[Gilbert1990]

R. G. Gilbert and S. C. Smith. Theory of Unimolecular and Recombination Reactions. Blackwell Sci. (1990).

[Baer1996]

T. Baer and W. L. Hase. Unimolecular Reaction Dynamics. Oxford University Press (1996).

[Holbrook1996]

K. A. Holbrook, M. J. Pilling, and S. H. Robertson. Unimolecular Reactions. Second Edition. John Wiley and Sons (1996).

[Forst2003]

W. Forst. Unimolecular Reactions: A Concise Introduction. Cambridge University Press (2003).

[Pilling2003]

M. J. Pilling and S. H. Robertson. Annu. Rev. Phys. Chem. 54, p. 245-275 (2003). doi:10.1146/annurev.physchem.54.011002.103822