Towards Modelling Mechanical Shaking Using Potential Energy Surfaces: A Toy Model Analysis

Abstract

In this work, we formalize the effect of mechanical shaking by using various forms of an externally exerted force, which may be constant or may be position-dependent, and we examine the changes in the potential energy surfaces that quantify the chemical reaction. We use a simple toy model to model the potential energy surfaces of a chemical reaction, and we study the effect of a constant or position-dependent externally exerted force for various forms of the force. As we demonstrate, the effect of the force can be quite dramatic on the potential energy surfaces, which acquire new stationary points and new Newton trajectories that are distinct from the original ones that were obtained in the absence of mechanochemical effects. We also introduce a new approach to mechanochemical interactions, using a dynamical systems approach for the Newton trajectories. As we show, the dynamical system attractor properties of the trajectories in the phase space are identical to the stationary points of the potential energy surfaces, but the phase space contains much more information regarding the possible evolution of the chemical reaction—information that is quantified by the existence of unstable or saddle fixed points in the phase space. We also discuss how an experimental method for a suitable symmetric liquid solution substance might formalize the effect of shaking via various forms of external force, even in the form of an extended coordinate-dependent force matrix. This approach may experimentally quantify the Epstein effect of shaking in chemical solutions via mechanochemistry methods.

1. Introduction

Mechanochemical reactions can be categorized into two types of reactions: isotropic and directed ones. In the case of directed mechanochemical reactions, the external forces are exerted in specific directions, usually aligned with the existing molecules’ symmetry axis, and isotropic mechanochemical reactions refer to the collective application of external forces in all directions in the chemical system [1,2]. Although mechanochemistry is not so widespread and popular in the industry, several mechanochemical phenomena exist in nature—for example, in living biological cells, such as in molecular motors, which convert chemical energy into mechanical work used to contract muscles and so on; see [3].

In general, there are various ways to accelerate or enhance the chemical reaction rates—for example, heating, light exposure, the use of catalysts or even the use of mechanochemical techniques. The latter are known to function for some time but are less studied compared to the other techniques, even though scientists such as Michael Faraday studied such techniques. Recently, however, there have been some developments in the theoretical understanding of mechanochemistry. The theoretical treatment of chemical reactions is achieved by studying the Born–Oppenheimer potential energy surfaces. These multidimensional surfaces incorporate the local energy minima of the atomical systems that are involved in the chemical compounds involved in a reaction. The chemical reaction indicates that an existing energy barrier in the potential energy surfaces is crossed, and, if it did not exist, the reactions would take place spontaneously. The minimum energy path between the reactant and the product goes through a saddle point, which is known as the transition state; thus, this causes the multidimensional problem to be reduced to the one-dimensional problem of obtaining these reaction path curves. Mechanochemistry aims to pinpoint the locations and energies of the transition states before and after the exertion of an external mechanical force in the chemical system. Essentially, the mechanical force modifies the activation energy of the system and hence the reaction rate itself is changed. It is thus vital to have knowledge of the potential energy surfaces of a chemical system, and, to this end, the synergy between experiments and theoretical modelling can yield interesting results for future theoretical modelling development; see, for example, ref. [4] for several experiments in this line of research.

In mechanochemistry, it is known that the exertion of an external force along the reaction coordinate(s) may enhance a chemical reaction. In fact, the theoretical description of the above procedure is obtained by studying the potential energy surfaces of molecular systems and how these surfaces are changed or affected by the exertion of external mechanical forces. This theoretical approach is of fundamental importance in mechanochemistry and may reveal information at a molecular level about chemical reactions and how mechanical forces may control chemical processes and reactions. Chemical reactions may be considered in terms of these potential energy surfaces, or as hypersurfaces for higher-dimensional chemical systems, since the potential energy surfaces map the saddle paths and valleys of the molecular potential energy, as the molecules that constitute the molecular system react. The minima, fixed points and several other structural properties of the potential energy surfaces, and how they are changed by the exertion of external force(s), reveal the reaction pathways of a chemical reaction [5].

Recently, the effect of shaking in solutions has been pointed out in Ref. [6], which is referred to as the Epstein effect in this paper. As was pointed out in [6], substances, when exposed to vibrations, acquire new structural characteristics and new physical properties. Clearly, this means that there is some mechanochemical effect that changes the structure of the solutions in some way. The mechanochemical effects of the shaking of aqueous solutions have also been discussed in [7]; see also [8]. In fact, long-lived luminescence was found to occur in deionized water, which was also saturated with atmospheric gases, after mechanical shaking. Moreover, the effects of intense mechanical shaking enhanced the bubbles in the solution, and it was observed that the liquid–gas interface area was significantly enhanced. In addition, it was observed that gas nanobubbles and macromolecules formed in the water after intense mechanical shaking, but the bubbles were absent in the freshly prepared deaerated water. It is known that the mechanical exertion of external forces may change the physicochemical properties of aqueous solutions, and these effects remain for a long time after the exertion of the mechanical stress or shaking. Currently, there is no sufficient explanation of these phenomena, and, in this work, we aim to connect the symmetric isotropic mechanochemically induced reactions with the study of potential energy surfaces in chemical systems. It is possible that these nanobubbles have a general effect in the initial chemical system, by directly altering the potential energy surfaces and the reaction paths and by lowering the energy barriers controlling the chemical reaction, acting as virtual catalysts to some extent [7,8].

