A Second-order Melnikov function

# Persistence of instanton connections in chemical reactions with time dependent rates

## Abstract

The evolution of a system of chemical reactions can be studied, in the eikonal approximation, by means of a Hamiltonian dynamical system. The fixed points of this dynamical system represent the different states in which the chemical system can be found, and the connections among them represent instantons or optimal paths linking these states. We study the relation between the phase portrait of the Hamiltonian system representing a set of chemical reactions with constant rates and the corresponding system when these rates vary in time. We show that the topology of the phase space is robust for small time-dependent perturbations in concrete examples and state general results when possible. This robustness allows us to apply some of the conclusions on the qualitative behavior of the autonomous system to the time-dependent situation.

###### pacs:
05.40.-a, 05.45.-a, 64.60.My, 82.20.-w

## I Introduction

Understanding the dynamics of chemical kinetics is a question of fundamental character in nonequilibrium statistical mechanics, apart from being of broad importance in applications to other sciences. Indeed, many models in chemistry (1), biochemistry (2), ecology (3), and biology (4) use stoichiometric relations as theoretical first principles describing some phenomenon. The simplest description of the time evolution of a set of reacting species is probably given by the mean-field equations, an dimensional dynamical system representing the concentrations, or total number of molecules of the reacting species. One of the advantages of this approach is that it allows us to use the powerful machinery of dynamical systems theory (5), raising the possibility of identifying stationary states with fixed points, periodic behavior with limit cycles, etc.

Of course, as with every theory in physics, the mean-field approximation has a range of validity. As it ignores fluctuations, its description of the chemical system might be accurate for short times; however, long-time dynamics will be affected, dramatically in some cases, by rare events. It is possible to study exactly the system evolution as a time-continuous Markov process, the probability distribution of which is given by the solution of an adequate master equation (6); (7). The full analytical solution of the master equation is, in most situations, a formidable problem, and it yields, usually, much more information than that needed in applications. The development of approximate theories is thus clearly justified.

One of the most common approximations is the Fokker-Planck equation, obtained from the master equation by means of a Kramers-Moyal or van Kampen size expansion (6); (7). This theory assumes that the number of reagents is very large, so we can consider the implicit stochasticity of the process as small Gaussian fluctuations around the mean-field behavior. Definitely, rare events do not belong to this category. If one wants to deal with fluctuations that are comparable with, or even greater than, the mean value, any approximation scheme must not rely on assumptions such as an asymptotically large value of the number of reagents. It is possible to construct such a theory, for instance, by taking advantage of the relation among chemical kinetics and quantum mechanics (8); (9). Once the master equation is formulated as a quantum problem, it is possible to develop eikonal (10); (11) or WKB approximations (12), or matched asymptotic expansions of the spectral formulation (13), able to tackle rare events. In this work we will concentrate on the eikonal approximation introduced in (11). One of the advantages of this approach is that it reduces the problem again to a dynamical system of dimensions for reacting species. The additional degrees of freedom are the conjugate ”momenta” corresponding to each of the concentrations and represent a measure of the size of the fluctuations. Due to the Hamiltonian symmetry of this system, it can be effectively reduced to a dynamical system on some Riemannian manifold. So, if there is only one chemical species, the case in which we will concentrate here, we have a reduced number of dynamical scenarios. In particular, only fixed points will be of physical relevance, and they will denote the possible stationary states in which the system can be found. Contrary to the mean-field situation, rare events can drive the system for one stationary state to another, a fact that is reflected by the existence of connections between the fixed points of the dynamical system. These connections are not the unique, but the optimal way in which a system evolves from one state to the other (10), and, as they are reminiscent of instantons in quantum mechanics (14), we will denote them as instanton connections. Its importance is huge: the web of connections, having the fixed points at their intersections, encodes the qualitative behavior of the chemical system (11), and its topology serves as a principle for the classification of nonequilibrium phase transitions in reaction-diffusion models (15).

All these works have been focused on the dynamics of chemical kinetics happening at constant rates. While this assumption is reasonable in many cases, there are situations in which we should go beyond it and consider the explicit time variation of the reaction rates. Periodically illuminated chemical reactions (16) or seasonal variation in population dynamics (17) serve as examples of nontrivial behavior generated by a temporal forcing. Not much attention seems to have been paid to eikonal approximations of time-dependent chemical kinetics, and when they are considered, nonautonomous perturbations are usually treated within the quasistationary approximation (18); (19). It is our goal to extend the existent approaches and consider arbitrary frequency, albeit small, perturbations. The concrete problem under study is how the phase portrait of the eikonal dynamical system is modified when time-dependent perturbations enter into play: do the instanton connections survive or do they disappear changing the qualitative behavior of the system? Giving a totally general answer to this question is a difficult task, but we will see that these connections seem to be very robust and persistent when they are subject to a small periodic forcing. This will be shown with the help of particular models, for which rigorous results (for the eikonal approximation) are easily proven, and then we will extend them to the general setting when possible.

