Revisiting the bonding evolution theory: a fresh perspective on the ammonia pyramidal inversion and bond dissociations in ethane and borazane

This work oﬀers a comprehensive and fresh perspective on the bonding evolution theory (BET) framework, originally proposed by Silvi and collaborators [X. Krokidis, S. Noury and B. Silvi, Characterization of elementary chemical processes by catastrophe theory, J. Phys. Chem. A , 1997 , 101, 7277–7282]. By underscoring Thom’s foundational work, we identify the parametric function characterizing bonding events along a reaction pathway through a three-step sequence to establish such association rigorously, namely: (a) computing the determinant of the Hessian matrix at all potentially degenerate critical points, (b) computing the relative distance between these points, and (c) assigning the unfolding based on these computations and considering the maximum number of critical points for each unfolding. In-depth examination of the ammonia inversion and the dissociation of ethane and ammonia borane molecules yields a striking discovery: no elliptic umbilic flag is detected along the reactive coordinate for any of the systems, contradicting previous reports. Our findings indicate that the core mechanisms of these chemical reactions can be understood using only two folds, the simplest polynomial of Thom’s theory, leading to considerable simplification. In contrast to previous reports, no signatures of the elliptic umbilic unfolding were detected in any of the systems examined. This finding dramatically simplifies the topological rationalization of electron rearrangements within the BET framework, opening new approaches for investigating complex reactions.