In this work, we aim to formalize the effect of shaking by using various forms of externally exerted force, which may be constant or position-dependent. We use a simple toy model to model the potential energy surfaces of a chemical solution, and we study the effect of the constant or position-dependent force for various forms of the force. As we demonstrate, the effect of the force can be dramatic on the potential energy surfaces, which acquire new stationary points and Newton trajectories. We also introduce a new approach to mechanochemical interactions, using a dynamical systems study approach for the Newton trajectories. As we show, the dynamical system attractor properties are identical to the stationary points of the potential energy surfaces, but the phase space contains much more information regarding the evolution of the chemical reaction—information that is quantified by the existence of unstable or saddle fixed points.

2. Potential Energy Surfaces and Mechanochemical Reactions: A Detailed Approach

The potential energy surface is of fundamental importance in theoretical chemistry, along with the reaction paths imagined on them. The reaction paths are essentially one-dimensional descriptions of a chemical reaction embedded in an N-dimensional configuration space of N molecules that take part in the reaction. The potential energy surfaces are effective surfaces that can be deformed due to externally exerted mechanical stress in a specific direction on the molecules, or even in many directions. The potential energy surfaces can be deformed due to this external mechanical stress and thus new minima and saddle points might arise in the new effective potential energy surface. The changes in the minima and the saddle points due to the exertion of the mechanical stress are described by the Newton trajectories in the original potential energy surface. Moreover, the barrier of a reaction is crossed when the norm of the gradient of the potential energy surface along a Newton trajectory is the maximum—a critical point that is known as the barrier breaking point. The importance of the Newtonian trajectories in identifying the reaction path is fundamental since, on a Newtonian curve, the gradient of the potential energy surface points in the same direction for every point on the curve [9].

A mechanochemical reaction is a chemical reaction that is generated or triggered by mechanical stress or energy, according to the IUPAC terminology. The formulation of a mechanochemical reaction is usually materialized by the first-order perturbations of the molecular system due to the exertion of an external force on the system, of the form [9]

$$V_f\left(\overrightarrow{r}\right)=V\left(\overrightarrow{r}\right)-f^T\delta \overrightarrow{r},$$

(1)

where $V_f\left(\overrightarrow{r}\right)$ describes the effective potential of the molecules’ interactions, including the exertion of the perturbation of the external constant force $\overrightarrow{f}$. The linear perturbation is a simplified approach and generalizations may apply [9]. In the case of a constant force, the direction of the force stays the same, and the new stationary points are located in different positions on the effective potential energy surfaces, where it holds that

$${\nabla }_r{V}_f\left(\overrightarrow{r}\right)=0=\overrightarrow{g}-\overrightarrow{f}=0,$$

(2)

where $\overrightarrow{g}={\nabla }_{r}V\left(\overrightarrow{r}\right)$. Essentially, the condition (2) is the condition that the internal forces of the system are equal to the external forces of the system, and it is a vector equation, so it can be analyzed in components. Essentially, the perturbation caused by the external force in the potential energy surfaces is assumed to generate a potential term $f^T\delta \overrightarrow{r}$, so the force acting on the system is essentially generated by the potential term $f^T\delta \overrightarrow{r}$, and we have the condition of equal forces (2), which yields the stationary points of the effective potential surfaces. One, however, may consider the holistic or collective effect of an external mechanical action to have some functional space-dependent form $F\left(\overrightarrow{r}\right)$, and then the effective potential energy surfaces are given by

$$V_f\left(\overrightarrow{r}\right)=V\left(\overrightarrow{r}\right)-F\left(\overrightarrow{r}\right),$$

(3)

and the new stationary points of the effective potential surfaces are given by the condition

$${\nabla }_r{V}_f\left(\overrightarrow{r}\right)=0=\overrightarrow{g}-{\nabla }_{r}F\left(\overrightarrow{r}\right)=0.$$

(4)

The generalized collective action potential may capture and quantify the initial idea of the supramolecular concept of Epstein and reflect extensive shaking, which is the abrupt exertion of an external force in the chemical system [6]. The exact form of $F\left(\overrightarrow{r}\right)$ for a specific chemical system can be obtained only via simulations, where the synergy of simulations and theoretical modelling may be needed.

Note that the conditions (2) and (4) are the conditions of equilibrium internal and external forces; thus, the external force is assumed to be generated by a potential of the form $f^T\delta \overrightarrow{r}$. It is this form that we aim to generalize by using a general potential $F\left(\overrightarrow{r}\right)$ in order to model external shaking, which is the external exertion of an extreme mechanical force.

The second issue that we will address in this work is related to the force-displaced stationary points. For constant forces, the reaction path is obtained by following the force-displaced stationary points, which are found in the solution of the following differential equation:

$$\frac{d\overrightarrow{r}}{dt}=\pm A\left(\overrightarrow{r}\right)g\left(\overrightarrow{r}\right),$$

(5)