## Ii Instanton Persistence

### ii.1 Branching and annihilation

For the shake of clarity, we will illustrate the problem of instanton persistence with a particular reaction, branching and annihilation of identical particles, but we will state general results for arbitrary reaction sets (20). Let us start considering a single species of identical particles , which annihilate in pairs and undergo binary branching

 A+A\lx@stackrelλ(t)⟶∅,A\lx@stackrelσ(t)⟶A+A. (1)

The master equation describing the probability distribution of having reagents at time reads

 dPn(t)dt=σ(t)[(n−1)Pn−1(t)−nPn(t)]+λ(t)2[(n+2)(n+1)Pn+2(t)−n(n−1)Pn(t)]. (2)

It can be represented by a partial differential equation (PDE) by introducing the generating function

 G(p,t)≡∞∑n=0pnPn(t), (3)

which provides us with the time-dependent Hamiltonian

 H(p,q,t)=σ(t)(p−1)pq+λ(t)2(1−p2)q2 (4)

and the imaginary time Schrödinger equation

 ∂tG=−H(p,−∂p,t)G. (5)

The eikonal approximation proposes the reduction of the problem to Hamilton equations (11), which in this case become the nonautonomous dynamical system (note, however, that we could reduce the problem to one time-dependent reaction rate by means of the time reparametrization )

 ˙p = −∂H∂q=σ(t)(1−p)p+λ(t)(p2−1)q, (6) ˙q = ∂H∂p=σ(t)(2p−1)q−λ(t)pq2. (7)

Suppose for a moment that both and are time independent. In this case the system exhibits three lines with zero energy: the invariant lines , , and

 q=2σp/λ1+p. (8)

Furthermore, we know that these three lines connect the fixed points , , and in the plane, see Fig.(1) (or Fig.(3) in (11)). Now we can try to understand what happens when we let the reaction rates vary in time. It is easy to see that the lines and remain invariant zero-energy lines of the system; however, the explicit time dependence of the Hamiltonian prevents the conservation of the energy and, in general, forces the disappearance of the third zero-energy line. We can generalize this fact for an arbitrary reaction Hamiltonian. Given any set of reactions with time-dependent rates, we know that the Hamiltonian necessarily fulfills

 H(p=1,q,t)=0 (9)

due to the conservation of probability (11); (15). So this means that for any set of (time-dependent) reaction rules, the line is an invariant, zero-energy line of the dynamical system. Since this line describes the mean-field dynamics of the system, we will call it the mean-field line (11); (19). Some systems possess an absorbing state when they contain zero particles; this happens if all reactions of the type

 ∅\lx@stackrelαn(t)⟶nA (10)

are absent in the dynamics. In this case, the Hamiltonian must obey the condition (11); (15)

 H(p,q=0,t)=0. (11)

So we can claim that any (nonautonomous) system without particle production from the vacuum has the invariant, zero-energy line . And thus, when it is present, we will call it the absorbing-state line. These properties can be clearly seen from the general form of the Hamiltonian term representing the reaction

 mA\lx@stackrelβmn(t)⟶nA; (12)

it is (15)

 Hmn=βmn(t)m!(pn−pm)qm. (13)

In the autonomous situation, the set of zero-energy lines determines the physics of the problem: the fixed points obtained when these lines cross represent the possible states in which the system can be found and the connections among them the possible transition paths. Together with the mean-field line (the global minimum of the action) and the absorbing-state line one finds other lines with zero energy: the instanton lines (11); (15); (19). We know that both the mean-field line and the absorbing-state line persist if we let the reaction rates vary in time, but however, it is not so obvious to see what happens with the instanton lines. These are defined in terms of zero energy if the system is not explicitly dependent on time, but when this is not the case, the definition loses its meaning since energy is no longer conserved. Due to this fact and because the physical role of the instanton lines is to be optimal paths between different states, what we would like to know at this point is if both fixed points and connections among them survive after the nonautonomous forcing is switched on. We can be sure that hyperbolic fixed points persist to a small periodic time-dependent forcing, after possible relocation of their position, but do the instanton connections persist?

To address this question we will use the method developed by Melnikov (21) (see also (5)). Consider an autonomous two-dimensional dynamical system

 ˙x=f(x), (14)

, , with two hyperbolic fixed points and linked by a heteroclinic connection, which is parametrized by the system solution for initial time . Assuming that the motion goes from to , this connection is formed by a branch of the unstable manifold of , which totally overlaps with a branch of the stable manifold of . Let us now consider the perturbed version of this problem,

 ˙x=f(x)+ϵg(x,t), (15)

