Stable chaos and delayed onset of statisticality in unimolecular dissociation reactions

Statistical models provide a powerful and useful class of approximations for calculating reaction rates by bypassing the need for detailed, and often difficult, dynamical considerations. Such approaches invariably invoke specific assumptions about the extent of intramolecular vibrational energy flow in the system. However, the nature of the transition to the statistical regime as a function of the molecular parameters is far from being completely understood. Here, we use tools from nonlinear dynamics to study the transition to statisticality in a model unimolecular reaction by explicitly visualizing the high dimensional classical phase space. We identify generic features in the phase space involving the intersection of two or more independent anharmonic resonances and show that the presence of correlated, but chaotic, intramolecular dynamics near such junctions leads to nonstatisticality. Interestingly, akin to the stability of asteroids in the Solar System, molecules can stay protected from dissociation at the junctions for several picoseconds due to the phenomenon of stable chaos.

ii ) 1/2 . Along q i (i = 1, 2) direction, q i ∈ [q i,min , q i,max ], where q i,min is obtained from the relation, V i (q i,min ) = E and q i,max is set to 5.5 au (considered maximum range). The allowed values for q 3 are then obtained from the constraint of energy conservation H ε (p, q 1 , q 2 ; q 3 ) = E. For the survival probability computations shown in Fig. 1b in the main text, a slice of the allowed energy shell is sampled as follows. First we choose a specific angle section θ θ θ s ≡ (θ 1s , θ 2s , θ 3s ). The allowed range of J i (i = 1, 2) is, J i ∈ [0, J i,max ], where J i,max = 2D i /ω i . The allowed values for J 3 are then obtained from H ε (θ θ θ s , J 1 , J 2 ; J 3 ) = E.
Suitably large grids along each of the phase space variables ({q, p} or {θ θ θ , J}) are placed and the energy shell is sampled uniformly. The grid sizes N q,p 1,2 along q i and p i (i = 1, 2) are chosen to be the same and along p 3 it is set to N p 3 . For the ensemble used in the results shown in the main text (Fig. 1a), (N q,p 1,2 , N p 3 ) = (20, 28). For Fig. 1b in the main text, the grid sizes are (N J 1 , N J 2 ) = (500, 500). The size of different ensembles for different ε's used in the main text for survival probability computation is listed in Supplementary Table 2. In addition, Supplementary Table 2 also provides information on the number of trajectories analyzed using wavelet timefrequency analysis (Fig. 4 of the main text). The initial conditions were then used to integrate the Hamilton's equations of motion for T = 40 ps using a 8 th order Runge-Kutta method 1 . The total energy of the trajectories were found to be conserved with sufficient accuracy. For example, the energy of a typical non-dissociated highly chaotic trajectory is conserved up to 8 th decimal place for the entire duration of the integration. The size of the grids were varied to create different ensemble sizes to check for convergence. Supplementary Fig. 1 shows the survival probability to be nearly converged. Slight variations at long times and large values of ε are due to the very small number of surviving trajectories. Note that the value of q i,max up to 7.5 au has been considered in order to check for convergence of the survival probabilities. The survival probabilities, as shown in Supplementary Fig. 1 (cyan diamonds), show slightly faster decay at initial times for small ε. However, the differences at later times are not significant. Moreover, the slopes of the survival probabilities do not depend much on the chosen value of q i,max . Therefore, for all the results shown in this work, the value of q i,max is set to 5.5 au.

Multi-exponential fits to the survival probability.
The non-exponential behaviour of the survival probabilities in Fig. 1 in the main text are fit to a general bi-exponential function 2 of the form, where ∑ 2 i=1 f i = 1 due to normalization. N(τ) is the number of undissociated trajectories at t = τ and N(0) is the size of the dissociating ensemble i.e., the total number of trajectories which dissociate within the integration time period of T = 40 ps. We have used MATLAB 3 to fit the numerical data up to 30 ps.
The resulting fits for different ε values are shown in Supplementary Table 3. In the table we also show the result of the short time (1 ps) single exponential fits for the survival probabilities shown in Fig. 1a of the main text. In Supplementary Table 4 the values of the fit parameters for Fig. 1b of the main text are shown. Note that in both cases a single exponential with rate constant k 1 ∼ 1.7 ps −1 is sufficient for ε = 1.0.
The accuracy of the fits are indicated in Supplementary Fig. 2a-b in terms of the residual defined as where S(τ) is computed from the original lifetime data set andS(τ) corresponds to the fitted data. Notice that in Supplementary Fig. 2a, the residual is very large for ε = 1.0 at longer times. This is expected since the fit is to a single exponential form and the number of trajectories contributing to the origin of the second time scale is extremely small. Table 3. Results of the fit to the survival probabilities in Fig. 1a of the main text. indicator (FLI) method. The FLI has been used to map the Arnold web for various dynamical systems and details can be found in the literature 5,6 . Here we provide a brief introduction to computing the FLI. For a general dynamical system,

