While hydrogen is a promising source of clean energy, the safety and optimization of hydrogen technologies rely on controlling ignition through explosion limits: pressure-temperature boundaries separating explosive behavior from comparatively slow burning. Here, we show that the emergent nonequilibrium chemistry of combustible mixtures can exhibit the quantitative features of a phase transition. With stochastic simulations of the chemical kinetics for a model mechanism of hydrogen combustion, we show that the boundaries marking explosive domains of kinetic behavior are nonequilibrium critical points. Near the pressure of the second explosion limit, these critical points terminate the transient coexistence of dynamical phases—one that autoignites and another that progresses slowly. Below the critical point temperature, the chemistry of these phases is indistinguishable. In the large system limit, the pseudo-critical temperature converges to the temperature of the second explosion limit derived from mass-action kinetics.

Heating a gaseous mixture of hydrogen and oxygen can ignite an explosion.1 Whether a mixture ignites spontaneously not only depends on the temperature but also on the pressure and composition. Together, these variables control the thermal and kinetic feedback mechanisms that can switch on and off across “explosion limits”—sharp boundaries that marble the thermodynamic state space. Knowing the precise location of explosion limits is critical to technological applications, including the safe use of hydrogen as an alternative fuel in combustion engines, in fuel-cell based electric vehicles, and as rocket propellant in space exploration.2–5 At the macroscopic level, these limits are understood, in part, through empirical models of the kinetics.6 However, these models assume mean-field rate laws that neglect even slight fluctuations in energy or radical numbers that can stimulate ignition. While there have been previous computational and theoretical studies,7–9 it remains unclear how these fluctuations affect explosion limits and how ignition phenomena emerge from finite-size systems. Understanding the role of fluctuations is also relevant to the prediction of ignition times from the relatively small systems accessible with molecular simulations10 and reactive force fields.11–15 Fluctuations in radical numbers and energy are also particularly important for the use of hydrogen as fuel because of its low ignition energy and susceptibility to autoignition, which impact both hydrogen safety and engine performance.

Explosion limits mark the macroscopic domains of pressure and temperature where a gaseous mixture can ignite. These boundaries can have an intricate structure, and for hydrogen oxidation, the limits have a well-known Z-shape.16 The branches of the “Z” are the first, second, and third explosion limits, in order of increasing pressure. Mixtures at pressures and temperatures to the right of these boundaries can spontaneously ignite. The sudden changes which mixtures of hydrogen and oxygen exhibit upon a small change in thermodynamic conditions around these dividing lines are reminiscent of phase transitions.17,18 Nonequilibrium phase transitions19,20 are accumulating interest recently in statistical physics, due in part to their connection to the glass transition21 and protein folding.22 In this article, we add another dimension to this work, which is largely focused on condensed phase problems under steady-state conditions, by demonstrating that the stochastic evolution of combusting mixtures in the gas phase is a transient nonequilibrium phase transition. We build on recent work advancing the statistical foundations of thermodynamics away from equilibrium that largely focus on relationships governing fluctuations in observables23 and their large deviations,24 which can determine the fate of both chemical and physical systems.

At the core of many efforts in nonequilibrium statistical mechanics is the idea of replacing the notion of “states” with that of “trajectories.”25,26 Here, we modify this perspective for a model of the gas-phase reactions near the second explosion limit of hydrogen, which is the reactive regime where the chemical heat release is small and a kinetic feedback mechanism stimulates ignition, not thermal runaway.1,16 This regime contrasts the second-extended limit, which is an extension of the second limit to higher pressures and temperatures and divides “strong” and “weak” ignition events. For hydrogen combustion at fixed temperature and pressure, we show that ignition is a phase transition in the space of trajectories, where each trajectory evolves through the space of mixture compositions. Though not analyzed from this perspective, similar phase transitions are known in prototypical models of autocatalytic chemical reactions under steady-state conditions.8,27,28 The present work, however, uses a microscopic model based on the chemical master equation to account for fluctuations in molecule (radical) numbers and experimental data to explicitly reproduce the ignition of hydrogen-oxygen mixtures. Through a finite-size scaling analysis of the fluctuations in trajectory observables, we show that the stochastic trajectories executing chemical reactions exhibit the signatures of both phase coexistence and criticality.

Hydrogen oxidation is an example of “complex chemistry.”29 The “complexity” is that there are roughly 8 chemical species, and while there is no universal mechanism, on the order of 20 elementary reactions involved in the overall reaction, 2H2 + O2 2H2O. Elementary reactions involve transient, free radicals that cancel out in the overall stoichiometry. Here, we study a population of N molecules initially divided stoichiometrically between hydrogen and oxygen. Over time, the mixture transforms through a set of elementary reactions into H2O and HO2. For pressures and temperatures near the second explosion limit, there is evidence that a reduced set of five elementary reactions is sufficient to estimate ignition-delay times and explosion limits.16,30,31 The mechanism we use consists of these five irreversible reactions