where the perturbation is time periodic with period , amplitude small enough and sufficiently regular. The dynamics of this nonautonomous system is given by the associated Poincaré map , which maps every initial condition point with the corresponding value of the solution after one period has elapsed (5). The hyperbolic fixed points of the unperturbed system, and , are hyperbolic fixed points of the Poincaré map . Since the map is a perturbation of , these points have a continuation, and , as hyperbolic fixed points of , and their invariant manifolds vary continuously with respect to . The perturbed system will not, in general, maintain the coincidence between the branches of the unstable and stable manifolds of and , respectively: now these branches might intersect, preserving the existence of the heteroclinic connection (see Fig.(2)) or might not, destroying it, like in Figs.(3) and (4). The distance between the stable and unstable branches is given by , with

 M(t0)=∫∞−∞f(xh(t−t0))∧g(xh(t−t0),t)dt, (16)

where

 f∧g=∣∣∣f1f2g1g2∣∣∣ (17)

denotes the wedge product of vectors and . Equation (16) defines the so-called Melnikov function, which yields the first-order approximation in of the distance between the stable and unstable manifolds measured along a direction that is perpendicular to the unperturbed connection at the point . A change of sign of means that there exists some such that , implying the existence of a solution of Eq. (15) defining a heteroclinic connection among the two hyperbolic fixed points and of the Poincaré map corresponding to Eq. (15), say,

 limt→−∞xϵh(t)=xϵa,limt→∞xϵh(t)=xϵb. (18)

Let us now consider the branching and annihilating system, Eqs. (6) and (7), with constant reaction rates. The instanton connection (the separatrix linking to ) is parametrized by the system solution

 xh(t−t0)=(ph(t−t0),qh(t−t0))=(11+eσ(t−t0),2σ/λ2+eσ(t−t0)). (19)

If we assume that the perturbation of the reaction rates is given by

 σ(t) = σ+ϵσ1(t), (20) λ(t) = λ+ϵλ1(t), (21)

and the new rates are periodic and regular enough, we obtain the system

 ˙p = σ(1−p)p+λ(p2−1)q+ϵ[σ1(t)(1−p)p+λ1(t)(p2−1)q], (22) ˙q = σ(2p−1)q−λpq2+ϵ[σ1(t)(2p−1)q−λ1(t)pq2], (23)

which is of type (15). The associated Poincaré map still has the origin as hyperbolic fixed point. Another hyperbolic fixed point of this map, corresponding to when , is , where is obtained through the solution of Eq. (23) for . This is nothing but a Bernoulli differential equation, which can be straightforwardly integrated to get

 q(t)=q(0)exp[∫t0σ(τ)dτ]1+q(0)∫t0exp[∫s0σ(τ)dτ]λ(s)ds. (24)

We just need to apply the condition to this solution to find

 qϵ=exp[∫T0σ(τ)dτ]−1∫T0exp[∫s0σ(τ)dτ]λ(s)ds. (25)

To determine if the unstable manifold of intersects the stable manifold of we substitute the corresponding values of and into the Melnikov function (16),

 M(t0)=∫∞−∞[λσ1(t)−σλ1(t)][(1−p)(1−p−p2)q2](xh(t−t0))dt= ∫∞−∞h(t)ϕ(t−t0)dt=∫∞−∞h(s/σ+t0)~ϕ(s)ds, (26)

after the change of variables and where

 h(t) = λσ1(t)−σλ1(t), (27) ϕ(t) = [(1−p)(1−p−p2)q2](xh(t)), (28) ~ϕ(s) = −4σ2λ2e2s[1+2sh(s)](1+es)3(2+es)2, (29)

the symbol standing for the hyperbolic sine. It is easy to see that has zero mean,

 ∫∞−∞~ϕ(s)ds=−4σ2λ2∫∞−∞e2s[1+2sh(s)](1+es)3(2+es)2ds=2σ2λ2[es(1+es)2(2+es)]∞−∞=0, (30)

a fact that will prove its usefulness later on. Since is periodic and continuous, we can expand it in Fourier series,

 h(t)=∞∑n=−∞ane2πint/T, (31)

and substitute it into the Melnikov function to obtain

 M(t0)=−4σ2λ2∞∑n=−∞ane2πint0/T∫∞−∞e2s[1+2sh(s)](1+es)3(2+es)2e2πins/(σT)ds, (32)

which is the Fourier decomposition of the Melnikov function. We can check the validity of the Fourier series (32) after contrasting that the norm

 ∫∞−∞∣∣∣e2s[1+2sh(s)](1+es)3(2+es)2∣∣∣ds=5√5−112 (33)