where $A\left(\overrightarrow{r}\right)=\det \left(H\right)H^{-1}$. H is the Hessian of the initial potential energy surface $H=\frac{{\partial }^{2}V\left(\overrightarrow{r}\right)}{\partial x_i\partial x_j}$, $i,j=1,\dots ,N$ and N is the dimension of the configuration space of the chemical molecular system under consideration. The curves satisfying the differential Equation (5) are called Newtonian trajectories. Essentially, the differential Equation (5) is an autonomous dynamical system, so the study of the fixed points and trajectories of this dynamical system may reveal the Newtonian trajectories of the potential energy surfaces. Recall that the changes in the minima and the saddle points caused by the external stress are described by the Newtonian trajectories in the original potential energy surface, and this is the information contained in the dynamical system (5). The stationary points for the above dynamical system are apparently the points in which $g\left(\overrightarrow{r}\right)=0$, but there are further fixed points in which $A\left(\overrightarrow{r}\right)g\left(\overrightarrow{r}\right)=0$, and these can easily be found, as we show for some toy model potential energy surfaces. An important feature to notice is that the disturbance of the stationary points of the effective potential energy surfaces is described by the Newtonian trajectories, in the case that only the magnitude of the external force changes but not its direction $I=\frac{\overrightarrow{f}}{\left|f\right|}$, for the linear perturbation models [5]. In the case of a non-constant force-induced potential, this scenario is described when the magnitude of the force ${\nabla }_{r}F\left(\overrightarrow{r}\right)$ has a norm that varies continuously, i.e., $|{\nabla }_{r}F\left(\overrightarrow{r}\right)|$ varies continuously, but the direction of the unit vector $I=\frac{{\nabla }_{r}F\left(\overrightarrow{r}\right)}{|{\nabla }_{r}F\left(\overrightarrow{r}\right)|}$ in the direction of the force ${\nabla }_{r}F\left(\overrightarrow{r}\right)$ remains constant. We shall consider an example of this type in a later section. Then, the Newton trajectories connect the stationary points of the potential energy surface [10]. Regarding the unperturbed molecular system, the reaction paths’ bifurcations are indicated by the existence of valley–ridge inflection points, and these are determined by using the reduced gradient method, which is utilized on the potential energy surfaces of the molecular system. The reduced gradient analysis has similarities with the Newtonian trajectories in the unperturbed system [10].

The minimum energy path is inherently connected to the reaction path of the potential energy surfaces, and this constitutes the main approach to the kinetics of chemical molecular systems. The reaction path is essentially the curve in the configuration space of the molecular system that connects the reactant and the product minimum, passing through the saddle point of a potential energy surface. The saddle points and the minima of the potential energy surface constitute the stationary points of a potential energy surface.

It is possible that reaction path branching may occur during a reaction. The points where this phenomenon occurs are called bifurcation points. These bifurcations of the reaction path may lead to equivalent final states and are related to the valley–ridge inflection points, which, by definition, are those points in the configuration space of the chemical system in which one main curvature of the potential energy surface becomes zero, orthogonally to the gradient. The reduced gradient method yields a family of curves that contain the stationary points of the potential energy surface, and the branching points of these curves are the valley–ridge inflection points of the potential energy surface; thus, the valley–ridge inflection points may be determined by the reduced gradient curves. The reduced gradient curves are obtained by solving the differential equation

$${\nabla }_{r}V\left(\overrightarrow{r}\right)=0,$$

(6)

for all dimensions of the configuration space of the molecular system. Finally, it is worth noting that the reduced gradient technique is strongly related to the Newtonian trajectories approach of Equation (5). Specifically, the reduced gradient approach yields the curves $\frac{\partial V\left(x_i\right)}{\partial x_i}=0$, $i=1,\dots ,N$, which are related to the Newtonian trajectories, because the stationary points of the potential energy surface (which are found by solving $\frac{\partial V\left(x_i\right)}{\partial x_i}=0$, $i=1,\dots ,N$) are limiting points of the Newtonian trajectories obtained by the differential Equation (5) for $g\left(\overrightarrow{r}\right)=0$. Of course, the Newtonian trajectories structure contains more stationary points, which are obtained for $A\left(\overrightarrow{r}\right)g\left(\overrightarrow{r}\right)=0$ with $g\left(\overrightarrow{r}\right)\ne 0$. In fact, these fixed points are the optimal bond breaking points [11].

Now, let us materialize the above formalism and the mechanically induced changes in the potential energy surfaces, by using a test surface for a model, instead of a real molecular potential energy surface. For example, consider a potential energy surface of the form

$$V(x,y)=2y+y^2+(y+0.4x^2)x^2.$$

(7)

Note that the above potential is a toy model and thus different choices might yield different results in the phase space. We choose this model as it is also used in the relevant literature [10], and it is simply a toy model used to demonstrate our ideas. Moreover, the forces are chosen in such a way as to cover most general scenarios. Our aim and motivation is to create sufficient scientific interest in our approach among experimentalists, who will perform experiments to find the potentials of some easy-to-handle substances and also to model and simulate the external forces. Thus, our approach is not an experimentally verified one but serves as a theoretical model in order to highlight the importance of our methods. These can be used by experimentalists to materialize our proposal; see also the discussion in a later section on this.

Now, let us analyze the valley–ridge inflection points and bifurcation points of this toy potential energy surface. Moreover, we shall consider several forms of external force applied in the chemical system and how it deforms the potential energy surfaces and the Newtonian trajectories. Regarding the Newtonian trajectories, we shall adopt an approach in which the dynamical system corresponding to the Newtonian trajectories is studied, using several linearization techniques. We shall reveal the fixed points in the phase space and we shall consider their stability and the connection of the dynamical system’s fixed points with the extremal points of the potential energy surfaces. Let us start with the stationary points and the form of the potential energy surfaces without any external force applied for the potential (7). The stationary points can be found by solving Equation (6), and these are the following, which can be compared to the ones obtained after the exertion of force is also taken into account:

$${\Phi }_0^1:(x,y)=(-1.82574,-2.66667),{\Phi }_0^2:(x,y)=(0,-1),{\Phi }_0^3:(x,y)=(1.82574,-2.66667).$$

(8)

We can also easily find the Newtonian trajectories by drawing the curves corresponding to $\frac{\partial V(x,y)}{\partial x}=0$ $\frac{\partial V(x,y)}{\partial y}=0$. The potential energy surfaces, the stationary points and the fixed points for the potential (7) are presented in Figure 1.

Now, let us proceed to examine in some detail the deformation of the potential energy surfaces and the corresponding stationary points and Newton trajectories in the case that some external force is applied, with the direction of the force being constant, even though the magnitude of the force may change. We start with a simple linear perturbation of the potential (7), in which case the potential has the following form

$$V(x,y)=x^2\left(0.4x^2+y\right)-25x+y^2-23y,$$

(9)

in which case the stationary points are

$${\Phi }_{*}^1:(x,y)=(1.05622,10.9422),$$

(10)

so the deformation of the new center of the interaction is apparent, by comparing the stationary point (10) with the stationary points of the unperturbed potential (8). In this case, the external force has the form

$$\overrightarrow{f}=-25\overrightarrow{i}-25\overrightarrow{j},$$

(11)

and it has a constant direction. Moreover, in Figure 2 we plot the potential energy surfaces, the stationary points and the Newton trajectories of the deformed potential (9).

Thus, we can see how the exertion of an external force deforms the potential energy surfaces and the stationary points. Thus, the chemical reaction is deformed due to a mechanochemical external effect. Let us quote several more examples to see the effect of various forms of the external force. At the end, we shall also consider forces that have no constant direction in space. Consider the following perturbation of the potential (7), in which case the potential has the following form

$$V(x,y)=x^2\left(0.4x^2+y\right)-0.25x+y^2+1.75y,$$

(12)

which is a rescaled form of the previously used potential. In this case, the stationary points are

$${\Phi }_{*}^1:(x,y)=(-1.63133,-2.20563),{\Phi }_{*}^2:(x,y)=(-0.143878,-0.88535),{\Phi }_{*}^3:(x,y)=(1.77521,-2.45069),$$

(13)

so the deformation of the new center of the interaction is apparent, by comparing the stationary point (13) with the stationary points of the unperturbed potential (8). In this case, the external force has the form

$$\overrightarrow{f}=-0.25\overrightarrow{i}-0.25\overrightarrow{j},$$

(14)

and it has a constant direction. It is apparent that, in this case, the simple rescaling of the external force has a quite distinct effect, different from the one in which the force is given by (11). In this case, we have three distinct stationary points; thus, simple rescaling has a profound effect on the potential energy surfaces. Moreover, in Figure 3, we plot the corresponding potential energy surfaces, the stationary points and the Newton trajectories for this case.

Thus, we can see how the exertion of an external force deforms the potential energy surfaces and the stationary points. It is worth noting that the simple deformation of the external force, i.e., rescaling, results in three different stationary points in the deformed potential energy surface. Consider now a subcase of the potential (9) perturbation of the potential (7), in which case the potential has the following form

$$V(x,y)=x^2\left(0.4x^2+y\right)-25x+y^2+2y,$$

(15)

in which case the stationary points are

$${\Phi }_{*}^1:(x,y)=(3.78648,-8.1687),$$

(16)

so the deformation of the new center of the interaction is apparent, by comparing the stationary point (16) with the stationary points of the unperturbed potential (8). In this case, the external force has the form

$$\overrightarrow{f}=-25\overrightarrow{i},$$

(17)

and it has a constant direction. Moreover, in Figure 4, we plot the potential energy surfaces, the stationary points and the Newton trajectories of the deformed potential (15).

Thus, a simple subcase of the initially exerted external force leads to quite different potential energy surfaces. In addition, let us simply change the direction of the force, to see the effects of this action on the potential energy surfaces. In this case, the potential has the following form

$$V(x,y)=x^2\left(0.4x^2+y\right)+25x+y^2+27y,$$

(18)

in which case the stationary points are

$${\Phi }_{*}^1:(x,y)=(-7.13046,-38.9217),{\Phi }_{*}^2:(x,y)=(0.944659,-13.9462),{\Phi }_{*}^3:(x,y)=(6.1858,-32.6321),$$

(19)

so the deformation of the new center of the interaction is apparent, by comparing the stationary points (19) with the stationary points of the unperturbed potential (8). In this case, the external force has the form

$$\overrightarrow{f}=25\overrightarrow{i}+25\overrightarrow{j},$$

(20)

and it has a constant direction opposite the one appearing in Equation (11). Moreover, in Figure 5, we plot the potential energy surfaces, the stationary points and the Newton trajectories of the deformed potential (18).

Let us present another interesting case before proceeding to the more generalized forms of the potential. In this case, we increase the magnitude of the force in the y-direction. In this case, the potential has the following form

