A catastrophe theory-based model for optimal control of chemical reactions by means of oriented electric fields

The effect of oriented external electric fields (OEEF) on chemical reactivity has been studied theoretically and computationally in the last decades. A central goal in this research area is to predict the orientation and the smallest amplitude electric field that renders a barrierless chemical process with the smallest possible strength. Recently, a model to find the optimal electric field has been proposed and described (Bofill JM et al., J. Chem. Theory Comput. 18:935, 2022). We here proof that this model is based on catastrophe and optimum control theories. Based on both theories a technical treatment of the model is given and applied to a two-dimensional generic example that provides insight into its nature and capability. Finally, the model is applied to determine the optimal OEEF for the trans-to-cis isomerization of a [3]cumulene derivative.

In the context of OEEF-driven catalysis, the direction in which an OEEF is applied is certainly one of the key variables that need to be considered. In a recent work, we introduced a theoretical framework that allows establishing the optimal direction in which an external electric field should be applied to accelerate a given reaction in the most efficient form [48]. In particular, the optimal OEEF obtained in our model is defined as the smallest amplitude field that renders a barrierless chemical process. Our theoretical framework is based on a generalization of the Newton trajectories (NT) model [49] and on the construction of an effective potential energy surface (PES) built from an unperturbed PES of the molecular system under study and the energy change due to the action of an OEEF on the dipole moment of the system. Given a chemical process and an external electric field with a fixed direction that is assumed to be adequate to promote the reaction, our ansatz provides theoretical tools to evaluate how the structures of the reactant and transition state configurations are distorted under the effect of the field. One of the main concepts of the proposed model is the so-called force displacement stationary points (FDSP) curve. There exists a different FDSP curve for each direction of an OEEF. Each point of a given FDSP curve corresponds to a stationary point of the effective potential generated by a particular intensity of the applied OEEF. For a well-behaved system, the OEEF will distort the reactants in such a way that their nuclear configuration becomes closer to the transition state (TS) configuration. By the same token, the TS configuration will become closer to the reactants. As the OEEF amplitude increases, reactants and TS become more distorted and their configurations resemble to a larger extent. There exists an OEEF amplitude at which both configurations coalesce. This specific FDSP point is called bond-breaking-point (BBP). Since each FDSP has its own BBP and each FDSP corresponds to a different OEEF direction, it follows that for each OEEF direction there exists a different BBP. For a given set of OEEFs, each one well suited to promote the reaction of interest, the optimal OEEF (i.e., the one that renders the reaction barrierless with the smallest possible amplitude) corresponds to the FDSP curve that passes through the socalled optimal BBP (oBBP). At this point, the gradient of the original or unperturbed PES is an eigenvector of zero eigenvalue of the Hessian matrix of the effective PES.
The main purpose of this article is to analyze key mathematical concepts introduced in our previously proposed model [48] for defining optimal OEEFs. In particular, we will show that the FDSP curve and BBP points are concepts rooted in catastrophe and optimal control theories. The article is structured as follows: Sect. 2 is devoted to discuss the mathematical nature and structure of the model and for proving the optimality of the OEEF derived in [48]. Section 3 shows that a detailed numerical analysis of a generic two-dimensional model system allows one to deeply understand the full nature of the proposed model. In Sect. 4 we apply the model to a real chemical system, namely, the trans-to-cis [3]cumulene isomerization. In Sect. 5 the main conclusions of the present study are provided.
2 Nature, structure, and treatment of the model

