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:

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 N reacting species is probably given by the mean-field equations, an N−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 2N dimensions for N reacting species. The additional N 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 (2N−1)−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.

Figure 1: Phase space for the branching and annihilating particle problem. The parameter values are σ=λ=1.

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 A, 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 n reagents at time t reads

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 ~t=∫t0σ(τ)dτ)

˙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 p=1, q=0, and

q=2σp/λ1+p.

(8)

Furthermore, we know that these three lines connect the fixed points (0,0), (1,0), and (1,σ/λ) in the (p,q) 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 p=1 and q=0 remain invariant zero-energy lines of the system; however, the explicit time dependence of the Hamiltonian prevents the conservation of the energy H 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 p=1 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 q=0. 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

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?

Figure 2: Sketch of the phase space for the branching and annihilating particle problem, assuming that the stable and unstable branches intersect. The heteroclinic connection goes through all the intersections of these branches.

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)

x=(x1,x2), f=(f1,f2), with two hyperbolic fixed points xa and xb linked by a heteroclinic connection, which is parametrized by the system solution xh(t−t0) for initial time t0. Assuming that the motion goes from xa to xb, this connection is formed by a branch of the unstable manifold of xa, which totally overlaps with a branch of the stable manifold of xb. Let us now consider the perturbed version of this problem,

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

(15)

where the perturbation g(x,t) is time periodic with period T, amplitude ϵ small enough and sufficiently regular. The dynamics of this nonautonomous system is given by the associated Poincaré map Pϵ, which maps every initial condition point x(0) with the corresponding value
of the solution x(T) after one period has elapsed (5). The hyperbolic fixed points of the unperturbed system, xa and xb, are hyperbolic fixed points of the Poincaré map P0. Since the map Pϵ is a perturbation of P0, these points have a continuation, xϵa and xϵb, as hyperbolic fixed points of Pϵ, 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 xa and
xb, 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 d(ϵ,t0)=ϵM(t0)+O(ϵ2), 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 f and g. 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 xh(t0). A change of sign of M(t0) means that there exists some t0 such that d(ϵ,t0)=0, implying the existence of a solution xϵh(t) of Eq. (15) defining a heteroclinic connection among the two hyperbolic fixed points xϵa and xϵb of the Poincaré map corresponding to Eq. (15), say,

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

(18)

Figure 3: Sketch of the phase space for the branching and annihilating particle problem, assuming that the stable and unstable branches do not intersect.

Let us now consider the branching and annihilating system, Eqs. (6) and (7), with constant reaction rates. The instanton connection (the separatrix linking (1,σ/λ) to (0,0)) 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 Pϵ still has the origin as hyperbolic fixed point. Another hyperbolic fixed point of this map, corresponding to (1,σ/λ) when ϵ=0, is (1,qϵ), where qϵ is obtained through the solution of Eq. (23) for p=1. 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 q(0)=q(T) to this solution to find

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

(25)

Figure 4: Sketch of the phase space for the branching and annihilating particle problem, assuming that the stable and unstable branches do not intersect.

To determine if the unstable manifold of (1,qϵ) intersects the stable manifold of (0,0) we substitute the corresponding values of f and g 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 s=σ(t−t0) 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 sh standing for the hyperbolic sine. It is easy to see that ~ϕ(s) has zero mean,

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 δn0 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 f, g, and xh), and with zero mean, implying that it either crosses zero or vanish identically. In the first case, we can claim that the point (1,qϵ) is connected to the origin, in the second, we cannot, due to the possible presence of small terms O(ϵ2). However, we can see that the second case only happens when h(t) is constant. Indeed, integrating by parts and changing variables x=es we obtain

and one can see that it only vanishes if n=0 or n→±∞, for σ and T positive finite real numbers. Hence, the Melnikov function is identically zero if and only if an=0 for all n≠0 in Eq. (31) or, what is the same, if and only if h(t) is constant.

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

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

(43)

for some constant c. If c=0, 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 σ1 and
λ1 are chosen constant, then for any value of c, 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 h(t)=c is a combination of these two, translation and reparameterization,

σ(t)

=

(σ+ϵσ1)[1+ϵϕ(t)],

(45)

λ(t)

=

(λ+ϵλ1)[1+ϵϕ(t)],

(46)

for some function ϕ(t), up to terms O(ϵ2). 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 βmn(t) are supposed to have the same period T (the case of perturbations with different periods will be discussed below). Assuming the form of the perturbation βmn(t)=β0mn+ϵβ1mn(t) allows us to split the Hamiltonian into two terms H=H0+ϵH1. 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

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 s=t−t0. 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 p(t) and q(t) must be continuously differentiable along the connection and with finite limits when t→±∞, which must coincide with their values at the fixed points. Furthermore, the derivatives ˙p and ˙q 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 β1mn(t), 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.

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)

so the reaction Hamiltonian reads

H=γ(p−1)2q.

(60)

The corresponding dynamical system is highly degenerate, with two invariant lines p=1 and q=0, with all their points being fixed. This means that if we start with an initial condition (q=n0,p=1), for some n0>0, 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 Pn=δn0. 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 W 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 ϕ(t) 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.

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.

where D2f and Dg denote the Hessian and Jacobian of f and g, respectively,
xh(t) is the solution parametrizing some given heteroclinic connection, and x1h 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 g and h 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 t0. But this solution appears quadratically in Eq. (65) and also multiplying the Jacobian of g, 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.