$$V(x,y)=x^2\left(0.4x^2+y\right)-25x+y^2-998y,$$

(21)

in which case the stationary point is

$${\Phi }_{*}^1:(x,y)=(0.025,498.999),$$

(22)

so the deformation of the new center of the interaction is apparent, by comparing the stationary point (22) with the stationary points of the unperturbed potential (8). In this case, the external force has the form

$$\overrightarrow{f}=-25\overrightarrow{i}-1000\overrightarrow{j},$$

(23)

and it has a constant direction. Moreover, in Figure 6, we plot the potential energy surfaces, the stationary points and the Newton trajectories of the deformed potential (21).

It apparent that the exertion of an external force that is stronger in the y-direction has a significant impact on the potential energy surfaces’ structure. Now, let us adopt an entirely different approach, in which the external force has no linear perturbation effect on the potential but is more general. Consider the case in which the potential has the following form

$$V(x,y)=x^2\left(0.4x^2+y\right)-x^4-2.5x^2-y^4-1.5y^2+2y,$$

(24)

in which case the stationary point is

$${\Phi }_{*}^1:(x,y)=(0,0.5),$$

(25)

so the deformation of the new center of the interaction is apparent, by comparing the stationary point (25) with the stationary points of the unperturbed potential (8). In this case, the external force has the non-constant magnitude form

$$\overrightarrow{f}=(-4x^3-5x)\overrightarrow{i}+(-4y^3-5y)\overrightarrow{j}.$$

(26)

Moreover, in Figure 7, we plot the potential energy surfaces, the stationary points and the generalized Newton trajectories of the deformed potential (24).

Now, consider a slightly different form of the potential (24), i.e., the following

$$V(x,y)=x^2\left(0.4x^2+y\right)-x^4+2.5x^2-y^4-1.5y^2+2y,$$

(27)

in which case the stationary points are greater in number than in the previous case, including the previous stationary point, and these are

$${\Phi }_{*}^1:(x,y)=(0,0.5),{\Phi }_{*}^2:(x,y)=(-1.66578,0.829794),{\Phi }_{*}^3:(x,y)=(1.66578,0.829794),$$

(28)

so the deformation of the new center of the interaction is apparent in this case too, by comparing the stationary point (28) with the stationary points of the unperturbed potential (8). In addition, a slight change in the potential produces more stationary points, as is apparent. In this case, the external force has the form

$$\overrightarrow{f}=(5x-4x^3)\overrightarrow{i}+(-4y^3-5y)\overrightarrow{j},$$

(29)

so it is the same in the y-direction as in the previous case, but different in the x-direction. Moreover, in Figure 8, we plot the potential energy surfaces, the stationary points and the Newton trajectories of the deformed potential (27).

Now, consider another slightly different form from the potential (24), in which case the potential has the following form

$$V(x,y)=x^2\left(0.4x^2+y\right)-x^4+2.5x^2-y^4+3.5y^2+2y,$$

(30)

in which case the stationary points are quite numerous and are the following

$$\begin{array}{cc}\hfill & {\Phi }_{*}^1:(x,y)=(0,-1.1462),{\Phi }_{*}^2:(x,y)=(-1.14614,-0.923631),{\Phi }_{*}^3:(x,y)=(1.14614,-0.923631),\hfill \\ \hfill & {\Phi }_{*}^4:(x,y)=(-1.22935,-0.68645),{\Phi }_{*}^5:(x,y)=(1.22935,-0.68645),{\Phi }_{*}^6:(x,y)=(0,-0.301352),\hfill \\ \hfill & {\Phi }_{*}^7:(x,y)=(0,1.44755),{\Phi }_{*}^8:(x,y)=(-1.85069,1.61008),{\Phi }_{*}^9:(x,y)=(1.85069,1.61008),\hfill \end{array}$$

(31)

so the deformation of the new center of the interaction is apparent, by comparing the stationary point (31) with the stationary points of the unperturbed potential (8). A simple comparison with the fixed points (25) also reveals great differences, caused by the slight deformation of the potential for the two cases. In this case, the external force has the form

$$\overrightarrow{f}=(5x-4x^3)\overrightarrow{i}+(5y-4y^3)\overrightarrow{j}.$$

(32)

Moreover, in Figure 9, we plot the potential energy surfaces, the stationary points and the generalized Newton trajectories of the deformed potential (30).

In conclusion, all examples considered point out one aspect: the mechanochemical disturbance of a molecular system can significantly affect the chemical reaction, since it significantly affects the potential energy surfaces, the stationary points and the Newton trajectories connecting the stationary points. Now, the challenge is to model, in an optimal way, the mechanochemical effect of shaking in order to quantify the Epstein effect of shaking. This could be a challenge, because the form of the force would probably be space-dependent and could potentially have the form studied in the above three examples. It is also notable that, in the three cases studied, the number of stationary points is significantly altered by a slight sign change in the external space-dependent force. This is quite interesting from a phenomenological point of view, since the actual mechanochemical external effect might actually change its sign along various directions. It is a challenge to model such effects as the Epstein effect; in this section, we saw qualitatively, in terms of some simple models, the strength of the effects of mechanical forces on the potential energy surfaces of a chemical system. Thus, the synergy between experiments and theory can yield some information about the actual mechanochemical effect of shaking in some molecular solutions. We shall discuss this issue further in the Discussion section.