is finite. This representation allows us to calculate the mean value

 1T∫T0M(t0)dt0=a0∫∞−∞~ϕ(s)ds=0, (34)

where we have used

 1T∫T0e2πint0/Tdt0=δn0, (35)

and denotes the Kronecker delta. So we know that the Melnikov function is periodic, continuous (from its very definition in Eq. (16) it only depends on , , and ), and with zero mean, implying that it either crosses zero or vanish identically. In the first case, we can claim that the point is connected to the origin, in the second, we cannot, due to the possible presence of small terms . However, we can see that the second case only happens when is constant. Indeed, integrating by parts and changing variables we obtain

 ∫∞−∞e2s[1+2sh(s)](1+es)3(2+es)2e2πins/(σT)ds=iπnσT∫∞−∞es(1+es)2(2+es)e2πins/(σT)ds= iπnσT∫∞0x2πin/(σT)(1+x)2(2+x)dx. (36)

This last integral can be obtained by means of contour integration in the complex plane. Noticing that for a complex variable

 z2πin/(σT)=e2πinln(z)/(σT), (37)

we can use the keyhole contour choosing the logarithm branch cut as the positive real axis, and by employing the residue theorem we obtain

 ∫∞0x2πin/(σT)(1+x)2(2+x)dx=2πi1−e−4π2n/(σT){Res[f(z),−1]+Res[f(z),−2]}, (38)

where

 f(z)=z2πin/(σT)(1+z)2(2+z). (39)

We can now compute the residues

 Res[f(z),−2] = exp[−2π2nσT+2πinσTln(2)], (40) Res[f(z),−1] = ddzz2πin/(σT)(2+z)∣∣∣−1=−(1+2πinσT)exp(−2π2nσT), (41)

to conclude

 ∫∞−∞e2s[1+2sh(s)](1+es)3(2+es)2e2πins/(σT)ds= 2π2n/(σT)1−exp[−4π2n/(σT)]exp(−2π2nσT)(1+2πinσT−exp[2πinσTln(2)]), (42)

and one can see that it only vanishes if or , for and positive finite real numbers. Hence, the Melnikov function is identically zero if and only if for all in Eq. (31) or, what is the same, if and only if is constant.

We still have to analyze what is the meaning of the condition

 λσ1(t)−σλ1(t)=c, (43)

for some constant . If , then the system can be reduced to the unperturbed one by means of a reparametrization of the time variable:

 u=∫t0[1+σ1(τ)σ]dτ. (44)

In this case the geometry of the phase space is exactly preserved; only the parametrization of the system solution along the separatrices might change. If and are chosen constant, then for any value of , the phase-space topology is the same, as a small retuning of the Hamiltonian parameters cannot modify it. In these two particular cases we know that the Melnikov function is identically zero because the connection is still given by the complete superposition of a branch of the stable manifold of one of the fixed points and a branch of the unstable manifold of the other. In the general case we know that a perturbation fulfilling is a combination of these two, translation and reparameterization,

 σ(t) = (σ+ϵσ1)[1+ϵϕ(t)], (45) λ(t) = (λ+ϵλ1)[1+ϵϕ(t)], (46)

for some function , up to terms . This means that, in order to know if the distance between the invariant manifolds of the fixed points is identically zero, we would have to use the extension of the Melnikov method to higher orders (23).

### ii.2 General setting

It would be highly desirable to extend two properties of this example to general systems, say, the possibility of expressing the Melnikov function as a Fourier series, and to show that it has zero mean. We will start with the second claim assuming that the first one is true, and we will check its validity afterwards. The most general reaction Hamiltonian can be built adding the generic terms (13),

 H=∑m≠nβmn(t)m!(pn−pm)qm,m,n=0,1,2,⋯, (47)

and it allows us to express the dynamical system as

 ˙p = ∑m≠nβmn(t)(m−1)!(pm−pn)qm−1,m,n=0,1,2,⋯, (48) ˙q = ∑m≠nβmn(t)m!(npn−1−mpm−1)qm,m,n=0,1,2,⋯, (49)

where all the rates are supposed to have the same period (the case of perturbations with different periods will be discussed below). Assuming the form of the perturbation allows us to split the Hamiltonian into two terms . So we can write the integrand of the Melnikov function as the set of Poisson brackets

 f(xh(t−t0))∧g(xh(t−t0),t)={H1(xh(t−t0),t),H0(xh(t−t0))} (50)

and obtain

 M(t0)=∫∞−∞{H1(xh(t−t0),t),H0(xh(t−t0))}dt= ∫∞−∞ddtH1(xh(t−t0),t)dt−∫∞−∞∂∂tH1(xh(t−t0),t)dt= −∫∞−∞∂∂tH1(xh(t−t0),t)dt, (51)