Introduction
The concept of a chemical bond has long been intensely debated.On the one hand, it encapsulates our understanding of matter's existence and stability.2][3][4] This dichotomy has led to various methodologies to decipher this elusive construct.7][8] Fifteen years later, Pauling 9 translated Lewis sharing idea into the emerging quantum mechanics language, building on Heitler and London's work. 10Bader and Henneker 11 used the electric field to explore a completely different interaction -the ionic bond.Gadre and co-workers 12 combined this field with topological concepts to describe various bonding situations in neutral and electrically charged molecular systems.][15][16][17][18][19][20] Bader's seminal work on the topological analysis of the electron density 21 introduced a new way to gain insights into the characterization of the chemical structure and bonding concepts.This method has also been applied to other local functions, including the electron localization function (ELF), 22 the nuclear potential, 23 the virial field, 24 the magnetically induced molecular current distribution, 25 and the spin density. 26Such type of analysis has been proposed to be classified as quantum chemical topology (QCT). 27QCT leverages the topographical analysis of the gradient vector field generated by scalar functions to extract valuable chemical information without explicitly mentioning the reference state.These approaches enable a definition of topological domains using partitioning schemes, which are parameter-free. 28Krokidis, Noury, and Silvi 29 explored for the first time the application of Thom's catastrophe theory (CT) 30 to the analysis of topographical changes of ELF 22 along a given coordinate, i.e., defining the so-called bonding evolution theory (BET) methodology.BET has been applied in an ever-increasing fashion over the last three decades to gain a deeper understanding of different covalent bonding situations, 29 C-H bonds activation, 31 proton transfer reactions, 32 cycloadditions, 2,[33][34][35][36][37] phase transition of the 14 group elements, 38 hemiaminals formation, 39,40 computing local and global properties of solids, 41 and deducing linear models for predicting activation enthalpies of organic and organometallic systems. 42he BET critical idea is modeling any reaction system as a gradient dynamical system (GDS) using the ELF as its potential.The gradient nullity condition applied to ELF yields four types of equilibria since this smooth function is defined on R 3 : attractors (ELF maxima, hereafter labeled as a), saddle points of index one and two (hereafter designated as s 1 and s 2 , respectively), and repellors (ELF minima, hereafter denoted as m).Following the changes of these solutions along a reaction path (e.g., the IRC 43 ), a rationalization of electronic events in terms of the so-called universal unfoldings 30 is carried out.5][46] Each function describes one of the seven types of degeneracy exhibited by ELF critical points (CPs).A CP is said to be degenerate if evaluating the determinant of the Hessian matrix yields 0; otherwise, it is called hyperbolic. 44,46,479][50][51][52][53] Poincare ´-Hopf theorem 54 imposes a strong constrain for all these changes since n 0 À n 1 + n 2 À n 3 = 1, where n 0 , n 1 , n 2 , and n 3 stand for the number of attractors, saddles of index one, saddles of index two, and repellors, respectively.Note that the subscript coincides with the index of the CP.This alternating sum yields zero in the context of periodic systems, such as crystalline structures.Herein we will follow the convention that the index represents the number of positive eigenvalues of the Hessian matrix at the point. 55Therefore, a CP with any positive eigenvalue is classified as a maximum of the function.In contrast, the critical point is indicative of a minimum on the potential landscape in cases where all its three eigenvalues are positive.The other two possible scenarios (on R 3 ) involving combinations of one and two positive eigenvalues lead to saddle points of indices one and two, respectively.On the other hand, there is a (big) set of points where the ELF gradient does not vanish, and more importantly, a departing trajectory of such points from the neighborhood of a Morse-type CP will never return.These intriguing points are the so-called wandering points. 45It is essential to emphasize that the BET has been applied to the study of a wide variety of chemical reactions ranging from the gas phase to the solid state; 56 however, despite its popularity, the identifying-assigning process of the correct unfolding to describe the formation/scissions of chemical bonds is not well understood.
Our group has described the electron reorganizations along a reactive coordinate in the ground 57,58 and electronically excited states [59][60][61][62][63] following the determinant of Hessian at all potentially degenerate CPs and their relative distance. 56This procedure guarantees a rigorous application of BET proposal by recovering the fundamentals underpinning Thom's works and leading to surprising findings. 64,65Furthermore, by combining this robust version of BET and other topological concepts, we have recently proposed a model for predicting the activation energy of organic and organometallic reacting systems at 0 K. 42 In the seminal work of Silvi et al. 29 the inversion of ammonia, the dissociation of ethane into two CH 3 radicals, and the N-B bond scission of BH 3 NH 3 molecules were used as prototypical examples to show the BET capabilities.The electronic preparation stages of ethane C-C bond scission were described through an elliptic umbilic and three simultaneous folds.The bond between the two carbons finally cracks following a cusp function.On the other hand, they rationalized the umbrella inversion via two elliptic umbilic polynomials.The topographical description of the dative bond breakage comprises one elliptic umbilic and four fold bifurcations.Despite the cryptic character of the assigning, this methodology is now a standard tool for studying different bonding situations and gaining deeper insights into chemical reactivity in a wide range of research fields.7][68][69][70][71] However, the elliptic umbilic contains three parameters, 30 which must be associated with real objects controlling the reaction.Typically, the distance between two reacting centers is employed as a parameter, 56,57 leaving the physical interpretation of the other coefficients unclear.
This paper aims to re-examine the three reactions discussed by Silvi and co-workers, 29 focusing on the assigning process of bifurcations via monitoring the Hessian determinant at potentially degenerate CPs and the relative distances between them, as recently proposed by Ayarde-Henrı ´quez and collaborators. 56he second goal is to show the need for consistency with the theoretical basis of CT since ignoring them would question the rigor of past and current applications of BET, casting strong doubts on how reactive processes have been rationalized via unfoldings over the last 30 years.