3. Newton Trajectories, Study of the Dynamical System and Proven Equivalence of the Two Approaches

In this section, we shall use a novel approach to study the Newton trajectories in the potential energy surfaces, and we shall demonstrate that the stationary points of the potential energy surfaces are essentially stable fixed points of the dynamical system that describes the Newton trajectories, namely the dynamical system of Equation (5). Thus, using a standard phase space analysis for the dynamical system (5), we shall provide some new insights into the structure of the potential energy surfaces. Regarding the dynamical system study, we shall explore the phase space using standard linearization techniques, which we will review shortly. The dynamical system (5) can take the following form for a two-dimensional space

$$\begin{array}{cc}\hfill & \frac{dx}{dt}=f_1(x,y),\hfill \\ \hfill & \frac{dy}{dt}=f_2(x,y),\hfill \end{array}$$

(33)

where $A\left(\overrightarrow{r}\right)g\left(\overrightarrow{r}\right)={(f_1(x,y),f_2(x,y))}^T$. Thus, we shall calculate the resulting vector $A\left(\overrightarrow{r}\right)g\left(\overrightarrow{r}\right)$ for each case of combined potential and applied force considered in the previous section, and we shall find the fixed points of the resulting dynamical system. Then, we shall investigate the stability of the fixed points and the behavior of the trajectories in the phase space; thus, we shall reveal the actual Newton trajectories in the phase space. Moreover, we shall show the strong relation of the stationary points of the potential energy surface and the corresponding stable fixed points of the phase space of the dynamical system of Equation (5), or equivalently that of Equation (33).

To start with, let us review in brief the linearization techniques that we shall consider. With the stability theory of the phase space of dynamical systems, one in fact reveals the stability of the fixed points and finds the attractors of the trajectories of the dynamical system. Specifically, the actual stability of the trajectories studies the behavior of these trajectories, subject to small perturbations in the initial conditions. In addition, one useful piece of information for the convergence of the solutions to the fixed point values is offered by the time variable, but we shall not consider this issue in detail here because it is out of context and would heavily rely on the specific model used. In all cases, the question posed is how the phase space orbits and trajectories behave in relation to the fixed points and whether they converge towards the stable fixed point, if any. Another question is whether some trajectories pass through unstable fixed points before they reach the stable fixed points, which we shall call attractors. Moreover, it is important to examine the above questions for various choices of initial conditions of the phase space variables. The fixed points of a dynamical system provide a concrete criterion that determines the structural stability of a given dynamical system. The only way to reliably reveal the stability of a fixed point of a dynamical system is to directly apply the Hartman–Grobman theorem, by appropriately linearizing the autonomous non-linear dynamical system. Then, by calculating the eigenvalues of the linearization matrix, which are real negative numbers at the fixed points or even complex numbers with negative real parts, the fixed point is deemed stable, and thus it is an attractor in the phase space. If the eigenvalues are positive and negative, the fixed point is a saddle; finally, if the eigenvalues are positive, the fixed point is unstable. We shall consider only hyperbolic fixed points, which means that the eigenvalues of these points have only non-zero values. As we will show, this is the case for all combinations of potential energy surfaces and applied forces. Now, let us briefly review the Hartman–Grobman theorem, which is a linearization theorem and in fact determines the stability and the phase space structures when only hyperbolic fixed points are considered. Let $\Phi \left(t\right)$ $ϵ$$R^n$ satisfy the following flow equation

$$\frac{\mathrm{d}\Phi }{\mathrm{d}t}=g(\Phi \left(t\right)),$$

(34)

where $g(\Phi (t\left)\right)$ is a continuous map that is also locally Lipschitz $g:R^n\to R^n$. Let ${\varphi }_{*}$ be all the fixed points of the dynamical system (34), and also the corresponding Jacobian matrix $\mathcal{J}\left(g\right)$ is defined to be equal to

$$\mathcal{J}=\sum _i\sum _j\left[\frac{\partial {\mathrm{f}}_{\mathrm{i}}}{\partial x_j}\right].$$

(35)

The Jacobian needs to be evaluated at all hyperbolic fixed points, and the eigenvalues $e_i$ must satisfy $\mathrm{Re}\left(e_i\right)\ne 0$. Let the spectrum of the eigenvalues of a matrix A be denoted as $\sigma \left(A\right)$; hence, a hyperbolic fixed point must satisfy $\mathrm{Re}\left(\sigma \left(\right(J\left)\right)\right)\ne 0$. The Hartman–Grobman theorem guarantees the existence of the homeomorphism $\mathcal{F}:U\to R^n$, with U being an open neighborhood of ${\varphi }_{*}$, which satisfies $\mathcal{F}\left({\varphi }_{*}\right)$. The homeomorphism generates a non-trivial flow $\frac{\mathrm{d}h\left(u\right)}{\mathrm{d}t}$, which satisfies

$$\frac{\mathrm{d}h\left(u\right)}{\mathrm{d}t}=\mathcal{J}h\left(u\right),$$

(36)

and this is a topologically conjugate flow of the flow in Equation (34). Note that the Hartman–Grobman theorem is applied only to autonomous dynamical systems, and this is the case studied in this paper. Thus, in view of the Hartman–Grobman theorem, the dynamical system of Equation (34) can be written in the following way