where we have used the fact that instanton lines link fixed points lying on zero-energy invariant lines. Since

 H1(xh(t−t0),t)=∑m≠nβ1mn(t)Jmn(t−t0),Jmn(t−t0)=1m!(p(t−t0)n−p(t−t0)m)q(t−t0)m, (52)

we can integrate by parts to get

 M(t0)=−∫∞−∞∑m≠nddtβ1mn(t)Jmn(t−t0)dt=∫∞−∞∑m≠nβ1mn(t)ddtJmn(t−t0)dt. (53)

This allows us to rewrite, at least formally, the Melnikov function as a Fourier series,

 M(t0)=∞∑k=−∞e2πikt0/T∫∞−∞∑m≠nβ1mn,kddsJmn(s)e2πiks/Tds, (54)

after the change of the integration variable . Its mean value is thus given by

 1T∫T0M(t0)dt0=∫∞−∞∑m≠nβ1mn,0ddsJmn(s)ds=0. (55)

In order to check the validity of Fourier series (54) we need to bound the norm

 ∫∞−∞∣∣∣dJmndt∣∣∣dt. (56)

This is not a difficult task because both and must be continuously differentiable along the connection and with finite limits when , which must coincide with their values at the fixed points. Furthermore, the derivatives and decrease to zero exponentially at infinity due to the hyperbolic character of the fixed points. This implies that the norm (56) is necessarily finite.

It is actually simpler to show that the Melnikov function has zero mean. Equation (53) can be cast into the form

 M(t0)=−ddt0∫∞−∞∑m≠nβ1mn(t+t0)Jmn(t)dt, (57)

which, together with the periodicity of the functions , directly provides the desired result. However, the Fourier series setting favors the identification of the conditions under which the Melnikov function vanishes, as Eq. (43) in the last section, which cannot be straightforwardly extracted from this expression. So we believe that the Fourier series representation will be more helpful in applications to concrete reaction schemes.

In conclusion, we have shown that the Melnikov function of a perturbed instanton connection can always be expanded in Fourier series, provided that it links hyperbolic fixed points. Furthermore, it is always a zero mean function, which implies that it either crosses zero or vanish identically. So for any reaction Hamiltonian we have the following result: given an instanton connection linking two hyperbolic fixed points, any small time-periodic perturbation of the rates will preserve the existence of the connection if the corresponding Melnikov function is not identically zero. In the other case, we would have to study the behavior of the Melnikov distance at higher orders, in order to be able to give a rigorous conclusion (23).

To finish, let us note the difficulties in extending this result and proving a general theorem about the persistence of the instanton connections even when the Melnikov function is identically zero. In this case, we could recall the expansion of the Melnikov distance in terms of higher-order Melnikov functions (23). One would be tempted to use an equivalent argument to that of the present section to try to show that an arbitrary-order Melnikov function either has zero mean or vanish identically. This would imply in turn that either some function in the expansion is nonzero, and thus the connection is preserved by means of a transversal crossing of the stable and unstable manifolds, or the whole series becomes identically zero. While this can suggest, at first sight, that the Melnikov distance is also identically zero in this case, we cannot rely on a rigorous argument to prove so, since the perturbative expansion is not analytic in the small parameter in general situations. Furthermore, we can always find a family of transformations (for instance, by continuing the chain of reparametrizations in Eqs. (45) and (46)) able to nullify an arbitrary order Melnikov function. So a full proof would need to handle these special perturbations, and this is complicated by the lack of analyticity of the series expansion, despite the probable fact that they are built by combinations of reparametrizations, constant shifts of the parameter values, or other trivial transformations preserving the phase-space topology. Also, proving that the higher-order Melnikov functions have zero mean is not a trivial fact. These imply nonlinear combinations of the external perturbations and the corresponding mixture of Fourier modes that complicates the development of a simple form like Eq. (54) in these cases. We illustrate this situation with the second-order Melnikov function in the Appendix. Anyway, perturbations requiring a higher-order Melnikov analysis do not constitute the generic case, and we believe that the majority of the physical situations could be analyzed within the first-order, or at most second-order, formalism.

## Iii Discussion

So far we have concentrated in showing that the instanton connections persist when the chemical system is temporally forced by a weak perturbation. We devote this section to explain the physical consequences of this fact. As we already argued, the topology of the phase space describes the physics of the system, so if its topology is preserved when the system is forced, this means that the qualitative behavior of the system is preserved too. The disappearance of an instanton connection would mean that this qualitative behavior is modified, but how strongly? These connections are optimal paths linking different physical states (given by the fixed points of the dynamical system), but the absence of one connection communicating two states does not mean that the system cannot go from one to the other. It only means that there is not an optimal way of going. We will show the importance of this fact using a toy model of biological relevance in plankton modeling. Suppose we have the reactions

 A\lx@stackrelγ⟶A+A,A\lx@stackrelγ⟶∅ (58)