Understanding the pyramidal inversion of ammonia
The process of an atom shifting from its original trivalent or tricoordinate (ground state) configuration through a planar formation and then returning to its initial geometric structure is what we call a pyramidal inversion.This fluxional process is an ideal example of dynamic stereoisomerism, a process where the arrangement of atoms changes without the formation or breakage of bonds. 72,73To understand this process better, we often refer to Linus Pauling's hybridization concept.The original molecule configuration (C 3v ) shifts to a planar configuration (D 3h ) during the transition state (TS).This process can be described by a sequence of hybridization transitions: sp 3sp 2sp 3 (see Scheme 1).
It is well known that the activation barrier (E a ) for the 15 group elements highly depends on the inversion center.The umbrella inversion of ammonia is almost barrierless, demanding only 5.8-5.9kcal mol À1 at room temperature (T = 298 K). 74,75 Similar spectroscopic approaches estimate significantly higher E a for phosphine (PH 3 ) and arsine (AsH 3 ): 17.1 and 32.1 kcal mol À1 , respectively. 76Indeed, a higher value of 31.6 kcal mol À1 has also been reported for PH 3 . 77Predicted activation enthalpies (DH ‡ ), calculated in the gas phase using range-separated corrected global-hybrid DFT functionals 78,79 match closely with the experimental values.Any discrepancies range from 0.17 to 0.78 kcal mol À1 for ammonia, 7.36 to 10.93 kcal mol À1 for phosphine, and 3.71 to 7.82 kcal mol À1 for arsine.These consistent results across various systems suggest that the methodologies used are appropriate for examining the ELF topography along the IRC, yielding consistent results across various systems.Furthermore, these results indicate the analysis performed is independent on the computational method, aligning with recent observations. 56Please refer to the ESI † for additional details.
Using the foundational BET framework, Silvi and collaborators 29 attempted to describe the inversion of ammonia through two elliptic umbilic functions.However, this model does not hold up when we examine the Hessian determinant at potentially involved critical points (CPs) and their relative distance. 56Fig. 1 illustrates the changes in the topographic map before and after the inversion.Nonetheless, it should be stressed that identifying the colliding CPs cannot be performed by a visual inspection but rather by using the mathematical principles supporting CT, which are crucial in selecting the right unfolding featuring the chemical process.Indeed, considering the Poincare ´-Hopf theorem, one can assert that only two CPs coalesce, having such points different indices.Yet, no unfolding can accurately capture this topographic event, as the involved CPs are saddle points. 44,45heme 1 Lewis-like structures depicting the prototypical pyramidal inversion of ammonia.The system commutes between C 3v configurations, (a) and (c) by passing through a D 3h transition state (b).
Fig. 1 The image depicts fragments of the ELF topographic map of the ammonia molecule before (panel a) and after merging two saddle points with differing indices (panel b).Following this merger, there is an appearance of a 2 and s 1d near the remaining saddles, which showcases a fold catastrophe Fig. 2, panel a, is essential to effectively apply Thom's functions 56 since it displays the absolute value of the Hessian determinant, |DH|.This is shown for the index-one saddles (s 1a , s 1b , and s 1c ) and the index-two saddle (s 2a ).It is worth mentioning that this step is not strictly necessary as the collision between saddles leads to topographical bifurcations rather than catastrophic ones.As the phase portrait is about to change, we can observe that the |DH| at s 1c decreases slightly faster than that of the other saddles of index one (s 1a and s 1b ).This is reflective of the impending change.A similar pattern can be seen in the relative distance s 1c -s 2a , indicating that s 1c and s 2a are approaching each other at the highest rate.It is also important to note that as the bifurcation is imminent, the separation symmetry between the s 2a and the three index-one saddles is lost.This transformation is depicted in Fig. 2, panel b.
After the disappearance of saddles of s 2a and s 1c , the remaining saddles, s 1a and s 1b , gradually increase their relative distance, leading to a corresponding increase in their |DH| value.This shift sets the stage for the appearance of a new attractor a 2 and another saddle s 1d near the original saddles, as depicted in Fig. 1, panel c.The emergence of this new attractor can be classified as a fold-type catastrophe, indicating that neither s 1a nor s 1b are degenerate.Within the original BET framework, considering the principles of Thom's work, this observation suggests that there are no signs of an elliptic umbilic yet.As the reaction progresses, on the opposite side of the plane that contains the nitrogen atom and orthogonal to the C 3 axis, the coalescence of a 1 and s 1g forms a wandering point; 7 as seen in Fig. 3, panels a and b.The subsequent loss of this lone pair is best described through a fold catastrophe, as illustrated in Fig. 3, panels d and e. Lastly, the emergence of saddles s 1h and s 2b near the remaining pair (s 1e , s 1f ) marks the final topographic alteration along the IRC.This is shown in Fig. 3, panel c.This analysis reveals that the inversion of the ammonia molecule only triggers flags for the most basic parametric function of Thom's catastrophe theory, contradicting a previous report. 29This significantly simplifies the topological rationalization of this reaction by eliminating the umbilic complexities because this unfolding is not a cuspoid, which means that it has more than one state variable, and more importantly, three (control) parameters are needed to remove its degenerations entirely.The parameter of the fold function typically has a straightforward physical interpretation.One suitable alternative for the chemical process at hand is the angular displacement (m) conveniently defined by the C 3 axis and the N-H bonds. 29See ref. 29 for further details on the critical values of this parameter.In this scenario, smooth changes in m lead to the creation/annihilation of (attractor, saddle) pairs when a critical value m * is reached, ''jumping'' the reacting system between regions (structural stability domains) with any and two CPs, as displayed in Fig. 4. In essence, within this succinct topographical description of the ammonia pyramidal inversion, m controls the changes in sp hybridization.However, it is crucial to underscore that giving a specific interpretation to each parameter in the context of the umbilic catastrophe remains a challenging task.

Dissociation of ethane (breaking of the C-C bond)
The dissociation of ethane involves the breaking of the C-C bond.Under this process, the high molecular symmetry of ethane leads to the production of two identical fragments.Each of these moieties retains an unpaired electron, representing Lewis' concept of shared electron pairs.In some Fig. 2 The diagram showcases two primary aspects of the ammonia umbrella inversion process: panel a represents the alteration in the modulus of the Hessian matrix at the points s 2a , s 1a , s 1b , and s 1c when we traverse the IRC.The inset graph shows the asymptotic behavior of |DH| near the bifurcation, and the dashed red line indicates the transformation occurring in the phase space.Panel b demonstrates the variations in the relative distances between the point pairs s 2a -s 1a , s 2a -s 1b , and s 2a -s 1c as we move along the IRC.This change in distances provides valuable insight into the evolving structure and interaction dynamics during the ethane molecule's dissociation.
This journal is © the Owner Societies 2023 instances, this dissociation may also result in two charged fragments. 80Remarkably, the dissociation of ethane is not a simple process.It demands a significant amount of energy.Experimental data estimates this to be around 89.9 kcal mol À1 at 298 K in the gas phase. 80In our studies using the complete active-space self-consistent field (CASSCF) and its second-order perturbation theory (CASSPT2) methods, we found a range of required energies deviating from 1.56 to 8.15 kcal mol À1 .For detailed information, refer to the ESI.† The first changes in the locality of each carbon take place as the relative separation between s 1a and the three saddles of index two (s 2a , s 2b , and s 2c ) decreases.Since the Poincare ´-Hopf formula 54 dictates both the number and type of CPs variations, only one of the three saddles of index two is involved in the bifurcation, as presented in Fig. 5, panels a and b.We again relied on the |DH| at CPs and their relative distance to determine that such a saddle is s 2a , meaning that this point and s 1a coalesce, and thus, the latter CP does not change its index but rather disappears through a bifurcation that cannot be rationalized using any unfolding; 44,45 see Fig. 5, panels c and d.
Fig. 3 The image features segments of the ELF topographic map of the ammonia molecule before (panel a) and after the coalescence of a 1 and s 1g (panel b).Post this merging, the saddles s 1h and s 2b emerge through a topographical bifurcation, an event that cannot be described in terms of any unfolding (panel c).Panel d visualizes the variation in the modulus of the Hessian matrix at points a 1 , s 1e , s 1f , and s 1g as we move along the IRC.The inner plot depicts the asymptotic behavior of |DH| near the bifurcation, and the dashed red line represents the transformation in the phase space.Panel e shows the change in the relative distances a 1 -s 1e , a 1 -s 1f , and a 1 -s 1g as we progress along the IRC.For ease of understanding, various features of the map are color-coded: attractors in purple, saddles of index one in orange, saddles of index two in yellow, and repellors in green.The map does not display gradient lines for simplicity.Additionally, Lewis-like structures are also provided to give a clearer idea of the associated molecular configurations.
Fig. 4 The image depicts a fragment of the ELF topographic map, showcasing a typical fold flag.This ''jump'' between regions with none and two CPs happens when a critical value of m * is achieved.Attractors are depicted in purple and the saddle points of index one are represented in orange.This image is reproduced from the work of Ayarde-Henrı ´quez et al. 56 As the CC bond further stretches, an index-two saddle (s 2d ) and one repellor (m 4 ) arise in the vicinity of the remaining points, s 2b and s 2c .Therefore, a fold catastrophe describes this topographic event since one of the new CPs is not a saddle, ''jumping'' the reaction system between topographic regions with a difference in the CPs number equal to two.We agree with Silvi and coworkers 29 that the simplest unfolding correctly describes the following three (simultaneous) changes since repellors m 1 , m 2 , and m 3 and saddles s 2b , s 2c , and s 2d collide and annihilate, respectively, as depicted in Fig. 6.
One of the most pivotal and intriguing events during this process is the division of the attractor a 1 into a 2 , a 3 , and s 1b .This event signals the beginning of the CC bond cleavage.For in-depth insights into the critical distances between carbons, we refer the reader to ref. 29.As demonstrated in Fig. 7, the topographic analysis of the ELF gradient offers a precise and intuitive depiction that aligns well with a chemist's expectations.
The absolute value of the Hessian determinant at a 1 diminishes as the distance between the carbon atoms grows.This culminates in the degeneration of this CP and its subsequent division into two new attractors (a 2 and a 3 ) and a saddle of index one (s 1b ).The maintained symmetry in both the |DH| at a 2 and a 3 and the relative distance between these attractors and s 1b along the CC reactive coordinate is of critical importance.These elements corroborate the occurrence of a cusp catastrophe, in line with the classification provided by Silvi and collaborators, 29 see Fig. 7, panels c and d.Although CT offers a straightforward methodology for identifying the correct unfolding to characterize meaningful electron fluxes along reactive coordinates, interpreting the parameter(s) of such a polynomial remains an open area of research, as this exceeds the scope of Thom's theory.By following the logic applied to the inversion of the ammonia molecule, the linear monomial coefficient of the cusp could be linked with the distance between the carbon atoms.Moreover, given that cusp flags have been observed in highly symmetrical reactive systems involving homolytic bond cleavage, 29,60 it seems plausible to associate the other parameter (bias 81 ) with the radical character of the components.The distance between the reacting centers and the biradical nature of the system are not independent variables; therefore, this relationship must be conceptually associated with the one characterizing the cusp parameters, i.e., a semi-cubic parabola.The function describing the unfolding coefficients thus requires us to associate them with physical quantities that are inherently correlated, further enhancing the comprehensiveness of BET applications.From a qualitative perspective, a complete correspondence between the cusp and the ELF geometric shape near the bifurcation requires the use of the dual of this unfolding (i.e., its negative version), as initially reported; 29 see Fig. 8.

Dissociation of ammonia borane (breaking of the B-N dative bond)
The so-called dative bond, less frequently referred to as a covalent bond, 82 a coordinate link, 83 and a semipolar double bond, 84 is an interesting bonding situation in which the lone pair of an atom (e.g., nitrogen) is transferred to the unoccupied orbital of other (e.g., boron).Usually, systems characterized by such a bond are called electron donor-acceptor complexes, 85 Lewis acid-base adducts, 86 or coordination compounds. 87mmonia borane is the textbook example used to illustrate this bonding situation due to its simplicity, 80,88,89 and due to this, it has been investigated through both theoretical 29,90,91 and experimental 92,93 approaches.Regardless that H 3 CCH 3 and H 3 BNH 3 are isoelectronic molecules exhibiting similar structural forms, they differ significantly in their physical and chemical properties, as well as in their dissociation processes.The dissociation of the BN bond in ammonia borane is heterolytic, which means it breaks in such a way that the shared electron pair is entirely taken by one of the atoms.Conversely, the CC bond in ethane breaks homolytically, where each carbon atom takes one electron from the shared pair. 80,88Specifically, the experimental dissociation energy of ammonia borane in the gas phase is about 31.1-37.5 kcal mol À1 at 298 K, 80,94 which is approximately a third lower than the energy required for ethane.In our calculations, using global hybrid functionals, the absolute error ranges from 2.33 to 3.71 kcal mol À1 .In addition, we investigated changes in the molecular graph of phosphine borane along the dissociation pathway.Unfortunately, no experimental dissociation energy values for this molecule were found in the existing literature.As expected, the ELF topological analysis of both systems provides consistent results across different levels of theory.For further details, please refer to the ESI.† The topological explanation for the dissociation of the dative bond between boron and nitrogen in ammonia borane is indeed compelling.Just as with the case of ethane that we discussed earlier, the breaking of the bond can be rationalized and visualized uniquely using topological methods.These methods provide valuable insights into the bond dissociation processes, adding depth to our understanding of the structural changes occurring during these reactions.In the case of the ammonia borane complex, while boron and nitrogen separate, the repellor m and the saddle point s 2d appear within the loop formed by the collection of saddles s 1i and s 2i , i = a, b, and c.Thus, this topographic change corresponds to a fold catastrophe as first identified by Silvi and co-workers; 29 see Fig. 9, panels a and b.The phase space between the reacting centers flattens due to the annihilation of the (m, s 2a ) pair through a fold upon further stretching of the B-N bond (Fig. 9 panels c  and d).Detailed information concerning the B-N critical values can be found in ref. 29.Finally, two pairs of saddles coalesce not simultaneously; see Fig. 9, panels d-f.Regardless that any topographical change resulting from merging/arising saddles cannot be describe via any unfolding, 44,45 the |DH| at all these CPs is provided for illustrating the assigning procedure; see Fig. 9 panel g.Relying on Thom's catastrophe theory, no umbilic fingerprint is detected, as no critical point changes its index.This indicates that a function with one parameter and one state variable is the correct polynomial for characterizing this reaction from a topological standpoint.This is contrary to the report by Silvi and collaborators. 29Additionally, the number of eigenvalues of |DH| tending to zero is a critical variable often overlooked when classifying relevant chemical processes along a reaction pathway via unfoldings.From a pure topological perspective, distinguishing between unfoldings with one and two state variables can only be done considering the Fig. 6 This graphic displays segments of the ELF topographic map of the number of null eigenvalues. 56Hence, we propose that the BH 3 NH 3 dissociation has to be characterized via a cuspoid (i.e., a one-state-variable function) and not in terms of an umbilic (i.e., a two-state-variable function), as only one eigenvalue of the CPs involved in the fold catastrophe tends to zero value.For this reaction, it is appropriate to consider the fold parameter controlling the changes in the ELF topography as the negative of the B-N enlargement since catastrophic bifurcations result from changes in this bond length.

Concluding remarks
This work focuses on the conceptual foundations of identifyingassigning elementary catastrophes within the so-called bonding evolution theory (BET) framework. 29In contrast to the earlier discussion concerning the bonding nature along the ammonia inversion and the dissociation of ethane and ammonia borane molecules, no signatures of the elliptic umbilic unfolding have been detected.Such a finding rest upon rigorous examination of the determinant of the Hessian at critical points and the relative distance between them in each of the reacting systems, a result which is also tested to be independent of the level of theory and substituents.Such a methodological framework leads to rigorously elucidating the nature of each catastrophe by being coherent with the formalities underpinning BET, suggesting that the identified Thom's functions are signatures of the investigated electron rearrangements.From a topological point of view, only two folds are required to rationalize the whole reaction mechanism of the ammonia inversion and ammonia borane dissociation.The CC cleavage in the ethane molecule is correctly characterized using the dual cusp polynomial, in agreement with Silvi and collaborators. 29These results constitute a substantial simplification because the fold is the one-parameter and onestate-variable function, i.e., the simplest polynomial of Thom's theory; in contrast, the elliptic umbilic contains three control parameters and two state variables.Furthermore, it was shown that the physical interpretation of the fold and cusp parameters is somewhat clear, and even though it is still an arbitrary process, we justified the proposed meaning of the latter unfolding coefficients based on their functional relationship, providing a higher degree of completeness to BET applications.Nonetheless, the rationalization of electron reorganizations along a reaction coordinate through higher-order parametric functions continues to be limited by the chemical/physical meaning of such coefficients and their correlations, remaining an open topic for further investigations.
Considering these findings, we here summarize a simple procedure for examining chemical bonding-related processes within the BET methodological framework, which aims to be consequent with the conceptual fundaments of Thom's works, thus unambiguously establishing the link between unfoldings and electron reorganizations, i.e., step (a): the determinant of the Hessian matrix must be computed at all potentially degenerate critical points, i.e., all approaching points along a reactive coordinate, unless they are all saddle points; step (b): the relative distance between potentially degenerate critical points must be computed (unless they are all saddle points) to have additional information, which is particularly useful when the potential energy surface is relatively flat; and step (c): the unfolding assignment must rely on results of steps (a), (b), and considering the maximum number of critical points of each unfolding.The first step is critical since it allows distinguishing between potentially and degenerate points as well as determining the number of eigenvalues tending to zero value, i.e., providing key knowledge concerning how many state variables should have the correct unfolding.It is customary to use the Poincare ´-Hopf theorem to further validate the results of steps (a) and (b).It should be underscored that special care must be considered when applying this methodical yet easy-to-follow receipt to investigate more intricate potential energy surfaces (e.g., almost barrierless rate processes).Reactions involving transition metals constitute typical examples of this critical bonding situation.In these cases, degenerate critical points (CPs) tend to exhibit a relatively low exposure factor (i.e., a metric that quantifies the CP appearance frequency in IRC frames 42 ) due to the flatness of the phase space portrait.To increase this metric, the reaction pathway must be computed using the smallest achievable step size, resulting in more Hessian data points.The ramifications of these findings are considerable.They challenge our present understanding of key chemical processes in terms of unfoldings and cast doubts on both historical and current applications of BET.As such, it is suggested that systems previously examined over the last three decades within the original BET framework should be reevaluated considering these new insights.

Theory and computational details
Within the BET framework, 29 a reaction becomes modeled as a gradient dynamical system, with the electron localization function (ELF) 22 serving as the time-independent ''potential function''. 45ELF acts as a quantum probe for visualizing the Pauli exclusion principle.It is well known that ELF maxima correlate with predictions from the valence-shell electron pair repulsion theory (VSEPR), 95 providing thus a solid connection to Lewis bonding ideas in terms of valence, bonds, core, and lone pairs concepts. 96,97The topographical analysis of ELF is independent of the theoretical level of theory, 56,98,99 and indeed, such a function has also been approximately obtained from electron densities derived from X-ray diffraction data. 100,101The BET approach enables us to identify distinct, neighboring structural stability domains (SSDs), each featuring a unique set of nondegenerate critical points.Every SSD can, in principle, be associated with a given bonding state of the evolving system.The abrupt change from one state to another can be traced back to a bifurcation catastrophe in the form of the potential originated by a small change in parameters on which the function depends.Considering Thom's conceptual foundation for these polynomials (unfoldings), 30 it is necessary to highlight that merely counting changes in the number and type of Morse critical points separating adjacent SSDs is insufficient for determining the associated parametric unfolding.
Let us briefly outline the basis of CT, 30 as comprehensive resources are readily available. 30,44,45,47The gradient nullity condition applied to an R 3 -potential, r r Z(r,m) = 0, where Z stands for the ELF, r represents the Cartesian coordinates, r = r(x, y, z), and r r = (q/qx)i + (q/qy)j + (q/qz)k yields four types of equilibria: attractors, saddles of index one and two, and repellors.The collection of these singular solutions is often referred to as the topographic map, molecular map, or phase-space portrait of the function or, more precisely, of its gradient. 21,102,103The phase space between equilibria, or critical points (CPs) of the function, flattens when the relative distance between two or more equilibria decreases.This change  in the space curvature can be directly computed via the determinant of the Hessian matrix, |DH|.This matrix, H, is a somewhat arbitrary construct for tabulating the second-order derivatives of a function, i.e., H ij = q 2 Z(r, m)/qr i qr j .It is significant to note that i, j = x, y, and z.A CP is termed degenerate if the |DH| evaluated at this point yields a zero value. 44,46,47CT is a robust mathematical program that succeeded in describing all distinct changes in the geometric shape of (generic) parametric families triggered by the creation/annihilation of CPs through seven unfoldings constrained to have four parameters at most.One-state-variable unfoldings are called cuspoids, and their common characteristic is to have only one null eigenvalue.In contrast, parametric functions containing two state variables are referred to as umbilics, which describe changes in the geometric shape of a function resulting from two null eigenvalues.Eigenvalues result from a systematic sequence of algebraic operations applied to a square matrix, such as H, through a process called diagonalization.From a geometric standpoint, this procedure entails matrix rotations aimed at simplifying its form by eliminating mixed spatial-derivative terms.Table 1 lists the adapted 44,47 Thom's unfoldings for suitable usage in bonding situations.The negative sign corresponds to the so-called dual function, whereas including a Morse-type function (y 2 ) in all cuspoids allows their application in 3D phenomena. 44,47ossible combinations of critical points featuring a catastrophe can be straightforwardly found by combining the Poincare ´-Hopf theorem and the maximum number of the unfolding CPs. 64Typical examples of fold in the ELF topography are the following: (attractor, saddle of index one) 2 wp and (repellor, saddle of index 2) 2 wp, where wp stands for a wandering point.For the cusp, one can have: (attractor, attractor, saddle of index one) 2 attractor and (repellor, repellor, saddle of index 2) 2 repellor.Recall that a wp is not a CP, and therefore it describes unbounded trajectories.For further clarification on how to use the information provided in Table 1, let us reconsider the ammonia borane dissociation.We already showed that the fold characterizes the electron rearrangements of this reaction by following the |DH| at potentially degenerate CPs and their relative distance.This information is further verified using Poincare ´-Hopf's relationship, yielding that one repellor and one saddle coalesced.Thus, this abrupt change should be described via a fold or an elliptic umbilic since these are the only functions with two CPs.However, the correct function is the fold since only one eigenvalue tends to zero value.Considering that the control parameter m was conveniently defined as m = d À d * (where d and d * correspond to the initial and bifurcation distances between reacting centers, respectively), the coalesce visualization and description is straightforward.Differentiating the fold with respect to the state variable x yields a quadratic expression with two real solutions, if and only if, m o 0. Indeed, this parameter is negative because the B-N bond is stretching.
All geometry optimizations were performed using the M06-2X 79 and oB97X-D 78 global hybrid functionals, each combined with the 6-31G(d), 6-31G(d,p), 6-31+G(df), 6-31+G(df,p) splitvalence basis sets and using the Gaussian 16 package of programs. 104The pyramidal inversion of ammonia, phosphine, and arsine was conducted following the IRC path with a step size equal to 0.001 amu 1/2 Bohr.In contrast, a relaxed scan with a step of 0.010 Bohr was employed for studying the dissociation of ethane, ammonia borane, and phosphine borane molecules.For the former system, the multiconfigurational character due to the spin multiplicity was assessed through the complete active-space self-consistent field (CASSCF) 105 employing the same number, six and eight, of electrons and virtual orbitals using ORCA 5.0. 106The energy was further corrected via a CAS second-order perturbation theory (CASPT2). 107,108The topological analysis of the ELF and the characterization of its phase space were carried out using the Multiwfn program. 109The images of the topographic map were generated using the VMD software. 110The fulfillment of the Poincare ´-Hopf theorem was verified for all topographic maps along the pathway.The classification of bifurcations was carried out by strictly examining the absolute value of the Hessian determinant evaluated at all potentially degenerate critical points and their relative distance near the degeneration region.

Fig. 5
Fig. 5 The image represents fragments of the ELF topographic map of the ethane molecule before (panel a) and after the merging of s 1a and s 2a (panel b).(Panel c) shows the change in the modulus of the Hessian matrix at s 1a , s 2a , s 2b , and s 2c along the reaction pathway.The inner plot indicates the asymptotic behavior of |DH| near the bifurcation, and the dashed red line marks the change in the phase space.(Panel d) displays the variation of the relative distances s 1a -s 2a , s 1a -s 2b , and s 1a -s 2c along the CC reactive coordinate.The various features of the map are color-coded: attractors are represented in purple, saddles of index one in orange, saddles of index two in yellow, and repellors in green.To simplify the visual, gradient lines are not shown.The dashed lines in the Lewis-like structures represent the bonds that are breaking.

Fig. 9
Fig. 9 ELF topography of ammonia borane before (a) and after s 2d and m arise through a fold (b).The repellor m and the saddle s 2a coalesce through a fold (c) and (d).Sequence of bifurcations that cannot be assigned to any unfolding since only saddles are involved (d)-(f).Change in the modulus of the Hessian matrix at m, s 1i , and s 2i (i = a, b, and c) altogether with the asymptotic behavior of |DH| near the bifurcation (inner plot), dashed red lines indicate changes in the phase space along the BN reactive coordinate (g).Attractors are represented in purple, saddles of index one in orange, saddles of index two in yellow, and repellors in green.Gradient lines are not shown for simplicity.The gradient lines are left out for simplicity, and the dashed lines in the Lewis-like structures indicate the breaking bonds during the reaction.