H2+O2k12OH,
(R1)
OH+H2k2H+H2O,
(R2)
H+O2k3OH+O,
(R3)
O+H2k4OH+H,
(R4)
H+O2+Mk5HO2+M,
(R5)

where M is an arbitrary collision partner. For each reaction, Ri, we model the temperature dependence of the rate constants, ki, with the modified Arrhenius equation, as is common for chain-branching reactions.

Quantum mechanical calculations coupled with transition state theory suggest that the initiation reaction (R1) is unlikely to proceed directly32 but could be a composite of two reactions included in some full models near the second explosion limit: the initiation reactions H2 + O2 H + HO2 and H + HO2 2OH.6 Termolecular reactions have been proposed as an alternative means of initiation (e.g., 2O2 + H2 2HO2 and 2HO2 2OH + O232) and are expected to be particularly important at high pressures.33 We will show data for the mechanism above, but we arrive at the same conclusions with other H2/O2 mechanisms (see the supplementary material). As with any mechanism, the reactions we use capture some of the basic features of chemical explosions but we do not expect this model to be applicable outside the temperatures and pressures of the second limit.

Massively parameterized chemical kinetic models form the foundation of our current understanding of chain-branching reactions, but, even with extensive experimental validation, these models may still be inaccurate, incomplete, or only work for a specific range of temperatures and pressures.7,16,36–50 Near the second explosion limit of hydrogen, chain branching reactions are largely responsible for stimulating ignition.1,7,16,30 Thus, we fix the temperature so that the energy released by these reactions as heat is effectively dissipated.1,16 The mechanism we adopt models the chemical reactions only at moderate temperatures and pressures. Consequently, the mechanism does not reproduce the first, second extended, or third explosion limits.

Deterministic rate equations are a mean-field approximation of the chemical master equation which assumes that the molecule numbers are self-averaging observables and decay in the infinite system-size scaling limit. To simulate these reactions and explicitly account for fluctuations due to the formation and destruction of molecules through the reactions above (Fig. S13 of the supplementary material), we use kinetic Monte Carlo and the Doob-Gillespie direct stochastic simulation algorithm.34,35 This algorithm generates trajectories consistent with the numerical solution of the chemical master equation. Assuming the system is well mixed, the time, τ, between reactions, Ri, is an exponentially distributed random variable with a probability p(τ,i)=aiea0τdτ. This probability depends on the mixture composition of (or number of molecules of each) chemical species, X(t)={X1,X2,,XM}, through the propensity of a single reaction, ai, and the total propensity, a0=iai. We transform the macroscopic rate constants, all of which are from experimental measurements, into stochastic rate constants through standard formulas for bimolecular and termolecular reactions.7,34,35 In the propensities, we include the collision efficiency for the termolecular reaction (R5) from Ref. 6 with M: 18.6 for H2O, 2.86 for H2, and 1.00 for both O2 and HO2.

Each simulated trajectory evolves a stoichiometric mixture of hydrogen and oxygen to a mixture of water and hydroperoxyl radical. Since the initial composition is a stoichiometric mixture of H2 and O2, and the final composition is a mixture of H2O and HO2, the mixture evolves away from chemical equilibrium during the course of the reaction. All reactions are treated as irreversible and so do not satisfy detailed balance. We do impose thermal equilibrium (isothermal reactor) and diffusive equilibrium since the system is “well mixed.” The composition of the mixture after the reaction is complete depends on the temperature, T, and pressure, P. Conditions of fixed NPT precisely define the explosion limit, which is a family of (T, P) points. Our principal conclusions, however, are not altered when working at NVT (Fig. S15 of the supplementary material), though these ensembles are not guaranteed to be equivalent for the finite-size systems we consider here. To fix the pressure, we update the volume after each reaction event through the ideal gas equation of state, instead of having an inflow of reactants or an outflow of products, for example.

Under fixed temperature and pressure conditions, the reduced mechanism for hydrogen oxidation does have a well-defined kinetic condition that marks the second explosion limit. This condition, 2k3 = k5[M], is based on the rate constants of the reactions that annihilate and create radical species and, so, only holds in the infinite-size scaling limit when the macroscopic, mass-action rate equations are valid. The temperature and pressure of the second explosion limit, 2k3 = k5P/RT, are fixed by the rate of the branching reaction, (R3), and the rate of the termination reaction, (R5),30 where R is the gas constant. Using the reduced mechanism, and the rate constants from Ref. 7, the explosion limit temperature is 920 K at 1 atm. For our stochastic simulations, we consider initially stoichiometric mixtures of finite numbers of molecules N at P = 1 atm, but, in the N limit, mixtures can spontaneously ignite when heated to temperatures above 920 K.