occurring at the same constant rate . This set of equations has been used to model plankton patchiness in some occasions (24); (25); (26). We can straightforwardly write the Schrödinger equation for the generating function (6); (27),

 ∂tG=γ(p−1)2∂pG, (59)

 H=γ(p−1)2q. (60)

The corresponding dynamical system is highly degenerate, with two invariant lines and , with all their points being fixed. This means that if we start with an initial condition , for some , then this will be the solution for all times. As is totally clear, there are no instanton lines linking the mean-field line with the absorbing-state line. Let us now concentrate directly in the exact PDE (59) instead of the approximate dynamical system. This is a linear, first-order PDE which can be easily solved along characteristics (6)

 G(p,t)=G0[1−1−p1+(1−p)γt], (61)

and due to normalization, we have the long-time asymptotics

 limt→∞G(p,t)=limt→∞G0[1−1−p1+(1−p)γt]=G0(1)=1, (62)

which corresponds exactly to the probability distribution . This means that in the infinite-time limit, regardless of the initial condition, the system will be in the absorbing state with zero particles. So, as we have seen, the absence of instanton connections does not stop the system from undergoing an extinction transition. One can figure out how this happens representing the process (58) with the Itô stochastic differential equation (25); (3)

 dρ=√2γρdW, (63)

where denotes a Wiener process. This equation describes Brownian motion with state-dependent diffusion and tells us that the transition to extinction is very different in this case: no optimal trajectory is chosen; instead, the state with zero particles is reached after performing a random walk from the initial condition. This will not necessarily happen in every situation; actually, state-dependent Brownian motion is one of the simplest possibilities. More complicated Hamiltonians with powers of the momentum above the second exhibit properties corresponding to non-Gaussian statistics (6). In general, the random walk will be more complex than simple Brownian motion, but all the cases will have in common the absence of an optimal way of going from one physical state to another.

The instanton connections are a useful tool that allows us to calculate the frequency of rare fluctuations that drive the system among different states. The mean transition time is obtained, at exponential order, computing the action along the instanton connection (11). If our perturbation is pure time reparametrization, as the second perturbation in Eqs. (45) and (46), and the function periodic, a straightforward computation reveals that the action along the heteroclinic orbit is exactly the same as for the unperturbed system. More general perturbations are not tractable analytically, and we would have to rely on a numerical treatment of the problem. The heteroclinic connection can be obtained, for instance, by means of a shooting method and the action computed as a numerical integral on it. In any case, we expect that small, well-behaved, periodic perturbations will yield transition times of the same order of magnitude as the unperturbed system. Instanton persistence is important as it allows approximating this sort of nonautonomous systems by their autonomous counterparts in the first instance and serves as a justification of the quasistationary approximation for slow signals. Perturbations that are not periodic might have a stronger effect on the system, and in particular, singular enough perturbations could change the phase-space topology; however, it is rather difficult to conceive physical situations that give rise to such a perturbation. Also, the mathematical framework changes considerably, since the absence of a Poincaré map does not allow the definition of the instanton as a connection among fixed points of this map, like in the periodic case. In any case, even if the instanton connections were broken, but the physical states did not drift apart form each other, transitions will be possible if we wait long enough, due to the uncorrelated fluctuations the system is subject to (28), but however, there will be no optimal paths linking those states.

## Iv Conclusions and Outlook

In this work we have studied the persistence of instanton connections in chemical systems, when we promote the reaction rates from constants to functions of time. A set of chemical reactions can be mapped, under the eikonal approximation, onto a dynamical system, the fixed points of which denote the possible states where the chemical system can be found. Connections among fixed points denote optimal transition paths communicating different states. The persistence of these instanton lines indicates the robustsness of the physical processes that are taking place in the reactor and shows us that the qualitative behavior of the system will be similar when subject to small nonautonomous perturbations. We have shown our results with one particular model, and we have stated them in a more general setting when it has been possible. The main theoretical tool that we have employed is the Melnikov function, which has proven itself very useful in dealing with this type of problems.