Basic nature and structure of the model
In this subsection, we will first summarize some of the key concepts and mathematical expressions of the model presented and discussed in [48]. We will then prove that this theory falls within the catastrophe theory. The action of an OEEF, = E n , with modulus E and direction n , modifies the original PES of a molecular system. It is assumed that the modified PES is given by the following ansatz [48] where V( ) is the original PES, P n ( , E) is the perturbation energy, ( ) is the dipole moment vector, is the vector of all M atoms' Cartesian coordinates, and T n = (e x , e y , e z ) is the three-dimensional normalized field direction vector. We recall that the dipole vector is a three-component vector, where each component is dependent (in a nonlinear manner) on the -vector, T ( ) = (d x ( ), d y ( ), d z ( )) . We assume that V( ) and d i ( ) are "well-behaved" and ( ) = ∇ V( ) . We note that − n ( )E = ∇ P n ( , E) , is the gradient of P n ( , E) with respect to . Let us consider a given normalized constant vector n , then the set of points ( , E) satisfying Eq. (2) is called force-displaced stationary point (FDSP) curve. At each point of this curve, ( (t), E(t)) , where t is the parameter that characterizes the curve, an effective PES, V n ( , E(t)) , is generated that satisfies stationary condition Eq. (2). Thus, for each normalized n -vector the corresponding FDSP curve generates a sequence of V n ( , E) effective potentials, and each effective potential has its own fixed E, being the coordinates the set of variables. Thus, the set of coordinates plays the rule of variables, whereas the intensity E the rule of parameter. The FDSP curve is a special NT curve where the amplitude of field E changes at each point to satisfy Eq. (2). At each point of a FDSP curve the parameter E coincides with the quotient A point outside the FDSP curve the value of the quotient q( ) does not coincide with E. Thus, at a stationary point of the effective PES, in particular at any point of a FDSP curve, the parameter E coincides with the quotient value between the norm of ( ) and the norm of n ( ) . It is interesting to note that the normalized vector n plays the rule of control axis. It then follows that the associated E is the control parameter. Thus, we emphasize that the effective energy potential, V n ( , E) , can be seen as a function of and the parameter E. When the control variable E has a fixed value the system settles into a structure where the variables stationarize (locally) the function V n ( , E) . In particular, the current point (t) of the FDSP curve corresponding to the control axis n is a stationary point of V n ( , E) since this point satisfies Eq. (2) for the current control parameter E(t). As the control variable changes, namely from E(t) to E(t + Δt) , a local stationary point can disappear and the variables jump suddenly to a different structure. In particular the current stationary point changes from (t) to a new point, (t + Δt) , satisfying Eq. (2) for the new E(t + Δt) . In a similar manner, for a given normalized n -vector the manifold of points satisfying det[ n ( , E)] = 0 varies as the control parameter E varies, where n ( , E) = ∇ ∇ T V n ( , E) is the Hessian matrix of V n ( , E) . For the E parameter associated with a given normalized n -vector the corresponding FDSP point is also a point belonging to the manifold det[ n ( , E)] = 0 , this point is a degenerate stationary point of the current V n ( , E) . The structure of V n ( , E) around this point is the object of analysis of catastrophe theory and has a shoulder form. In the present context, we call this point bond-breaking-point (BBP) ( BBP , E BBP ) associated with the control axis n . The explanation given above is represented schematically in Fig. 1. Finally, if we assume that a reaction path, (s) , is a parameterization of a PES through the so-called reaction coordinate s, then the effect of an OEEF on a chemical reaction system according to this model is schematically explained in Fig. 2.