Over the course of each stochastic trajectory, an initially stoichiometric mixture of hydrogen and oxygen transforms into the product water. On short time scales during the reaction, the intermediate concentrations might be treated using mass-action kinetics as being at quasi-steady-state. Here, we simulate the fluctuations in the number of each radical species. At the end of the overall reaction, some byproducts can also remain, mostly HO2, since the reaction mechanism we use treats it as unreactive on the time scales relevant to ignition.6 To characterize these dynamical histories through the space of compositions, we choose an order parameter measuring the dynamical activity, K. We define the K of a trajectory as the number of reactions occurring in a mixture of N molecules at a pressure P and temperature T over an observation time 0ttobs, K[x(tobs)]=i1(Δti). The indicator function 1(Δti)=1 if a reaction event occurs and is zero otherwise. We divide tobs uniformly into intervals, Δt, apply the indicator function in each, and sum over time intervals to get the activity.

For systems in nonequilibrium steady-states, the activity is linearly extensive in both time and system size.21,22 However, in the transient process of hydrogen combustion, the activity is not extensive in time; the activity does grow during the course of the reaction, but it reaches a stationary value once the reaction is complete. We have found this behavior to be typical of several hydrogen combustion mechanisms, and in this particular reduced model, it is a result of neither H2O nor HO2 being able to transform into other chemical species. The composition of these final products varies with temperature [Fig. S5(b) of the supplementary material]. As a consequence, the time is not a property whose infinite scaling limit can reveal a phase transition, and, in fact, the limit tobs would prevent an analysis of transient ignition processes that occur after finite ignition-delay times. Consequently, we treat time as an intensive control variable, which is an approach that is different from, but complements, the traditional application of the thermodynamics of trajectories.51 

During the reaction, the temporal profile of the activity depends qualitatively on whether the conditions imposed on a mixture place it above or below the explosion limit (Fig. S12 of the supplementary material). Below the explosion limit temperature, the growth in the activity is steady while the reaction progresses. But, above the explosion limit temperature, the activity undergoes a rapid jump that signals an ignition event, after a time delay that depends on the macroscopic conditions and the stochastic trajectory. In our simulations of the reduced model, the activity reaches a stationary value after the reaction is complete. The stationary activity value is extensive in the initial number of molecules, N [Fig. S7(b) of the supplementary material shows the linear scaling]. The rate of linear growth is dependent upon the temperature, with more rapid divergence at higher temperatures. These dependencies motivate our use of inverse temperature and the time as the conjugate fields to the activity density, K/N.

For all system sizes, we define the zero of time to be when the mixture is stoichiometric and no reactions have taken place, tobs = ti. We define the final time as the time at which the overall reaction is complete, tobs = tf. As a consequence of these temporal boundary conditions, K(ti) = 0 and the activity density is K(tf)N1 (for this model). For the 106 trajectories at each T and P, we can shift the zero of time to that of the first (initiation) reaction. We have found that this definition gives good agreement between the ignition time from extensive molecular dynamics simulations and the predictions of GRI Mech 3.0 from CHEMKIN.10 With this definition, the clock of every simulation starts at the first reaction event and the subsequent variation in the ignition delay is a result of the ensuing chemistry. However, choosing this time zero does not impact our conclusions below.

Near the combination of thermodynamic fields that define the second explosion limit, we observe the signatures of a dynamical phase transition (Fig. 1). These signatures appear in the integrated activity of 106 simulated trajectories over temperatures ranging from 800 K to 1000 K. From the trajectory ensemble, we focus on the time evolution of the activity distribution. When the system has the potential to ignite at a fixed T and P, there is a sudden change in this distribution that signals a dynamical phase transition: At sufficiently high temperatures, for example, the activity distribution is bimodal at a time t* [Fig. 1(c)], when the maxima of the peaks in the activity distribution have equal heights. This time depends on the temperature, but at much earlier times or later times the distribution is unimodal (movies in Sec. S4 of the supplementary material).

FIG. 1.

(a) At sufficiently high temperatures, the probability distributions of dynamical activity per molecule, K/N, are bimodal over a range of times during the combustion of hydrogen. Bimodality reflects the separation of trajectories into active and inactive dynamical phases. (b) Dynamical phase diagram using temperature and time as control variables and the most probable activity as the order parameter. A solid black line marks where active and inactive trajectories coexist in equal proportion. Dashed lines mark the coexistence region of the finite system size, when the inactive (active) phase is 95% of all sampled trajectories. Coexistence ends at a pseudo-critical point, roughly when the two peaks are unresolvable. (c) Probability distributions at four times and temperatures when dynamical phases coexist in equal proportion (e.g., at 900 K, t*=3 s). As temperature decreases from 900 K to 840 K, the peaks merge, indicating the approach to a pseudo-critical point with temperature, Tc(N). These representative data are for a system with N = 999 molecules, P = 1 atm, and Tc(N=999)840 K.

FIG. 1.

(a) At sufficiently high temperatures, the probability distributions of dynamical activity per molecule, K/N, are bimodal over a range of times during the combustion of hydrogen. Bimodality reflects the separation of trajectories into active and inactive dynamical phases. (b) Dynamical phase diagram using temperature and time as control variables and the most probable activity as the order parameter. A solid black line marks where active and inactive trajectories coexist in equal proportion. Dashed lines mark the coexistence region of the finite system size, when the inactive (active) phase is 95% of all sampled trajectories. Coexistence ends at a pseudo-critical point, roughly when the two peaks are unresolvable. (c) Probability distributions at four times and temperatures when dynamical phases coexist in equal proportion (e.g., at 900 K, t*=3 s). As temperature decreases from 900 K to 840 K, the peaks merge, indicating the approach to a pseudo-critical point with temperature, Tc(N). These representative data are for a system with N = 999 molecules, P = 1 atm, and Tc(N=999)840 K.

Close modal

At t*, each set of trajectories populating a distinct band of activity coincides with a dynamical phase. Within the sampled trajectories, there is an autoigniting “active” phase associated with the maximal peak height at Ka=Nκa and a slow burning “inactive” phase associated with the lower activity peak located at Ki=Nκi. Two peaks in the activity distribution reflect the transient bistability of trajectory space and the transient coexistence of these dynamical phases. While we show data for a single reduced mechanism, our main results appear to be robust to the choice of H2/O2 mechanism (see the supplementary material). These data suggest that the coexistence of dynamical phases is an alternative perspective of the well-known competition between branching and termination of chain carriers—the key elements of combustion reaction mechanisms.

Each dynamical phase has a distinct composition of chemical species and reactions. For example, mixtures at t* that are undergoing high activity trajectories are largely composed of products, HO2 and H2O. Those undergoing low activity trajectories are largely composed of reactants, H2 and O2. As one would expect, the mole fractions of each species in the active and inactive trajectories (Fig. S14 of the supplementary material) reflect the frequency of reactions in the reduced mechanism. Analyzing the composition of reactions in each phase, the low activity (inactive) phase trajectories are dominated by termination reactions (R5) that annihilate radicals, while the high activity (active) phase consists of trajectories dominated by branching reactions (R3) and (R4) that create radicals (Sec. S5 of the supplementary material). These data suggest that the branching and termination reactions are the most critical elements of the mechanism for bimodal activity distributions. It is important to emphasize that both phases at the time t* are undergoing combustion, but the active phase consists of mixtures that have an initial proliferation of radicals that precipitate ignition.

The time t* of dynamical phase coexistence, marked by the bimodal distribution function for the order parameter, K, depends strongly on temperature. Using temperature and time as control variables, we map the most probable activity per molecule, Kmp/N. The resulting phase diagram in Fig. 1(b) shows a jump from κi to κa in the most probable activity per molecule at t* (solid black line). Because the long-time limit of the activity is linearly extensive in the number of molecules, the density Kmp/N has an approximate range of [0,1] (Fig. S7 of the supplementary material). Above the phase boundary, all trajectories are in the active phase, and below the boundary, all trajectories are in the inactive phase. An increase in time or temperature that crosses the coexistence boundary causes the most abundant dynamical phase to switch from inactive to active. Also marked on the phase diagram is the coexistence region, which has a finite width because of the finite system size. We find that the width of this region, however, decays as 1N at temperatures above and below the threshold for coexistence (Fig. S6 of the supplementary material). This finding is supported by the fact that each trajectory is consistent with the chemical master equation, which has fluctuations in the mixture composition that decay similarly and yield mean-field rate equations in the infinite size limit. Proofs of this fact, however, are only possible for simpler reaction schemes.34,35

At temperatures and times along the phase boundary, the active and inactive peaks are of equal height. A similar diagram results, however, using the “equal area rule.”52 Moving along the coexistence line in the direction of increasing time, the distance between the peak heights, ΔK=KaKi at t* measures the activity change that accompanies the sudden burst of reactions upon ignition. This quantity measures the “latent activity” associated with the transformation of a gaseous mixture over a small time window. In a first order phase transition, there is generally a radical change in structure. Here, the radical change is in the composition of the mixture.

First order transitions often (but not always) end up at a critical point where there is a second order phase transition. In Fig. 1(c), the distance between peaks’ heights shrinks upon lowering the temperature, which suggests that the coexistence terminates at a pseudo-critical point. A rough estimate of the pseudo-critical temperature, Tc(N), from the activity distributions is the temperature below which only one peak is observable for all simulation times. The chemistry of the trajectories in the branching rich active phase and the termination rich inactive phase is distinguishable at temperatures above the critical temperature. But, they become increasingly similar until the critical point where they ultimately merge. We will make a more precise estimate of the critical temperature, however, through finite-size scaling, which has had success both in statistical mechanics53,54 and quantum mechanics.55,56

Coexistence between two dynamical phases suggests a first-order phase transition in the trajectory space of igniting mixtures. To quantitatively establish that the dynamical phase change is a first order coexistence line terminating at a critical point, we analyze the magnitude of the activity jump, ΔK, as a function of temperature for a range of system sizes from N = 999 to 666, 666 (Fig. 2). These system sizes are necessary for the initial mixture to have stoichiometric amounts of H2 and O2. The crossover in this order parameter as a function of temperature sharpens with increasing system size, suggesting a finite discontinuity in the infinite size limit [Fig. 2(a)]. This finding is additional evidence that the trajectory observables reflect a phase change in the combustion of finite-size systems and that the onset of ignition becomes sharp with a well-defined temperature and time t* in the infinite-size limit.

FIG. 2.

(a) The magnitude of the activity jump, ΔKN=κaκi, as a function of temperature for system sizes ranging from 999 to 666, 666 molecules initially. The crossover in this order parameter as a function of temperature sharpens with increasing system size. (b) Estimates of the pseudo-critical temperature, Tc(N) (in Kelvin), using the temperature at which ΔK/N exceeds 0.02 (chosen to be larger than our statistical errors in Ka and Ki). The pseudo-critical temperature converges to 920 K in the large N limit, the mass-action rate law prediction of the explosion temperature. Inset: the P-T second explosion limit, T = k5P/2k3 = 920 K, at 1 atm for an initially stoichiometric mixture. (c) The pseudo-critical time, tc(N), as a function of system size. Inset: magnification showing the convergence of tc(N) to 0.78 seconds at N = 666, 666.

FIG. 2.

(a) The magnitude of the activity jump, ΔKN=κaκi, as a function of temperature for system sizes ranging from 999 to 666, 666 molecules initially. The crossover in this order parameter as a function of temperature sharpens with increasing system size. (b) Estimates of the pseudo-critical temperature, Tc(N) (in Kelvin), using the temperature at which ΔK/N exceeds 0.02 (chosen to be larger than our statistical errors in Ka and Ki). The pseudo-critical temperature converges to 920 K in the large N limit, the mass-action rate law prediction of the explosion temperature. Inset: the P-T second explosion limit, T = k5P/2k3 = 920 K, at 1 atm for an initially stoichiometric mixture. (c) The pseudo-critical time, tc(N), as a function of system size. Inset: magnification showing the convergence of tc(N) to 0.78 seconds at N = 666, 666.

Close modal

To determine the critical temperature, Tc(N)=Tc, we again use the finite-size scaling of the “latent activity” of ignition, ΔK, at the time t*. For a given system size, we identify Tc(N) with the temperature at which ΔK/N exceeds 0.02. At this value, active and inactive phases are resolvable in the activity distributions within our statistical errors (Sec. S7 of the supplementary material). Despite a dependence on the total number of molecules, N, this pseudo-critical temperature is still an estimate of the lowest temperature at which there is a range of times with a bimodal activity distribution. Using this metric, we find that the pseudo-critical temperature converges to 920 K in the large N limit [Fig. 2(b)]. Remarkably, this critical temperature coincides with the mass-action rate law prediction of the explosion temperature, T = k5P/2k3 = 920 K, at 1 atm for an initially stoichiometric mixture.6 In the large-N limit, the traditional explosion limit temperature and mass-action kinetics should be a valid description. For reference, the inset of Fig. 2(b) shows the second explosion limit in (P, T)-space. The time associated with the critical point for finite-size systems, the pseudo-critical time tc(N), also converges to a finite value for the largest systems we simulated. This critical time is the longest time where the overall reaction completes and there is a detectable jump in activity.

Together with the decay of the coexistence width (as 1N), these results are evidence that as N the separation of active and inactive phases becomes a sharp boundary ending in a critical point (Sec. S6 of the supplementary material). The finite discontinuity in the most probable activity vanishes at temperatures approaching the putative critical point, and the onset of the discontinuity (seen as a change in the slope of ΔK) marks the critical temperature. Above the critical temperature, when there is phase coexistence, the active phase consists of trajectories with more reactions between initiations than those in the inactive phase (Fig. S5 of the supplementary material). Below the critical temperature, there is homogeneity in the types of reactions that occur over time and we cannot distinguish between active and inactive phases. In other words, at TTc or ttc, the inactive phase transforms continuously into the active phase.

At equilibrium, phase transitions are not possible for one-dimensional systems with short-range interactions. But, out of equilibrium, there are well-known examples of nonequilibrium phase transitions in one spatial dimension.57 The branching process is one example of a zero-dimensional model—one lacking spatial structure and correlations. The model here, and the nonequilibrium phase transition it exhibits, is related in that there is no spatial dependence and, therefore, no correlation length.58 There are, however, divergent quantities moving down the coexistence line marked by T* and t*. As the peaks in activity merge, approaching the critical point from above ΔKN0, and so NΔK diverges (in analogy with Δρ1 along the liquid-gas line). Another divergent quantity is the ignition time, as can be derived in the present case from the mass-action rate equations. The characteristic (branching) time, 1(2k3k5[M]), is defined in terms of the explosion limit condition. Clearly this time diverges at the T and P conditions of the explosion limit. The ignition time, which depends on the branching time, diverges logarithmically (see Sec. S8 of the supplementary material).