$$\frac{\mathrm{d}\Phi }{\mathrm{d}t}=\mathcal{J}\left(g\right)(\Phi ){|}_{\Phi ={\varphi }_{*}}(\Phi -{\varphi }_{*})+\mathcal{S}({\varphi }_{*},t),$$

(37)

with $\mathcal{S}(\varphi ,t)$ being a smooth map $[0,\infty )\times R^n$. Therefore, if the Jacobian matrix $\mathcal{J}\left(g\right)$ satisfies the relation $\mathcal{Re}\left(\sigma \left(\mathcal{J}\right(g\left)\right)\right)<0$ and, in addition, if

$$\underset{\Phi \to {\varphi }_{*}}{lim}\frac{\left|\mathcal{S}\right(\varphi ,t\left)\right|}{|\Phi -{\varphi }_{*}|}\to 0,$$

(38)

the fixed point ${\varphi }_{*}$ of the dynamical flow $\frac{\mathrm{d}\Phi }{\mathrm{d}t}=\mathcal{J}\left(g\right)(\Phi ){|}_{\Phi ={\varphi }_{*}}(\Phi -{\varphi }_{*})$ is also a fixed point of the flow in Equation (37), and it is asymptotically stable. Now, let us apply the Hartman–Grobman theorem for the dynamical system (33) for each of the combination of forces and potentials studied in the previous section. In all cases, the potential energy surface before the application of the external force is that of Equation (9).

For the purposes of this work, the choices for the dynamical system study are again based on a toy model approach. We choose two cases, a force-free case and a case in which a force is applied, for simplicity. Our aim is to show that there exist structures in the phase space and that important information can be gained by this type of study. It is notable that an actual model derived from experiments is needed, so our approach is purely demonstrational, aiming to motivate studies by experimentalists in this line of research. Let us consider first the force-free case, so the functions $f_1(x,y)$ and $f_2(x,y)$ appearing in Equation (33) are

$$\begin{array}{cc}\hfill & f_1=1.2x^3-4x,\hfill \\ \hfill & f_2=x^2(7.6y+9.6)+1.6x^4+y(4y+4),\hfill \end{array}$$

(39)

and thus the matrix $\mathcal{J}={\sum }_i{\sum }_j\left[\frac{\partial {\mathrm{f}}_{\mathrm{i}}}{\partial x_j}\right]$ in this case equals

$$\mathcal{J}=\left(\begin{array}{cc}3.6x^2-4& 0\\ 6.4x^3+2x(7.6y+9.6)& 7.6x^2+8y+4\end{array}\right).$$

(40)

The fixed points of the dynamical system (33), for the case at hand, are

$$\begin{array}{cc}\hfill & {\varphi }_{*}^1=(-1.82574,-4.66667),{\varphi }_{*}^2=(-1.82574,-2.66667),{\varphi }_{*}^3=(0,-1),\hfill \\ \hfill & {\varphi }_{*}^4=(0,0),{\varphi }_{*}^5=(1.82574,-4.66667),{\varphi }_{*}^6=(1.82574,-2.66667).\hfill \end{array}$$

(41)

Note that the fixed point ${\varphi }_{*}^3$ is a stationary point on the potential energy surface, appearing in Equation (8), denoted as ${\Phi }_0^2$. Now, regarding the stability of the fixed points, the eigenvalues of the linearization matrix evaluated for each fixed point are

$$\begin{array}{cc}\hfill & {\varphi }_{*}^1\to (8,-8),\hfill \\ \hfill & {\varphi }_{*}^2\to (8,8),\hfill \\ \hfill & {\varphi }_{*}^3\to (-4,-4),\hfill \\ \hfill & {\varphi }_{*}^4\to (-4,4),\hfill \\ \hfill & {\varphi }_{*}^5\to (8,-8),\hfill \\ \hfill & {\varphi }_{*}^6\to (8,8),\hfill \\ \hfill \end{array}$$

(42)

thus, ${\varphi }_{*}^2$ and ${\varphi }_{*}^6$ are unstable, ${\varphi }_{*}^3$ is stable and the rest are saddles. We can easily evaluate numerically the trajectories in the phase space for various initial conditions. The result of our numerical analysis for the case at hand is depicted in Figure 10. The stable fixed points are denoted with magenta dots, the saddles with red dots and the unstable fixed points with blue dots.

As can be seen in Figure 10, all trajectories are attracted to the stable fixed point of the dynamical system (39), which is identical to the stationary point of the potential energy surface.

Now, let us consider the case in which the force is given by Equation (14), so the functions $f_1(x,y)$ and $f_2(x,y)$ appearing in Equation (33) are

$$\begin{array}{cc}\hfill & f_1=\frac{\left(0.214286x^3-0.625x-0.0892857\right)\left(5.6x^2+4y+0\right)}{1x^2+0.714286y+0},\hfill \\ \hfill & f_2=\left(5.6x^2+4y\right)\left(\frac{\left(x^2+2y+1.75\right)\left(4.8x^2+2y\right)}{5.6x^2+4y}-\frac{0.571429x\left(1x^3+1.25xy-0.15625\right)}{1x^2+0.714286y}\right).\hfill \end{array}$$

(43)

The fixed points of the dynamical system (33), for the case at hand, are