One of the directions in which we would like to extend the theory is the problematic of connections linking nonhyperbolic, or one nonhyperbolic and one hyperbolic, fixed points. The physical motivation comes from the recent classification of phase transitions in reaction-diffusion models, which identifies critical points of the physical system with the nonhyperbolic fixed points at a bifurcation threshold of the dynamical system (15). While it is very difficult to see what happens in the general situation, we know that the mean value of the Melnikov function is zero if the mean value of all the external perturbations is zero, as is shown by Eq. (55), provided it is well defined in this case. This suggests to us that in a broad number of situations the topology of the phase space will remain unchanged if the mean value of all the external perturbations vanishes, even when some of its fixed points were nonhyperbolic. It is, however, very difficult to prove a result in that direction, since perturbating a system in a critical state could have unpredictable consequences. Furthermore, the existence of the Fourier series expansion of the Melnikov function is no longer guaranteed, as the loss of hyperbolicity implies that we do not have an exponential decay at infinity of the absolute value in Eq. (56).

If some of the perturbations have a nonzero mean, then the problem becomes much more difficult. The system will be moved out of criticality at some of its points, and its behavior will be more similar to that of the (partially) noncritical phase corresponding to the autonomous system with the parameters retuned according to the mean value of the external perturbations. Unfortunately, we cannot map this situation to that described in Sec. II, because in this case we have no control on the amplitude of the perturbation, which may be comparable to the distance to some of the critical points. So the problem becomes genuinely nonperturbational, and we can no longer employ a technique like the Melnikov function. In this case it is difficult to deal even with very slow perturbations, since this type of settings favors the appearance of soft modes, which rules out the possibility of treating the external signals adiabatically (22).

Another problem we would like to deal with is the case of several perturbations with different periods. The simplest case is that in which the ratio of the periods of all possible different pairs of perturbations that are affecting the system is a rational number. In this case, we say that the perturbations are commensurable, and we can find a period which is common to all of them. Once obtained, we have reduced our problem to the one studied in Sec II. The case of incommensurable perturbations is, of course, more complex, as it implies that the dynamics of the perturbed system cannot be reduced to the dynamics of a periodic Poincaré map. The fixed points of the unperturbed system give rise to quasiperiodic solutions of the perturbed problem, with their invariant manifolds being quasiperiodic as well (29).

There are many other questions that remain to be answered. One is the possible appearance of chaotic behavior induced by internal stochasticity. As noted in (12), for two reacting species we have a four-dimensional eikonal Hamiltonian, which will give rise, in general, to a three-dimensional dynamical system on some Riemannian manifold. In this case, chaos, unlike in the two-dimensional mean-field dynamical system, is indeed possible. But chaos is also possible in the case of one reacting species obeying reaction rules with time-dependent rates. In this situation, ”energy” is no longer conserved and the nonautonomous dynamical system, without integrals of motion, becomes effectively three dimensional. Indeed, situations like the one plotted in Fig. (2), that is, when the heteroclinic connection is preserved due to the intersection of the stable and unstable manifolds, give rise to complex dynamical scenarios. These include the appearance of Smale horseshoes and even the presence of strange attractors related to the creation and destruction of such horseshoes (30); (31); (32). Although the generation of chaotic behavior of this type can be attractive from a nonlinear dynamics point of view, we are not aware, at this point, of its possible physical meaning.

Of course, studying multispecies reactions, both autonomous and nonautonomous, is a very interesting problem that deserves further efforts. The nonautonomous situation implies the technical difficulty of extending the Melnikov function theory to a higher dimensionality. The ideas developed for three-dimensional systems (33); (34) could perhaps be adapted for studying four- or higher-dimensional problems and to try to obtain in this way the results that we would need to understand the more general multispecies reactions. Also, as we have already pointed out, the Melnikov function is a perturbative result. It allows us to treat arbitrary frequency but necessarily small perturbations. To fully understand the dynamics of nonautonomous chemical reactions, even in the one species case, we would need a nonperturbational result. With it we could try to address questions such as the forcing of systems near criticality or the appearance of soft modes. As we can see, chemical kinetics presents a very complex phenomenology, and it is challenging from a mathematical point of view. A deep understanding of the physics of these nonequilibrium systems will presumably imply the parallel development of powerful methodological techniques.

## Acknowledgments

The authors gratefully acknowledge input from Ernest Fontich and Carles Simó. This work has been partially supported by the Ministerio de Educación y Ciencia (Spain) through Projects Nos. EX2005-0976, FIS2005-01729, and MTM2005-02094.

*

## Appendix A Second-order Melnikov function

Consider the two-dimensional dynamical system

 ˙x=f(x)+ϵg(x,t)+ϵ2h(x,t), (64)

where , , , and are sufficiently regular functions. The second-order Melnikov function reads

 M2(t0)=∫∞−∞f(xh(t))∧[2−1D2f(xh(t)){x1h(t+t0,t0)}2 +Dg(xh(t),t+t0)x1h(t+t0,t0)+h(xh(t),t+t0)]dt, (65)