Mathematical treatment
As shown in [48], the integration of the FDSP curve is an important ingredient of the model. Here, we will explain in detail some mathematical tricks regarding this integration that were not presented in [48].
Basically, the differential equation that gives us the tangent of the FDSP curve is obtained by applying the directional derivative to the FDSP curve condition, Eq. (2), Equation (5) can be written in a more compact way as where the explicit form of this Hessian matrix reads being ( ) = ∇ ∇ T V( ) the Hessian matrix of the original PES, V( ) , and is the Hessian matrix of P n ( , E) with respect to . We note that the Hessian matrix, n ( , E) , is symmetric since ( ) and −E⟨ ( ) n ⟩ are symmetric matrices. Equation (6) tells us how the planes tangent to the iso-contours, V( ) = and P n ( , E) = , with parallel normals, ( ) and −E n ( ) , respectively, are transported through the FDSP curve associated with the field n .
The practical problem to integrate the FDSP curve under a given normalized constant n -vector implies solving Eq. (6). For this purpose we first consider a predictor step. If det[ n ( , E)] ≠ 0 , the tangent, d ∕dt , is computed through the expression, If we are in a point of the FDSP curve, then E is computed from Eq. (4). We note that the value of dE/dt is only a scaling factor on the resulting predictor direction. Normally, a fixed step length through the normalized tangent criterion is used, namely, (d ∕dt) T (d ∕dt) = 1 . From this normalization we compute first, det( n ( , E i )) = 0 , that is to say, the determinant of the second derivatives with respect to of the function V n ( , E i ) vanishes. The point of the FDSP curve that belongs to this manifold is a degenerate stationary point labeled as ( BBP , E BBP ) . Each control axis, n , has a different ( BBP , E BBP ) point. We note that a point not belonging to the FDSP, the corresponding quotient q( ) defined in Eq. (4), does not coincide with the E parameter We discuss now the initial value problem to integrate a FDSP curve. Once the original PES, V( ) , and the dipole function, ( ) , are given, we need as initial values a fixed normalized field vector, n , and a point in configuration space, 0 . With this information, the electric field amplitude E 0 is computed through Eq. (4). Applying Eq. (6) to this initial point, ( T 0 , E 0 ) , a new point is found as explained in the previous paragraph. We proceed in this way until we find a point where the absolute value of E is below a given threshold.
It is usually convenient to consider that the initial point 0 is associated with the reactants and that the initial amplitude of external field action is null. That is, ( T 0 , E 0 ) = ( T min , 0) , where min is a minimum of the original PES, V( ) , associated with the reactants, i.e., ( min ) = and det[ ( min )] > 0 . Since ( min ) = then q( min ) = 0 and E 0 = 0 due to Eq. (4). Thus, given the normalized n -vector, the Hessian matrix ( ) at = min is computed. Since E 0 = 0 then ( min ) = n ( min , 0) . With these elements at hand and the vector, n ( min ) , we can apply Eqs. (10) and (9) and a new point is obtained. With this new point a new value E is computed through Eq. (4). We proceed in this way until the BBP of this FDSP curve is located. The corresponding E BBP is the maximum amplitude of the given external field, n , to be applied on this system. The integration process is continued until a point where E is below to a given threshold in absolute value is reached.