$$\begin{array}{cc}\hfill & {\varphi }_{*}^1=(-1.63133,-3.72575),{\varphi }_{*}^2=(-1.63133,-2.20563),{\varphi }_{*}^3=(-0.143878,-0.88535),\hfill \\ \hfill & {\varphi }_{*}^4=(-0.143878,-0.0289814),{\varphi }_{*}^5=(1.77521,-4.41193),{\varphi }_{*}^6=(1.77521,-2.45069).\hfill \end{array}$$

(44)

Note that the fixed point ${\varphi }_{*}^3$ is a stationary point on the potential energy surface, appearing in Equation (12), denoted as ${\Phi }_{*}^2$, and ${\varphi }_{*}^2$ corresponds to ${\Phi }_{*}^1$ and ${\varphi }_{*}^6$ to ${\Phi }_{*}^3$. Now, regarding the stability of the fixed points, the eigenvalues of the linearization matrix evaluated for each fixed point are

$$\begin{array}{cc}\hfill & {\varphi }_{*}^1\to (5.58051,-5.58051),\hfill \\ \hfill & {\varphi }_{*}^2\to (6.58051,5.58051),\hfill \\ \hfill & {\varphi }_{*}^3\to (-3.92548,-2.92548),\hfill \\ \hfill & {\varphi }_{*}^4\to (-3.92548,3.92548),\hfill \\ \hfill & {\varphi }_{*}^5\to (7.34497,-7.34497),\hfill \\ \hfill & {\varphi }_{*}^6\to (8.34497,7.34497),\hfill \\ \hfill \end{array}$$

(45)

thus, ${\varphi }_{*}^2$ and ${\varphi }_{*}^6$ are unstable, ${\varphi }_{*}^3$ is stable and the rest are saddles. We can easily evaluate numerically the trajectories in the phase space for various initial conditions. The result of our numerical analysis for the case at hand is depicted in Figure 11. The stable fixed points are denoted with magenta dots, the saddles with red dots and the unstable fixed points with blue dots.

As can be seen in Figure 11, all trajectories are attracted to the stable fixed point of the dynamical system (43), which is identical to the stationary point of the potential energy surface. As is apparent, the phase space structure is richer in structures in comparison to the potential energy surface, and thus might provide interesting new insights regarding the direction of the chemical reaction and how it is affected by the external force.

4. Proposal for an Experimental Method

Now, the mechanochemistry approach for shaking in aqueous chemical solutions might be used to theoretically explain the Epstein effect. In fact, the synergy between suitably chosen experiments for suitably chosen substances and theoretical modelling in terms of mechanochemistry may provide a force matrix that can model the Epstein effect of extensive shaking in those substances under consideration. Indeed, consider that the effect of shaking is not the simple exertion of an external force in the symmetry axis of a substance, but that it might be the exertion of multiple forces in various directions of the substance, and, in some cases, this effect can be position-dependent. Thus, if a suitable substance with a high degree of symmetry in its molecules, and at a specific low temperature, is shaken considerably in a specific way, several molecular experimental techniques may reveal the ways in which shaking acts on the molecules of the substance. From this procedure, the overall force matrix may be revealed for this specific substance. Thus, with the initial potential energy surface, one can use the force matrix obtained and calculate in detail the new stationary points of the potential energy surfaces with the inclusion of the force matrix. This procedure is feasible, although technically demanding, but we discuss this issue here given the importance of the synergy between theory and experiments in the above research context. Note that, from the cases studied, it is evident that these are parts of a toy model. We have revealed, for very simple cases, the rich structure of the phase space and the potential energy surfaces; however, for experiments, one needs a realistic potential and realistic forms of the forces. Then, one can obtain true and realistic information about several substances and their behavior after external forces are applied in the molecules in an aqueous solution. Thus, experiments are needed to provide the realistic potential and the forms of the forces.

5. Conclusions

There are many studies to be done at multiple levels, both theoretical and experimental, regarding the effects of mechanical shaking in liquid solutions and the effects of shaking on the reactions and physical properties of the solutions. Apart from modelling mechanical shaking using a supramolecular matrix approach for the isotropic exertion of external forces, one must also obtain deeper knowledge of how liquids function, and this is not a simple task. For example, it is known that liquids support solid-like modes of vibration with wavelengths that extend to the shortest distance, which is comparable to interatomic separation [12]. The way that collective excitations propagate in liquids is not yet fully understood, and these in fact propagate to the smallest wavelengths and do not decay or even become damped. This is a remarkable property of liquids that is also found in solids, and, recently, it has been pointed out that liquids and solids are of a dual nature [12]. Liquids behave similarly to gases as they also flow, but they also behave as solids as the intermolecular forces are strong and displacements are large [12]. It is surprising that the specific heat problem for liquids is not addressed even in standard textbooks. Thus, regarding liquid solutions, future work can model the propagation of the external vibrations caused by shaking, which is essentially the isotropic exertion of external forces in the system. The question of how these collective exertions of external forces propagate in a liquid solution may be answered by studying the effective potential energy surfaces, but this is not an easy task. A synergy between experiments and theory is required in order to first quantitatively model, in a collective way, the effect of shaking by using a supramolecular force matrix [6], probably by modelling the exertion of a force in multiple directions, and to study the effects on the molecules of the liquid solution. Future research in the above areas would certainly provide us with several surprising results.