Supplementary
with X = {x 1 , x 2 , . . . , x i , . . .} being the dynamical variables, the FLI 7 is defined as where, v is a tangent vector and However for technical reasons the alternative definition 5 is usually preferred. Consequently, in this work we use Supplementary Equation 6 to compute the FLI. The variation of FLI with time for regular and chaotic trajectories is distinctly different and hence different types of dynamical behaviour can be identified, usually within timescales t f much shorter compared to those required for computing the Lyapunov exponents. In Supplementary Fig. 3 we show an example for the behaviour of FLI for different types of trajectories. For regular trajectories, the FLI saturates quickly whereas for the chaotic trajectory the FLI increases linearly. As can be seen from Supplementary Fig. 3, the distinction between chaotic and regular trajectories is already clear by ∼ 10 ps. For comparison, in Supplementary  Fig. 3b-d we also show the variations in the mode frequencies with time for the three different trajectories using wavelet-based time-frequency analysis (see below). For the regular trajectories, the mode frequencies are nearly constant. In contrast, for the chaotic trajectory the mode frequencies exhibit drifts and rapid fluctuations. Note that even for the chaotic trajectory the frequency of the bending mode, being modeled as a harmonic oscillator, does not fluctuate as much as the stretching frequencies.
In this work we compute the Arnold web in the action space as follows. For every initial condition on a specific slice of the action-angle space, as described in the previous section, we determine the associated FLI value using Supplementary Equation 6 with a final time of T = 40 ps. However, as discussed below, for trajectories that dissociate before T = 40 ps the total integration time T is set by the lifetime t lt and hence the FLI value is noted at t lt . Note that an arbitrary choice of the initial tangent vector can lead to spurious structures 8 on the web. To avoid this, we follow the prescription by Barrio

Time-frequency analysis of the intramolecular dynamics.
To compute the local frequencies we use the wavelet-based time-frequency approach 9 of Arevalo and Wiggins. In this approach local frequencies are determined by a continuous wavelet transform of the time series corresponding to the trajectories. Specifically, the wavelet transform provides the joint time-frequency information of a signal or time series. For a system in the near-integrable regime, according to the Kolmogorov-Arnold-Moser (KAM) theory, almost all of the unperturbed tori remain intact and the frequencies associated with the KAM tori can be approximated by Ω Ω Ω 0 (J) ≡ ∂ H 0 (J)/∂ J. However with increasing the perturbation, the KAM tori are destroyed and both the regular and chaotic regions coexist in the phase space. In such cases, the frequencies associated with a chaotic trajectory can not be approximated by Ω ≈ Ω 0 . Although, for a chaotic trajectory there are several frequency components ( Supplementary Fig. 11), the joint time-frequency analysis can be used to gain information on the dominant frequencies in an appropriate time-window. Several studies [9][10][11] have established that the dominant frequencies can be associated with the chaotic trajectory visiting specific set of phase space structures. In the molecular context, the dominant frequencies can be correlated to specific intramolecular motions with certain vibrational modes being more active than others in a given time interval. We give here a brief description of the time-frequency analysis since details can be found in previous works 9, 12 .
The time series z k (t) = q k (t) + ip k (t), obtained from the dynamics, for the three modes (k = 1, 2, 3) are subjected to a continuous wavelet transform where, a > 0 and b is real and is the Morlet-Grossman mother wavelet. λ and σ are the parameters which can be tuned for better accuracy. For our calculation, we have used λ = 1 and σ = 2. Other values for these parameters yield similar results with slightly differing accuracy, depending on whether the mode of interest is a low or high frequency one. Using Supplementary Equation 7 the local frequency associated with z k (t) over a small interval of time around t = b, related to the inverse of the scale factor a, can be determined. There are two ways of analyzing the transform. First is to extract the dominant frequency at a given time by computing the maximum of the modulus of the wavelet transform i.e., max a |L ψ z z z(a, b)|. The ratios of the dominant frequencies can then be used to analyze the dynamics in the frequency ratio space (FRS). A second more comprehensive way 11 is to compute the full scalogram |L ψ z k (a, b)| 2 to look at the ridge maxima and dominant frequency variations over time. Although we have examined the scalograms of several individual trajectories, in this work we have considered only the approach based on the dominant frequencies. The reason being that working with the scalogram for the large ensemble sizes considered in this work is rather impractical.
As a check for the accuracy of the wavelet transform, Supplementary Fig. 5 shows a representative example in the case when the system is uncoupled and the total energy of E = 25 kcal mol −1 is partitioned into the three modes. The energy of the each mode as a function of time is shown in Supplementary Fig. 5a. Note that due to the modes being uncoupled, there is no IVR and the individual mode energies are conserved. The frequencies computed using the wavelet analysis for this example are shown in Supplementary Fig. 5b. The calculated frequencies, using Supplementary Equation 17, of the Morse oscillators are 1003.9 cm −1 and 828.5 cm −1 respectively, which are very close to the numerically computed frequencies (1006.7 cm −1 , 830.6 cm −1 ). In the presence of the resonant coupling terms, the system becomes non-integrable and the frequencies show nontrivial variations with time. The frequency variation depends on the nature of the trajectory. In Supplementary Fig. 3d and Supplementary Fig. 11, examples for such frequency variations can be seen.

Supplementary Note 1
A complication in interpreting the FLI results for unbounded systems, as studied in the present work, is worth noting at this point. For bounded systems a fixed, appropriately chosen time T is sufficient to discriminate between different types of dynamics. However, in our case the trajectories can dissociate with different lifetimes t lt . Thus, for trajectories with very short ( 5 ps) lifetimes the FLI values are expected to be small and hence indistinguishable from regular trajectories. In other words, the dissociation event sets the total integration timescale for each trajectory. However, note that such ambiguities occur only for trajectories with rather short lifetimes. As examples for such behaviour, in Supplementary Fig. 7 we show the comparison of FLI for different trajectories with varying lifetimes. For trajectories with small t lt ( Supplementary  Fig. 7a), the FLI is very close to that of a regular trajectory. However, the value of t lt helps in this regard to distinguish these two class of trajectories because of very large t lt for regular trajectories. For relatively long t lt (Supplementary Fig. 7b), the FLI is clearly capable of distinguishing between the different dynamical behaviours.

Supplementary Note 2
Here we provide a summary of the action-angle transformations and strengths of the nonlinear resonances extracted from the analysis. For the bending mode the harmonic oscillator Hamiltonian used in the main text is of the where, q 3 and p 3 are the conjugate coordinate and momentum and ω 3 is the frequency. The action-angle transformation is given by, where J 3 and θ 3 are the action and angle variables respectively. In the action-angle representation the transformed Hamiltonian can be expressed as As expected, the frequency ∂ H with α i and D i (i = 1, 2) being the Morse length scale parameter and the dissociation energy respectively. The action-angle transformation 13 can be expressed as ii D i and the action variable J i . In the action-angle representation the Morse Hamiltonian transforms to with the energy dependent nonlinear frequency

11/25
Based on the transformations above, the Hamiltonian in the main text (eq. 3,4), H(J, θ ) ≡ H 0 (J) + εV (J, θ θ θ ) has an integrable zeroth-order part and the coupling term which, using the known Fourier expansion 13 of the momentum p j = , can be explicitly written as a Fourier expansion where we have denoted The coefficients in Supplementary Equation 19 can now be explicitly expressed as f (12) l,k (J 1 , for i = 1, 2 and B (3) (J 3 ) = 2ω 3 J 3 /G (0) 33 .

Supplementary Note 3
Another way of visualizing a subspace of the phase space is by constructing the zero momentum surface (ZMS) 14 . The ZMS is defined as the surface with p i = 0, for all i. To map the ZMS, a 500 × 500 grid is constructed along q 1 − q 2 keeping p 1 = p 2 = p 3 = 0. The allowed value of q 3 is then obtained from energy conservation. All the trajectories with initial conditions on the ZMS are integrated for 40 ps, using the same integrator mentioned above, and the resulting FLI is projected on the q 1 − q 2 plane. In Supplementary Fig. 8 we show the ZMS for varying ε and the corresponding lifetime plots. Note that the different regular regions in the ZMS are clearly identified. Specifically, on comparing to the Arnold webs shown in Fig. 2 of the main text, specific structures on the ZMS can be related to the different resonance and the junctions formed by them. The corresponding lifetime plots also shows one-to-one correspondence with the ZMS features.
Note that the results in Supplementary Fig. 8 indicates persistence of regular regions near the potential minima with undissociated trajectories up until ε ∼ 1. However, the survival probability results shown in Fig. 1 of the main text do not include such initial conditions. Thus, the second longer timescales in the survival probabilities are not due to these regions. Instead, the long lifetime trajectories observed near the periphery of these regions and the partially broken 1:1 resonance regions give rise to the tail in the lifetime distributions. More importantly, and related to the results shown in Fig. 4 of the main text, the observation of significantly large lifetimes near the M 1,0,−1 0,1,−1 junction regions at q 1 ≈ q 2 ∼ 3.1 in Supplementary Fig. 8 are key to the non-statisticality. We also point out that the regular region seen at (q 1 , q 2 ) ≈ (2.1, 3.2) in Supplementary Fig. 8 for ε = 0.4 does not manifest on the Arnold web slice shown in Fig. 2 of the main text. Clearly, a combination of the (J 1 , J 2 ) Arnold web slice and the ZMS yields fairly detailed information regarding the phase space structures that lead to strongly correlated intramolecular dynamics.

13/25 6 Supplementary Note 4
Here we briefly sketch possible estimates of the critical ε for the onset of widespread chaos using the idea based on the Chirikov criterion 15 .
Thinking along the lines of Oxtoby and Rice 16 , one can estimate the onset of extensive overlap and hence RRKM behavior can be predicted from the Arnold webs shown in the main paper. Note that a straightforward generalization of the overlap criterion to the present case is not obvious. In f = 2 the zeroth order energy surface is a line in the (J 1 , J 2 ) action space. The various resonances puncture this zeroth order energy surface and the size of the holes created scales with the coupling strength. In our case the zeroth-order energy surface is in the (J 1 , J 2 ) plane with the various resonances forming a dense network. Thus, it is not immediately apparent as to what set of resonances need to be considered for the estimate of critical ε. A natural choice would be the resonances of lowest order, since they have the largest widths. Still, which set of lowest order resonances is far from obvious. In the following, we use the fact that our computations indicate the 1:1 stretch-bend resonances as the prominent ones to consider.
We now sketch a calculation to obtain an approximate estimate. Note that we follow the notation in Lichtenberg and Lieberman 17 . Consider the Bunker Hamiltonian expressed in actionangle coordinates and focus on the resonance mθ 1 =θ 3 i.e., mΩ 1 (J) = Ω 3 (J). This can be done by performing a canonical transformation (θ θ θ , J) → (θ θ θ ,Ĵ) via the following generating function Transforming the original Hamiltonian using the identities and averaging overθ 3 we obtain a Hamiltonian corresponding to an integrable Hamiltonian that is relevant to the resonance of interest. In the above, the coefficients g (13) m are known from the Morse action-angle transformation and the zeroth-order Hamiltonian is given bȳ Since the averaged Hamiltonian is ignorable inθ 2 andθ 3 , the conjugate actionsĴ 2 = J 2 and J 3 = (J 1 /m) + J 3 (polyad for the specific resonance) are conserved. As usual, one now expands the above Hamiltonian about the fixed points to obtain a pendulum like Hamiltonian which has a full widthŴ = 4(|F/G|) 1/2 . In the original variables the full width of the mΩ 1 (J) = Ω 3 (J) resonance is estimated as

14/25
where we have denoted In the above expression r ± ≡ 1 ± (ω 3 /mω 1 ), α 1 is the exponent in the Morse oscillator potential function, and we have taken J 2 = 0 to estimate the value of J r 3 using the zeroth-order Hamiltonian. A very crude estimate for the threshold coupling strength ε by invoking only the dominant 1:1 stretch-bend resonance is to set the full width of the resonance to be nearly the entire bound range of the action J 1 with J 2 = 0. In other words, using parameters appropriate for Bunker's model # 6 and at E = 34 kcal mol −1 , we set W 1:1 ≈ 16.3 √ ε ∼ 13. This yields the threshold value to be ε max ∼ 0.64. There are two factors worth mentioning at this stage. Firstly, Chirikov analysis is independent of the angles i.e., the estimate just made is independent of the specific angle slice being chosen. For such large values of ε max this is not correct. Secondly, it is well known that the critical ε value estimated using a Chirikov analysis 18 is not very accurate. In fact, typically, the actual values obtained from dynamics are about half that of the Chirikov estimate. However, this estimate can be made sharper and more accurate in a number of ways using the Greene residue theorem 19, 20 or the renormalization approach 18, 21 . We do not undertake this task here since the improved estimates will only reduce the ε max value. Thus, as argued next, a more reasonable estimate is ε max < 0.32.
From the analysis of the zeroth-order Hamiltonian we estimate the resonant actions for the 1:1, 2:1, and 3:1 stretch bend-resonances to be around J r 1 ∼ 6.5, 10.8, and 12.2 respectively. These resonances can be clearly identified in the Arnold web shown in Supplementary Fig. 4b. Denoting the resonant actions by J r m:1 , the condition for the overlap of the prominent 1:1 with the higher order stretch bend resonances can be written as In the above we have invoked the so called "two-thirds rule" 22 for heuristically correcting the simple Chirikov estimate for the role of higher order resonances. Using the above criterion with values appropriate for Bunker's model # 6 at E = 34 kcal mol −1 , we find ε ∼ 0.05 and 0.1 for m = 2 and 3 respectively. Comparing to the webs computed using the FLI approach, the estimates are quite good. Thus, large scale chaos is predicted, and confirmed by computations, to occur already by ε ∼ 0.2 − 0.3. Figure 10. As in Supplementary Fig. 9. Here the angles θ 2 and θ 3 are fixed at 0 and θ 1 is varied. (a) 0 (b) π/5 (c) 2π/5 (d) 3π/5 (e) 4π/5 (f) π.

17/25 9 Supplementary Note 5
In the main text the trapping due to the M 1,0,−1 0,1,−1 junction were highlighted. Nevertheless, the dynamics for other intermediate, but beyond the Chirikov overlap regime, values of ε exhibits trapping due to several other junctions. Here we provide some examples of such trapped intramolecular motion. In Supplementary Fig. 11 and Supplementary Fig. 12 we show the time-frequency analysis and the action space projection of the example long lifetime trajectories for ε = 0.2 and 0.3. Note that the action space projections are shown on a specific angle slice for convenience. In particular, the various actions at different times do not correspond to the initial angle slice. For example, the trajectory in Supplementary Fig. 11a  As mentioned in the main text, trajectories trapped at ( f 1 , f 2 ) junction in the FRS obey an approximately conserved polyad Note that P d is conserved i.e., the Poisson bracket {H, P d } ≈ 0 for the duration that the trajectory is trapped near that specific junction. The same trajectory can then further get trapped near different junctions before dissociating. In Supplementary Fig. 12, the approximately conserved polyads are shown. Note that in Supplementary Fig. 12B, the correlated dynamics near M 1,0,−1 0,1,−1 is established by the fact that the quantities J 1 + J 3 and J 2 + J 3 are not conserved in comparison to the full polyad.  Supplementary Fig. 11. Different trapping events near junctions are shown in cyan and red (as indicated in Supplementary Fig. 11). The portion of the trajectories shortly before dissociation are shown in green. Panels A and B show the variation of the zeroth order actions corresponding to the trajectories (a) and (d) respectively. J 1 (brown), J 2 (magenta) and J 3 (orange). Also shown are the appropriate approximately conserved polyads for trajectories A and B. The polyads (shown in blue) are (J 1 + 3J 2 /2 + J 3 ) for A and (J 1 + J 2 + J 3 ) for B. Shown in green is the variation of (J 1 + J 3 ) for both cases, while the violet line for B shows the variation of (J 2 + J 3 ). The time period of polyad conservation is indicated by a cyan arrow.

Supplementary Note 6
The FRS density map in Fig. 4 of the main text indicates enhanced density near various junctions. For quantitatively determining the extent of locking near the junctions, we perform further analysis of the time-frequency data. The analysis is not straightforward since the resonances, and hence the junctions, are dense on the Arnold web. Strictly speaking, the width of a given resonance r · Ω Ω Ω = 0 of order O r = |r 1 | + |r 2 | + |r 3 | with r ≡ (r 1 , r 2 , r 3 ) is given by Thus, at the intersection of two independent resonances of orders O r and O r one needs to consider a W r × W r rectangle in order to determine the trapping time statistics. Note that the size of the rectangular region also scales with the coupling strength ε. This approach, however, is not very practical since it requires a careful consideration of the total order, related to the timescales of dynamical interest, up to which resonances need to be considered. At the same time, for junctions involving large order resonances one requires significantly larger number of trajectories to obtain reliable statistics. In this work we take a simple, but an approximate approach since we are interested in demonstrating that the trajectories spend time trapped in the vicinity of the junctions. In particular, we do not attempt here to do a detailed characterization of the trapped dynamics near junctions in terms of the crossing between specific order resonances. We have considered several key resonance junctions on the basis of number of visitations of trajectories in FRS density map 19/25 Supplementary Table 5. Key resonance junctions ( f 1 , f 2 ) considered (indicated by " ") to compute locking time distribution (Fig. 4e-f in the main text). and demarcated a square region around them with a tolerance ( f 1 ± δ , f 2 ± δ ) with δ = 0.125. In Supplementary Fig. 6, we show junctions (1, 3/2), (1, 1), and (3/2, 1) as examples with the demarcated square region. The full list of the junctions considered in this work are shown in Supplementary Table 5. Note that too tiny a δ would not provide enough statistics and too large a value may include several nearby junctions. However, we did the calculations with a reasonably smaller δ value, and the results do not differ qualitatively.
The model used in the current study uses the equilibrium G-matrix elements. This then allows us to make a detailed action-angle analysis for the system. For the original Bunker model, however, the action-angle transformations are not available. Therefore, an important issue is whether our results are generic enough. We believe there are reasons, apart from the fact that resonance junctions are ubiquitous in any nonlinear Hamiltonian system with f ≥ 3, to expect the influence of the resonance junctions in a broader class of systems. For example, previous studies on the OCS molecule, with a Hamiltonian essentially of the form considered by Bunker, do indicate the presence of junctions that influence the IVR dynamics. It would be interesting to study the dissociation dynamics of OCS and link the deviations from RRKM behaviour with the extent of trapping near junctions. In addition, our preliminary studies on some of the other Bunker models suggest that the stabilisation near resonance junctions should be a key factor in understanding the origins of non-statistical reaction dynamics in more general systems as well. The purpose of this note is to summarize our preliminary results on a few of the other Bunker models and the full Bunker case for model 6.
• For reasons noted above, for the full Bunker case we cannot present the Arnold webs in the action space. However, the ZMS representation (Supplementary Note 3) can still be utilized and we expect that any indication of long lifetimes due to the influence of the resonance junctions should be discernible. Thus, the ZMS of the full Bunker model 23 with coordinate dependent G-matrix elements is shown in Supplementary Fig. 14a along with the corresponding lifetimes at ε = 0.1 and energy, E = 34 kcal mol −1 . Note that the nontrivial coupling due to G 33 (q) makes the system nonintegrable even at ε = 0. In particular, the ε values cannot be compared for the equilibrium G-matrix case and the coordinate dependent G-matrix case. Nevertheless, the wavelet time-frequency results shown in Supplementary Fig. 14 indicate the influence of resonance junctions in the system. One can clearly observe trapping near M 1,0,−2 0,1,−1 junction and M 1,0,−1 0,1,−1 junction. Note that for the full system as well the M 1,0,−1 0,1,−1 junction is important.
• Here we present the web for two other Bunker models 24 , model 8A and 8D. For both the cases, the equilibrium value of angle, q 0 3 = π radians leads to G (0) i3 = 0 (i = 1, 2), i.e. the system becomes an effectively two degrees of freedom system. We therefore artificially choose q 0 3 = 2.0 radians, resulting in three degrees of freedom models. The artificial value of q 0 3 can be thought of as the case wherein the molecule is in some excited electronic state with a bent geometry. Both the system have been studied at a total energy E = 70 kcal mol −1 . At this energy, the second stretching mode dissociates (dissociation energy, D 2 = 60 kcal mol −1 ). We note that the original versions of these models were shown to be non-RRKM by Bunker.
The evolution of the web with increasing coupling strength ε for the two model are shown in Supplementary Fig. 15a-c. Although the webs indicate the appearance of widespread chaos with increasing ε, note that a prominent junction around (J 1 , J 2 ) ≈ (1.5, 7.0) (for model 8A) and (J 1 , J 2 ) ≈ (1.0, 13.0) (for model 8D), indicated by arrows, persist for fairly large coupling strengths. The lifetime data for ε = 1.0 confirm that the junctions continue to influence the dissociation dynamics. An additional junction structure, with long lifetimes, can also be observed for model 8D. The angle slice is (θ 1 , θ 2 , θ 3 ) = (π/2, π/2, 0).