Optimal nature of the model
Having presented the OEEF model as well as a way to solve it computationally, a remaining question that needs to be addressed is: is there an OEEF with normalized direction n where ̇ = d ∕dt . Equation (12) tells us that at the BBP the tangent ̇ is an eigenvector of the matrix, n ( , E) , resulting from the sum of the two matrices, ( ) and −E⟨ ( ) n ⟩ . Thus the expectation value of these two matrices, ̇ T ( )̇ ∕(̇ Ṫ ) , and ̇ T [ ( ) − E⟨ ( ) n ⟩]̇ ∕(̇ Ṫ ) , should be the same but opposite sign at the BBP, ( T BBP , E BBP ) . We conclude that the functions V( ) and P n ( , E) at the BBP coincide in the first and second derivatives with respect to , but with the opposite sign. Thus the BPP point is the contact point of order two of the functions V( ) and P n ( , E) for E = E BBP . Two functions are said to be contact order two in a point if the first and second derivatives of these functions coincide in this point.
At any point of a FDSP curve, such curve crosses the isocontours of V( ) = and P n ( , E) = , through the tangent planes of these iso-contours with parallel normals, ( ) and −E n ( ) , respectively, forming with the tangent ̇ an angle. At a BBP, both functions, V( ) and P n ( , E) , are contact of order two, which means that ̇ T n ( BBP , E BBP )̇ ∕(̇ Ṫ ) = 0 (see the discussion of the previous Eq. (12)). Now we consider the case that the tangent vector ̇ forms a zero angle with the parallel normals ( ) and −E n ( ) in the BBP. We note that in this case the relation , being this equality a particular case occurring when Eq. (13) is satisfied at the BBP. Finally, if we differentiate at the BBP the quotient q( ) that coincides with the amplitude of the electric field E, see Eq. (4), with respect to and we impose condition Eq. (13), we obtain, In the derivation of Eqs. (17), (4), (13) and (7) have been used. Thus the control axis n that generates a FDSP curve such that the corresponding BBP satisfies Eq. (13) enables the transformation of Eq. (12) into Eqs. (15) and (16) in addition to Eq. (17). This control axis is labeled optimal OEEF, n , and the BBP associated is the optimal BBP (oBBP). To analyze more deeply Eq. (17) we apply the resolution of identity in this expression obtaining, Thus, according to Eq. (18), the projection of ∇ q( ) in the space spanned by ̇ -vector results in ̇ T ∇ q( ) = dq( )∕dt = 0 coinciding with dE∕dt = 0 , which is a condition satisfied by any BBP, see Eq. (11a). In addition the projection of ∇ q( ) in the outer space orthogonal to ̇ -vector should be the zero vector, . These two stationarity conditions constitute the optimal conditions that an optimal OEEF should satisfy, see Fig. 3. Algorithms to locate the oBBP and the optimal OEEF are reported and described in Refs. [48,50,51]. These results are a generalization of the corollary described in Ref. [52] for mechanochemistry. In this case, the second derivative with respect to variables of corresponding perturbed energy term of Eq. (1) is equal zero since it depends linearly with respect to variables; thus, the right-hand side term of Eq. (15) is equal zero and the FDSP curve at the oBBP coincides with a gradient extremal point [53] of V( ) function. (18)

A Generic Example
In this section, we will show, through a generic two-dimensional example, the features of the model explained in previous Sect. 2 and more specifically in Sects. 2.1 and 2.3. The original PES, V( ) , used to illustrate the model was taken from Ref. [52]. Notice that in the present case = (x, y) .

The expression of this function is
Concerning the function P n ( , E) = P n (x, y, E) , we have selected the nonlinear dipole vector field reported in [48] such that this function takes the final form where The first and second derivatives of V(x, y) and P n (x, y, E) with respect to x and y, given in Eqs. (19) and (20), respectively, are reported in Appendix 1.
We will now analyze two different FDSPs: one passing through the oBBP and the other one passing through a nonoptimal BBP. Both curves start at (x, y) = (1.0, 0.0) , which is a minimum energy configuration of V(x, y), assumed to be the reactant species. The curve passing through the oBBP corresponds to an optimal external field or control axis, whereas the curve passing through a BBP corresponds to a nonoptimal external field. Since in both cases we start at the minimum point of V(x, y), then as discussed in Sect. 2.2, E 0 = 0 and the tangent vector is computed by the expres- n (x, y) are evaluated using Eqs. (28) and (30), respectively, once the field (e x , e y ) is given. The term dE∕dt| x=1,y=0,E=0 is determined through Eq. (10) taking into account that n (x, y, E)| x=1,y=0,E=0 = (x, y)| x=1,y=0 . Proceeding as explained in Sect. 2.2, the FDSP curves are computed until their corresponding BBP (the point where det( n (x, y, E)) = 0 and dE∕dt = 0 ) is reached. The FDSP curve generated by the optimal field or control axis reaches the oBBP. All units of the quantities appearing in this section are expressed in arbitrary units.
The normalized nonoptimal field used for the present discussion is n = (−0.373, −0.928) T , whereas the optimal one is n = (−0.549, −0.836) T . The optimal normalized field n was found solving Eq. (14b) by the procedure briefly discussed in Ref. [48]. We recall that the set of conditions given in Eq. (14) is the unique one to be satisfied by the field or control axis to be optimal. These conditions differentiate the optimal field from the nonoptimal fields as discussed in Sect. 2.3.
In Fig. 4 we represent the effective potential V n (x, y, E) for different values of the E-parameter from E = E 0 = 0 until E = E BBP generated during the integration of the FDSP curve corresponding to the nonoptimal OEEF with normalized direction n = (−0.373, −0.928) T . At E = E BBP = 3.97 we found the maximum value of the E-parameter being located at the point (1.470, −0.017) , which is the BBP of this FDSP curve generated by the normalized control axis n = (−0.373, −0.928) T . On the other hand, in Fig. 5 we represent the effective potential V n (x, y, E) for different values of the E-parameter from E = E 0 = 0 until E = E oBBP Fig. 3 Optimal OEEF control of a FDSP curve. A scheme illustrating our framework to solve the problem of FDSP curve generated by an OEEF associated with chemical reaction by optimizing the E-parameter. The spatio-external electric field profile is defined by { , n } , where is the set of variables and n the control axis. The chemical transformation along the FDSP involves the evolution of the reactant at the position 0 and field, n , to a position BBP with the same field n in a finite E value. The final E value is E BBP or E oBBP generated during the integration of the FDSP curve corresponding to the optimal OEEF with normalized direction n = (−0.549, −0.836) T . At E = E oBBP = 3.77 we found the maximum value of the E-parameter being located at the point (1.476, 0.072), which corresponds to the oBBP, since the FDSP curve at this point satisfies the conditions given in Eq. (14), numerically demonstrated below. Notice that E oBBP = 3.77 < E BBP = 3.97 , which can be seen as a numerical proof of Eq. (18) and schematically represented in Fig. 3.
At the oBBP, located at (1.476, 0.072), one of the expectation values of the Hessian matrix associated with V(x, y) function coincides with one of the expectation values of the Hessian matrix associated with P n (x, y, 3.77) function with n = (−0.549, −0.836) T . Taking into account the results discussed in Sect. 2.3, the oBBP, as well as any other BBP, is the point where the functions V(x, y) and P n (x, y, E) given in Eqs. (19) and (20), respectively, are of contact order two. The tangent vector of the nonoptimal FSDP at (x, y) BBP = (1.470, −0.017) is an eigenvector of Hessian matrix, n (x, y, E) , with null eigenvalue, being this Hessian function of Hessian matrices of the functions V(x, y) and P n (x, y, E) with respect to x, y, see Eqs. (11b) and (7). The tangent of this FDSP curve at the BBP does not coincide with the gradient of V(x, y) and P n (x, y, E) with respect to x, y. On the other hand, for the optimal FDSP curve the tangent at (x, y) oBBP = (1.476, 0.072) coincides with the gradients of V(x, y) and P n (x, y, E) with respect to x, y. Notice that in both cases, the gradients of V(x, y) and P n (x, y, E) are collinear since Eq. (2) should be satisfied at any point of any FDSP curve. Now we proof numerically the results described in Eqs. (12), (15), and (16) for these nonoptimal and optimal FDSP curves. For the nonoptimal FDSP, the expectation values of (x, y) and [∇ (x,y) ∇ T (x,y) P n (x, y, E)| E=E BBP ] matrices To analyze numerically these results, we first show in Fig. 6 the function V(x, y) given in Eq. (19) and the function P n (x, y, E) given in Eqs. (20) and (21) for the optimal n (x, y) . The normalized gradients (x, y) and n (x, y) and the normalized tangent of this FDSP curve at the BBP point form an angle of 6.5 • . This angle is very small which means that the optimal OEEF is close to the present field and also the BBP is close to oBBP. Equation (15) is not satisfied since we are in a BBP. In the present case the vector, differs numerically from the vector, These two vectors are depicted in yellow in Fig. 6a n (x, y) = −0.0089 −0.3071 . Fig. 6 Representations of original PES, V(x, y), given in Eq. (19) and perturbation, P n (x, y, E) , given in Eqs. (20) and (21) for the nonoptimal, figure b, and optimal, figure c, OEEFs. The E-parameter is taken the corresponding E max value ▸ Thus, neither (x, y) nor n (x, y) are eigenvectors of (x, y) and [∇ (x,y) ∇ T (x,y) P n (x, y, E oBBP )] , respectively, at oBBP point. Finally, the norm of the vector of Eq. (14b) is ≃ 10 −5 indicating that this point is an oBBP. Finally, taking Eq. (18) and the numerical results of Eqs. (23) and (26) we see that ∇ (x,y) E = (0, 0) T , proving that the normalized axis-control, n = (−0.549, −0.836) T , is the optimal OEEF. With these results we have proved numerically the analytical expressions derived and discussed in Sect. 2.

Application of the model to a real chemical example: trans-to-cis isomerization of a [3]cumulene derivative
We now show the application of the model to a real system. The aim of this example is to show how the optimal external electric field annihilates the energy barrier of a chemical process. We have chosen the trans-to-cis isomerization of a [3]cumulene derivative because recent experiments have shown that the OEEFs accelerate isomerizations in this type of systems [54] (see Fig. 7 for the molecular structure of the n (x, y) . system under consideration). We have described the isomerization process considering a single variable, namely, the dihedral angle characterized by the C 1 -C 2 -C 3 -C 4 atoms, see Fig. 7. The potential energy profile with respect to this dihedral angle was computed by means of constrained geometry optimizations in which the mentioned dihedral was fixed and the remaining degrees of freedom were optimized. We scanned the dihedral from 0 • to 180 • with a step of 10 • . By convention the s-cis isomer is described by a 0 • dihedral angle, while the s-trans is described by a 180 • dihedral angle. The calculations reported in this section were carried out within the framework of ab initio single-determinant spinunrestricted Hartree-Fock (UHF) method [55,56] with the 6-31G* basis set [57] employing the ORCA 5.0.3 suite of programs [58][59][60][61][62]. The computational approach relied on the broken-symmetry formalism; in this case, we calculated the set of orbitals associated to the lowest triplet state and used them as a guess to perform the calculations of the singlet state.
It is now possible to locate the oBBP on the potential energy profile and to compute the optimal OEEF. Figure 8 shows the potential energy profile and the location of the oBBP. The oBBP is the point at which the dihedral angle is equal to 140 • . The energy difference between the transition state and the oBBP is equal to 2.37 kcal/mol; the energy difference between the oBBP and the s-trans minimum is equal to 1.47 kcal/mol. The molecular structure of the oBBP and the optimal OEEF is reported in Fig. 9. In Fig. 10 we report the potential energy profile calculated under the effect of the optimal OEEF. It nicely shows how  the energy barrier is completely removed. This real example shows that conformational control by OEEF provides a plausible mechanism for electrostatic molecular motors.

Conclusions
With the present study we have proven analytically and numerically that the model presented in [48] is based on the catastrophes and optimal control theories. We have proved that the optimality of the field is due to a type of control axis that leads to the optimal BBP (oBBP). The oBBP point has a unique mathematical property: its E parameter or amplitude of electric field coincides with the quotient q( ) , which is stationary in the direction of the tangent of the optimal FDSP curve and also stationary in the outer subspace. We think that catastrophes and optimal control theories provide an important tool to analyze and understand the complexity of chemical process and to seek possible control of chemical mechanisms. Fig. 10 Potential energy curve of the s-cis/s-trans isomerization. In blue the original potential energy profile that it is calculated without the effect of the external electric field, = E n = . In orange the potential energy profile calculated with the effect of the optimal OEEF. The values are relative to the energy minimum without the external electric field. The red dot represents the oBBP point located at dihedral angle 140 •