From observations of phase coexistence and jumps in activity, the distinctive feature of a first-order phase transition, we anticipate this discontinuous change in the order parameter to translate into an associated delta-function behavior in the susceptibility. Finite-size fluctuations in composition affect the location of the critical temperature and the size of the coexistence region. To quantify the scaling behavior of activity fluctuations, we use the dynamical susceptibility, χK(t)=K2K2, where indicates an average over trajectories. For a given T and N, the susceptibility has a maximum at a time near where the active and inactive phases have equal probability density. The finite-size scaling of χK(t) has the quantitative features of a first-order phase transition above the explosion limit temperature (Fig. 3). The maximum of the susceptibility, χmax, grows with N. Specifically, the growth is linear below Tc(N) and as N2 above. The full width at half maximum (FWHM) of this susceptibility peak scales as 1/N and the time of the maximum variance, tχ, scales as 1/N, tχ(N)=tccNβ, with β1 above the explosion limit and β0 below, with a smooth transition at intermediate temperatures. Near the critical temperature, the growth of χmax with N shows changes of slope that we attribute to the critical slowing down of the statistics in trajectory space. Overall, the scaling of the susceptibility, χK(t), identifies the phase boundary as a first-order coexistence line that terminates at a critical point in the infinite-size limit.

FIG. 3.

Finite-size scaling of the susceptibility, χK(t)=K2K2, between 860 K and 1000 K. (a) The maximum, χmax, scales as N below and N2 above 920 K. (b) The time, tχ, at which the susceptibility reaches a maximum. (c) The full width at half maximum of the susceptibility as a function of system size. All quantities obey first-order scaling above the critical temperature. (d) The scaling exponents of the maximum susceptibility, its temporal location, tχ, and the full width at half maximum, FWHM, as a function of temperature. Exponents are from linear fits of scaling data in (a)–(c). The exponents are consistent with a first-order phase transition above a critical temperature that agrees with the known explosion temperature, 920 K at 1 atm, in the thermodynamic limit. Below 920 K, χmax is linearly extensive, while tχ and the FWHM are intensive. Statistical errors are smaller than the symbol size.

FIG. 3.

Finite-size scaling of the susceptibility, χK(t)=K2K2, between 860 K and 1000 K. (a) The maximum, χmax, scales as N below and N2 above 920 K. (b) The time, tχ, at which the susceptibility reaches a maximum. (c) The full width at half maximum of the susceptibility as a function of system size. All quantities obey first-order scaling above the critical temperature. (d) The scaling exponents of the maximum susceptibility, its temporal location, tχ, and the full width at half maximum, FWHM, as a function of temperature. Exponents are from linear fits of scaling data in (a)–(c). The exponents are consistent with a first-order phase transition above a critical temperature that agrees with the known explosion temperature, 920 K at 1 atm, in the thermodynamic limit. Below 920 K, χmax is linearly extensive, while tχ and the FWHM are intensive. Statistical errors are smaller than the symbol size.

Close modal

When solving deterministic rate equations, which we do not do here, it is common to add a rate equation for temperature to model an adiabatic reactor and account for self-heating.1 Here, we instead fix temperature and pressure to show how important features of combustion, such as explosion limits, manifest from the perspective of phase transitions and critical phenomena. The finite-size scaling of activity at fixed temperature shows that the pseudo-critical temperature converges to the explosion limit temperature in the large N limit, the same temperature derived exactly from the mass-action rate equations. Having established that explosion limit temperatures are critical temperatures, it will be interesting to see how energetic fluctuations and self-heating affect the chemistry and the quantitative features of coexistence and criticality.

In summary, two phases coexist in the trajectories of igniting hydrogen-oxygen mixtures: a set of trajectories that ignites by the time t* and another set that does not. The coexistence of these phases in well-defined bands of temperature, pressure, and time exhibits the characteristics of a first-order, dynamical phase transition above a critical temperature. A surprising, and important, finding is that the critical temperature below which the dynamical phases are indistinguishable, and ignition does not occur, coincides with the explosion limit temperature. The temperature of the explosion limit is paramount to practical applications that require the safe storage of hydrogen. While fluctuations of finite systems blur the macroscopically sharp explosion limits, the finite-size scaling of activity fluctuations shows that there is a well-defined macroscopic limit with the coexistence line terminating at this critical point.

Our analysis focuses on the second explosion limit for stoichiometric hydrogen-oxygen mixtures, but it seems plausible that other explosion limits, dominated by other sets of elementary reactions and other chemical feedback mechanisms, are also critical points—especially, given that hydrogen oxidation accompanies the combustion of hydrocarbons. The universality of this view of ignition as a dynamical phase transition remains to be tested (e.g., how this framework might distinguish “strong” and “weak” ignition events near the second extended explosion limit). Overall, these results suggest that combustion kinetics might be recast as phase transitions and critical phenomena, mapping the nonequilibrium problem of ignition in the space of compositions or configurations onto an equilibrium problem in the space of trajectories. Treating the coexisting dynamical phases in a combusting mixture this way is an alternative to the existing paradigm of mean-field rate equations and opens up the possibility for further analogies with implications for practical applications.

See supplementary material for more information on the chemical reaction mechanism, derivations, finite-size scaling results at NVT, and estimation methods.