where and denote the Hessian and Jacobian of and , respectively, is the solution parametrizing some given heteroclinic connection, and is the solution to the variational equation

 ˙x1h(t,t0)=Df(xh(t−t0))x1h(t,t0)+g(xh(t−t0),t). (66)

We now assume that this dynamical system is a reaction Hamiltonian system and both and come from small nonautonomous perturbations of the Hamiltonian parameters, just like in Sec. II. Then, because Eq. (66) is linear, its solution will depend linearly on the Fourier modes of the time-dependent parameters evaluated at . But this solution appears quadratically in Eq. (65) and also multiplying the Jacobian of , which depends linearly, by assumption, on the same Fourier modes. This quadratic dependence complicates recasting the second-order Melnikov function into a form similar to Eq. (54), which would help enormously for calculating its mean value.

### References

1. E. A. Mastny, E. L. Haseltine, and J. B. Rawlings, J. Chem. Phys. 125, 194715 (2006).
2. T. B. Kepler and T. C. Elston, Biophys. J. 81, 3116 (2001).
3. C. Escudero, J. Buceta, F. J. de la Rubia, and K. Lindenberg, Phys. Rev. E 69, 021908 (2004).
4. M. Scalerandi, B. C. Sansone, C. Benati, and C. A. Condat, Phys. Rev. E 65, 051918 (2002).
5. J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields (Springer-Verlag, New York, 1983).
6. C. W. Gardiner, Handbook of Stochastic Methods (Springer-Verlag, Berlin, 1996).
7. N. G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amsterdam, 1981).
8. M. Doi, J. Phys. A 9, 1479 (1976).
9. L. Peliti, J. Phys. (Paris) 46, 1469 (1985).
10. M. I. Dykman, E. Mori, J. Ross, and P. M. Hunt, J. Chem. Phys. 100, 5735 (1994).
11. V. Elgart and A. Kamenev, Phys. Rev. E 70, 041106 (2004).
12. M. Assaf and B. Meerson, Phys. Rev. E 74, 041115 (2006).
13. M. Assaf and B. Meerson, Phys. Rev. Lett. 97, 200602 (2006).
14. J. Zinn-Justin, Quantum Field Theory and Critical Phenomena (Claredon Press, Oxford, 1989).
15. V. Elgart and A. Kamenev, Phys. Rev. E 74, 041101 (2006).
16. A. K. Horváth, M. Dolnik, A. P. Muñuzuri, A. M. Zhabotinsky, and I. R. Epstein, Phys. Rev. Lett. 83, 2950 (1999).
17. J. Buceta, C. Escudero, F. J. de la Rubia, and K. Lindenberg, Phys. Rev. E 69, 021906 (2004).
18. M. I. Dykman, T. Horita, and J. Ross, J. Chem. Phys. 103, 966 (1995).
19. C. Escudero, Phys. Rev. E 74, 010103(R) (2006).
20. D. A. McQuarrie, C. J. Jachimowski, and M. E. Russell, J. Chem. Phys. 40, 2914 (1964).
21. V. Melnikov, Trans. Moscow Math. Soc. 12, 3 (1963).
22. D. Ryvkine, M. I. Dykman, and B. Golding, Phys. Rev. E 69, 061102 (2004).
23. Z. Zhi-Fen and L. Bao-Yi, Annali di Matematica pura ed applicata CLXI, 181 (1992).
24. Y-C. Zhang, M. Serva, and M. Polikarpov, J. Stat. Phys. 58, 849 (1990).
25. R. Adler, Monte Carlo simulation in oceanography, Proceedings of the 10th ’Aha Huliko’a Hawaiian Winter Workshop, University of Hawaii at Manoa (1997).
26. W. R. Young, A. J. Roberts, and G. Stuhne, Nature (London) 412, 328 (2001).
27. F. Baumann, M. Henkel, M. Pleimling, and J. Richert, J. Phys. A: Math. Gen. 38, 6623 (2005).
28. K. Lindenberg and B. J. West, J. Stat. Phys. 42, 201 (1986).
29. S. Chow, J. Hale, and J. Mallet-Paret, J. Diff. Equations 37, 3 (1980).
30. S. Smale, Bull. Amer. Math. Soc. 73, 747 (1967).
31. L. Mora and M. Viana, Acta Math. 171, 1 (1993).
32. A. Pumariño and J. Á. Rodríguez, Coexistence and Persistence of Strange Attractors, Lecture Notes in Mathematics, Vol. 1658 (Springer-Verlag, Berlin, 1997).
33. J. Á. Rodríguez, Arch. Ration. Mech. Anal. 93, 81 (1986).
34. F. Costal and J. Á. Rodríguez, Arch. Ration. Mech. Anal. 104, 185 (1988).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters