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.

## INTRODUCTION

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 simulations^{10} 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 transitions^{19,20} are accumulating interest recently in statistical physics, due in part to their connection to the glass transition^{21} 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 observables^{23} 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.

## MODEL AND SIMULATION METHODS

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, 2H_{2} + O_{2} $\u2192$ 2H_{2}O. 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 H_{2}O and HO_{2}. 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

where M is an arbitrary collision partner. For each reaction, R_{i}, we model the temperature dependence of the rate constants, *k*_{i}, 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 directly^{32} but could be a composite of two reactions included in some full models near the second explosion limit: the initiation reactions H_{2} + O_{2} $\u21cb$ H + HO_{2} and H + HO_{2} $\u21cb$ 2OH.^{6} Termolecular reactions have been proposed as an alternative means of initiation (e.g., 2O_{2} + H_{2} $\u21cb$ 2HO_{2} and $2HO2\u21cb$ 2OH + O_{2}^{32}) 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 H_{2}/O_{2} 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, $\tau $, between reactions, R_{i}, is an exponentially distributed random variable with a probability $p(\tau ,i)=aie\u2212a0\tau d\tau $. This probability depends on the mixture composition of (or number of molecules of each) chemical species, $X(t)={X1,X2,\u2026,XM}$, through the propensity of a single reaction, *a*_{i}, and the total propensity, $a0=\u2211iai$. 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 H_{2}O, 2.86 for H_{2}, and 1.00 for both O_{2} and HO_{2}.

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 H_{2} and O_{2}, and the final composition is a mixture of H_{2}O and HO_{2}, 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, 2*k*_{3} = *k*_{5}[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, 2*k*_{3} = *k*_{5}*P*/*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\u2192\u221e$ limit, mixtures can spontaneously ignite when heated to temperatures above 920 K.

## RESULTS AND DISCUSSION

### Signatures of a nonequilibrium phase transition

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 HO_{2}, 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 $0\u2264t\u2264tobs$, $K[x(tobs)]=\u2211i1(\Delta ti)$. The indicator function $1(\Delta ti)=1$ if a reaction event occurs and is zero otherwise. We divide *t*_{obs} uniformly into intervals, $\Delta 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 H_{2}O nor HO_{2} 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\u2192\u221e$ 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, *t*_{obs} = *t*_{i}. We define the final time as the time at which the overall reaction is complete, *t*_{obs} = *t*_{f}. As a consequence of these temporal boundary conditions, *K*(*t*_{i}) = 0 and the activity density is $K(tf)\u2215N\u223c1$ (for this model). For the 10^{6} 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 10^{6} 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).

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\kappa a$ and a slow burning “inactive” phase associated with the lower activity peak located at $Ki=N\kappa 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 H_{2}/O_{2} 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, HO_{2} and H_{2}O. Those undergoing low activity trajectories are largely composed of reactants, H_{2} and O_{2}. 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, *K*_{mp}/*N*. The resulting phase diagram in Fig. 1(b) shows a jump from $\kappa i$ to $\kappa 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 *K*_{mp}/*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 $1\u2215N$ 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, $\Delta K=Ka\u2212Ki$ 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, *T*_{c}(*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 mechanics^{53,54} and quantum mechanics.^{55,56}

### Finite-size scaling

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, $\Delta 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 H_{2} and O_{2}. 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.

To determine the critical temperature, $Tc(N\u2192\u221e)=Tc$, we again use the finite-size scaling of the “latent activity” of ignition, $\Delta K$, at the time $t*$. For a given system size, we identify *T*_{c}(*N*) with the temperature at which $\Delta $*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* = *k*_{5}*P*/2*k*_{3} = 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 *t*_{c}(*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 $1\u2215N$), these results are evidence that as $N\u2192\u221e$ 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 $\Delta 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 $T\u2264Tc$ or $t\u2265tc$, 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 $\Delta K\u2215N\u21920$, and so $N\u2215\Delta K$ diverges (in analogy with $\Delta \rho \u22121$ 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, $\u221d1\u2215(2k3\u2212k5[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, $\chi K(t)=\u27e8K2\u27e9\u2212K2$, where $\u2026$ 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 $\chi K(t)$ has the quantitative features of a first-order phase transition above the explosion limit temperature (Fig. 3). The maximum of the susceptibility, $\chi max$, grows with *N*. Specifically, the growth is linear below *T*_{c}(*N*) and as $\u223cN2$ above. The full width at half maximum (FWHM) of this susceptibility peak scales as 1/*N* and the time of the maximum variance, $t\chi $, scales as 1/*N*, $t\chi (N)=tc\u2212cN\beta $, with $\beta \u2248\u22121$ above the explosion limit and $\beta \u22480$ below, with a smooth transition at intermediate temperatures. Near the critical temperature, the growth of $\chi 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, $\chi K(t)$, identifies the phase boundary as a first-order coexistence line that terminates at a critical point in the infinite-size limit.

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.

## CONCLUSIONS

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.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.