This material is based on work supported by the U.S. Army Research Laboratory and the U.S. Army Research Office under Grant No. W911NF-14-1-0359. We acknowledge the use of the supercomputing facilities managed by the Research Computing Group at the University of Massachusetts Boston as well as the University of Massachusetts Green High Performance Computing Cluster. The authors thank Jesus J. del Pozo for useful conversations and comments on the manuscript.

1.
C. K.
Law
,
Combustion Physics
, 1st ed. (
Cambridge University Press
,
2006
).
2.
C. M.
White
,
R. R.
Steeper
, and
A. E.
Lutz
,
Int. J. Hydrogen Energy
31
,
1292
(
2005
).
3.
I.
Ray
,
T.
Chakraborty
,
D.
Roy
,
A.
Datta
, and
B. K.
Mandal
,
Int. J. Emerging Technol. Adv. Eng.
3
,
119
(
2013
), available at http://citeseerx.ist.psu.edu/viewdoc/citations?doi=10.1.1.414.953.
4.
S.
Verhelst
,
Int. J. Hydrogen Energy
39
,
1071
(
2014
).
5.
S.
Verhelst
and
T.
Wallner
,
Prog. Energy Combust. Sci.
35
,
490–527
(
2009
).
6.
M. P.
Burke
,
M.
Chaos
,
Y.
Ju
,
F. L.
Dryer
, and
S. J.
Klippenstein
,
Int. J. Chem. Kinet.
44
,
444
(
2011
).
7.
D. G.
Vlachos
,
Chem. Eng. Sci.
53
,
157
(
1998
).
8.
G.
Nicolis
,
F.
Baras
, and
M. M.
Mansour
,
NATO ASI
120
,
171
(
1983
).
9.
A.
Fernández
and
H.
Rabitz
,
Ber. Bunsengesellschaft Phys. Chem.
92
,
754
(
1988
).
10.
M.
Alaghemandi
,
L. B.
Newcomb
, and
J. R.
Green
,
J. Phys. Chem. A
121
,
1686
(
2017
).
11.
K.
Chenoweth
,
A. C. T.
van Duin
, and
W. A.
Goddard
,
J. Phys. Chem. A
112
,
1040
(
2008
).
12.
S.
Agrawalla
and
A. C. T.
van Duin
,
J. Phys. Chem. A
115
,
960
(
2011
).
13.
T.
Cheng
,
A.
Jaramillo-Botero
,
W. A.
Goddard
, and
H.
Sun
,
J. Am. Chem. Soc.
136
,
9434
(
2014
).
14.
M.
Alaghemandi
and
J. R.
Green
,
Phys. Chem. Chem. Phys.
18
,
2810
(
2016
).
15.
S. B.
Nicholson
,
M.
Alaghemandi
, and
J. R.
Green
,
J. Chem. Phys.
145
,
084112
(
2016
).
16.
X.
Wang
and
C. K.
Law
,
J. Chem. Phys.
138
,
134305
(
2013
).
17.
N.
Goldenfeld
,
Lectures on Phase Transitions and the Renormalization Group
, 1st ed. (
Westview Press
,
1992
).
18.
J. M.
Yeomans
,
Statistical Mechanics of Phase Transitions
, 1st ed. (
Clarendon Press
,
1992
).
19.
20.
F.
van Wijland
,
Phys. Rev. Lett.
89
,
190602
(
2002
).
21.
J. P.
Garrahan
,
R. L.
Jack
,
V.
Lecomte
,
E.
Pitard
,
K.
Van Duijvendijk
, and
F.
Van Wijland
,
Phys. Rev. Lett.
98
,
195702
(
2007
).
22.
J. K.
Weber
,
R. L.
Jack
, and
V. S.
Pande
,
J. Am. Chem. Soc.
135
,
5501
(
2013
).
23.
C.
Jarzynski
,
Annu. Rev. Condens. Matter Phys.
2
,
329
(
2011
).
25.
S.
Pressé
,
K.
Ghosh
,
J.
Lee
, and
K. A.
Dill
,
Rev. Mod. Phys.
85
,
1115
(
2013
).
26.
H.
Touchette
and
R. J.
Harris
, “
Large deviation approach to nonequilibrium systems
,” in
Nonequilibrium Statistical Physics of Small Systems: Fluctuation Relations and Beyond
, 1st ed., edited by
R.
Klages
,
W.
Just
, and
C.
Jarzynski
(
John Wiley & Sons, Inc.
,
2013
), pp.
335
360
.
27.
M.
Frankowicz
,
M. M.
Mansour
, and
G.
Nicolis
,
Physica A
125
,
237
(
1984
).
28.
P.
Grassberger
,
Z. Phys. B: Condens. Matter
47
,
365
(
1982
).
29.
M. J.
Pilling
and
P. W.
Seakins
,
Reaction Kinetics
(
Oxford Science Publications
,
1995
).
30.
F. A.
Williams
,
J. Loss Prev. Process Ind.
21
,
131
(
2008
).
31.
P.
Bovin
,
C.
Jiménez
,
A. L.
Sánchez
, and
F. A.
Williams
,
Proc. Combust. Inst.
33
(
1
),
517
(
2010
).
32.
S. P.
Karkach
and
V. I.
Osherov
,
J. Chem. Phys.
110
,
11918
(
1999
).
33.
M.
Monge-Palacios
and
H.
Rafatijo
,
Phys. Chem. Chem. Phys.
19
,
2175
(
2017
).
34.
D. T.
Gillespie
,
J. Phys. Chem.
81
,
2340
(
1977
).
35.
36.
J. W.
Moore
and
R. G.
Pearson
,
Kinetics and Mechanism
, 3rd ed. (
John Wiley & Sons
,
New York
,
1981
), p.
480
.
37.
C. K.
Law
,
Combustion Physics
, 1st ed. (
Cambridge University Press
,
Cambridge
,
2006
).
38.
G.
Mittal
,
C. J.
Sung
, and
R. A.
Yetter
,
Int. J. Chem. Kinet.
38
,
516
(
2006
).
39.
G.
Marin
and
G. S.
Yablonsky
,
Kinetics of Chemical Reactions
, 1st ed. (
Wiley-VCH
,
Singapore
,
2011
), p.
446
.
40.
Z.
Zhao
,
Z.
Chen
, and
S.
Chen
,
Chin. Sci. Bull.
56
,
215
(
2011
).
41.
A. L.
Sánchez
,
E.
Fernández-Tarrazo
,
P.
Bovin
,
A.
Liñán
, and
F. A.
Williams
,
C. R. Mec.
340
,
882
(
2012
).
42.
Y.
Zhang
,
Z.
Huang
,
L.
Wei
,
J.
Zhang
, and
C. K.
Law
,
Combust. Flame
159
,
918
(
2012
).
43.
E.
Fernández-Tarrazo
,
A. L.
Sánchez
, and
F. A.
Williams
,
Combust. Flame
160
,
1981
(
2013
).
44.
A.
Kéromnès
,
W.
Metcalfe
,
K.
Heufer
,
N.
Donohoe
,
A.
Das
,
C.-J.
Sung
,
J.
Herzler
,
C.
Naumann
,
P.
Griebel
,
O.
Mathieu
 et al,
Combust. Flame
160
,
995
(
2013
).
45.
R. V.
de Vijver
,
N. M.
Vandewiele
,
P. L.
Bhoorasingh
,
B. L.
Slakman
,
F. S.
Khanshan
,
H.-H.
Carstensen
,
M. F.
Reyniers
,
G. B.
Marin
,
R. H.
West
, and
K. M. V.
Geem
,
Int. J. Chem. Kinet.
47
,
199
(
2015
).
46.
C. W.
Gao
,
J. W.
Allen
,
W. H.
Green
, and
R. H.
West
,
Comput. Phys. Commun.
203
,
212
(
2015
).
47.
A. K.
Burnham
,
X.
Zhou
, and
L. J.
Broadbelt
,
Energy Fuels
29
,
2906
(
2015
).
48.
N.
Donohoe
,
K. A.
Heufer
,
C. J.
Aul
,
E. L.
Petersen
,
G.
Bourque
,
R.
Gordon
, and
H. J.
Curran
,
Combust. Flame
162
,
1126
(
2015
).
49.
T.
Varga
,
T.
Nagy
,
C.
Olm
,
I. G.
Zsély
,
R.
Pálvölgyi
,
É.
Valkó
,
G.
Vincze
,
M.
Cserháti
,
H. J.
Curran
, and
T.
Turányi
,
Proc. Combust. Inst.
35
,
589
(
2015
).
50.
A. L.
Wagner
,
P. E.
Yelvington
,
J.
Cai
, and
W. H.
Green
,
J. Propul. Power
33
,
350
(
2017
).
51.
C.
Beck
and
F.
Schlögl
,
Thermodynamics of Chaotic Systems: An Introduction
(
Cambridge University Press
,
1995
).
52.
C.
Borgs
and
R.
Kotecký
,
J. Stat. Phys.
61
,
79
(
1990
).
53.
K.
Binder
,
Phys. Rev. Lett.
47
,
693
(
1981
).
54.
K.
Binder
and
D. P.
Landau
,
Phys. Rev. B
30
,
1477
(
1984
).
55.
J. P.
Neirotti
,
P.
Serra
, and
S.
Kais
,
Phys. Rev. Lett.
79
,
3142
(
1997
).
56.
P.
Serra
,
J. P.
Neirotti
, and
S.
Kais
,
Phys. Rev. Lett.
80
,
5293
(
1998
).
57.
V.
Privman
,
Nonequilibrium Statistical Mechanics in One Dimension
(
Cambridge University Press
,
Cambridge, MA
,
1997
).
58.
J.
Marro
and
R.
Dickman
,
Nonequilibrium Phase Transitions in Lattice Models
(
Cambridge University Press
,
Cambridge, MA
,
1999
).

Supplementary Material