text:
# Dynamics of FitzHugh-Nagumo excitable systems with delayed coupling

## 1 Introduction

Excitability is a common property of many physical and biological systems. Since the work of Hodgkin and Huxley , and the development of the basic mathematical model by FitzHugh  and Nagumo  the reported research on the subject has grown enormously. As for a general review we cite just the classic references , and the references , for examples of a recent physical, and ,  for neuro-biological applications. For instance, a single neuron displays excitable behavior, in the sense that a small perturbation away from its quiescent state, i.e. a stable stationary value of the cross membrane potential, can result in a large excursion of its potential before returning to quiescent. Such generation of a single spike in the electrical potential across the neuron membrane is a typical example of the excitable behavior. Many other cells, besides neurons, are known to generates potential spikes across their membrane. Such excitable units usually appear as constitutive elements of complex systems, and can transmit excitations between them. The dynamics of the complex system depends on the properties of each of the units and on their interactions. In biological, as well as physical, applications the transmission of excitations is certainly not instantaneous, and the representation by non-local and instantaneous interactions should be considered only as a very crude approximation. For example, significant delays of more than $`4\%`$ of the characteristic period of the $`40Hz`$ frequency oscillations of the brain neurons, occur during the nerve conduction between the neurons less then $`1mm`$ apart ,.

This paper is devoted to an analyzes of a small lattice of a particular type of excitable systems, with a finite non-zero duration of the transfer of the excitations between the neighboring units. Despite its relevance and a large amount of related research ( to be summarized and discussed in the last section ) excitable systems with time-delayed coupling have not been sufficiently studied. We shall be particularly interested in the bifurcations that turn on and turn off the oscillatory behavior as the coupling constant and the small time-lag are varied.

In the introduction we formulate the model that is to be analyzed, and then briefly preview our results and discuss the context of our work.

As a model of each of the excitable units we shall use the paradigmatic example of the FitzHugh-Nagumo system in the form, and for the parameter range, when the system displays the excitable behavior. The dynamical equations of the single uncoupled excitable unit are :

$`\dot{x}`$ $`=`$ $`x^3+(a+1)x^2axy,`$
$`\dot{y}`$ $`=`$ $`bx\gamma y.`$ (1)
where $`a,b`$ and $`\gamma `$ are positive parameters. In the original interpretation of (1), as a model of the neuronal excitability, $`x`$ represents the trans-membrane voltage and the variable $`y`$ should model the time dependence of several physical quantities related to electrical conductances of the relevant ion currents across the membrane. In the model $`x`$ behaves as an excitable variable and $`y`$ is the slow refractory variable.

The particular form (1) of the FitzHugh-Nagumo model does not admit periodic solutions for any values of the parameters. Furthermore we shall restrict our analysis to the range of the parameter values where the system exhibits excitability, with only one attractor in the form of a stable fixed point at the origin. For this to be the case, $`b`$ and $`\gamma `$ should be of the same order of magnitude and considerably smaller than $`a`$ (see section 2). We refer to the system (1) in this range of parameters as the excitable FitzHugh-Nagumo model. On the other hand, the minimal modification of (1) that renders a system which could have a stable limit cycle is obtained by adding to the first equation an external constant current $`I`$ of a prescribed intensity. We shall refer to such a system with the stable limit cycle as the FitzHugh-Nagumo oscillator as opposed to the excitable FitzHugh-Nagumo (1).

The full system is a one-dimensional lattice of N identical excitable units of the form (1), given by the equations of the following type:

$`\dot{x}_i`$ $`=`$ $`x_i^3+(a+1)x_i^2ax_iy_i+cF(x_{i1}^\tau ,x_i,x_{i+1}^\tau ),`$
$`\dot{y}_i`$ $`=`$ $`bx_i\gamma y_i,i=1,\mathrm{}N,`$ (2)
where

$$x_i^\tau (t)x_i(t\tau ),$$

and $`\tau `$ is a fixed time lag and $`c`$ is the coupling constant. General form of the coupling term will be specified later.

Local stability near the rest state of (2), and global dynamics like existence of stable limit cycles, and the properties of the oscillations on such a cycle, do depend on the coupling function. However, we shall see that local properties and even the global dynamics are qualitatively the same for a large class of coupling functions, that are dependent only on the voltages of the neighbors, like for example:

$$F(x_{i1}^\tau ,x_i,x_{i+1}^\tau )=f(x_{i1}^\tau )+f(x_{i+1}^\tau ),f(x)=\mathrm{tan}^1(x).$$

On the other hand, diffusive coupling, i.e. proportional to $`x_i(t)x_{i1}(t\tau )`$ implies different properties of the global dynamics. Furthermore, important dynamical phenomena that occur for $`N=2`$, happen also for $`N>2`$. In fact, most of our results will be derived by considering first the system with only two coupled units, and then checking the conclusions in the case of medium $`N>2`$ by numerical computations.

It is well known, and often used, fact that the time-delay could destabilize a stationary point and introduces oscillatory behavior. Also, networks of oscillatory units with delayed coupling have been analyzed before. The studied oscillatory systems could be roughly divided into those where the oscillatory units are general limit cycle oscillators, say near the Hopf bifurcation, (for example: ,, ), phase-coupled phase oscillators ,, ,,,), or the relaxation oscillators (for example: ,) typical in the neuro-biological applications ,,. In the later case the form of the coupling takes, more or less, into the account the properties of real synaptic interactions between the neurons ,.

In the last section we shall more systematically compare, the system (2) and our results with several similar or related models. In the introduction, we should like to point out that the major part of our analysis deals with the system of coupled excitable units, and the system of FitzHugh-Nagumo oscillators with the same coupling is mentioned only in order to stress the differences. On the other hand, a sufficiently strong instantaneous coupling (time lag equal to zero) between the excitable (not oscillatory) units can introduce the oscillatory solutions. This phenomenon has been known already to Turing  and was studied by Smale  and Johnson . As we shall see, for such sufficiently strong coupling, a time-lag which is small on the scale of the interspikes or refractory period, induces drastic qualitative changes in the dynamics. Phenomena like death of oscillations, bi-stable excitability, and transitions between symmetric in-phase and non-symmetric, phase-shifted asymptotic oscillations, all occur in the system (2) as the time-delay is varied. On the other hand, dynamics of the coupled FitzHugh-Nagumo oscillators with the same type of coupling is quite different.

The results of our study are presented as follows. Sections 2 and 3 are concerned with the system with just two excitable units. Analytic results about the codimension 1 bifurcations of the stationary solutions are given in detail, in section 2, for a specific common type of coupling such that there are difference between coupled excitable and coupled oscillatory units. Other types of coupling are briefly discussed. Numerical analyzes of global dynamics and in particular of the periodic solutions and their bifurcations are presented in section 3. Here we also point out some of the differences between coupled excitable systems vs. the oscillators. In section 4 we demonstrate, by direct numerical computations, that the phenomena analyzed in sections 2 and 3 for $`N=2`$ occur also in a similar way in the system consisting of $`N>2`$ identical units. Conclusions, discussion and comparison with related works are given in section 5.

## 2  Two coupled units: Local stability and bifurcations

In this section we study stability and bifurcations of the zero stationary point of only two coupled identical FitzHugh-Nagumo excitable systems, given by the following equations:

$`\dot{x}_1`$ $`=`$ $`x_1^3+(a+1)x_1^2ax_1y_1+cf(x_2^\tau ),`$
$`\dot{y}_1`$ $`=`$ $`bx_1\gamma y_1,`$
$`\dot{x}_2`$ $`=`$ $`x_2^3+(a+1)x_2^2ax_2y_2+cf(x_1^\tau ),`$
$`\dot{y}_2`$ $`=`$ $`bx_2\gamma y_2,`$ (3)
where the coupling function satisfies

$$f_0=0,f_0^{}=\delta >0,$$
(4)
and the subscript $`0`$ denotes that the function is evaluated at $`(x_1,x_2,y_1,y_2)=(0,0,0,0)`$. If fact, the first condition is not crucial, and is introduced only for convenience.

A single neuron:

Consider first one of the units in the case of zero coupling. Point $`(x,y)=(0,0)`$ is an intersection of the qubic $`\dot{x}`$ nullcline and the linear $`\dot{y}`$ nullcline for any value of the parameters $`a,b,\gamma `$, so that it is always a stationary point. Furthermore it is always a stable stationary point, that could be a node, if $`(a\gamma )>2\sqrt{b}`$, or a focus, when $`(a\gamma )<2\sqrt{b}`$. There could be one more (non-generic case) and two more stationary points, but we shall restrict our attention to the case when $`(0,0)`$ is the only stationary solution. This is the case if $`4b/\gamma >(a1)^2`$. We shall make no further assumptions as to the nature of the stable stationary point $`(0,0)`$, but, as we shall see, some of the typical behavior of the delayed coupled systems will be lost in the singular limit $`b0,\gamma 0`$, and difficult to observe very close to this limit. The particular form of the FitzHugh-Nagumo model, with no external current and $`(0,0)`$ stationary point does not possess periodic solutions for any values of the parameters. However, there are solutions that start in a small neighborhood of $`(0,0)`$, quite rapidly go relatively far away, and then approach back on to the stationary point, (see figure 5a). Such solutions represents typical excitable behavior. The excitability that is displayed by the FitzHugh-Nagumo system is of, the so called, type II , in the sense that there is no clear-cut threshold in the phase space between the excitable orbits and the orbits that return quickly, and directly, to the rest state. In fact, there are orbits which continuously interpolate between the two types of orbits. However, as the parameters $`b`$ and $`g`$ are decreased, as compared with a fixed $`a`$, the excitable behaviour quite rapidly (but continuously) become dominant. We shall always use such values of the parameters that the excitable behaviour is clearly demonstrated, for example $`b>\gamma ^2`$ and $`ab,a\gamma `$.

In order to stress typical properties of the excitable, but not oscillatory, systems, we shall also need a convenient system with a stable oscillatory behaviour. Such a system is obtained by adding an external (say, constant) current $`I`$ to the $`\dot{x}`$ equation (or to the other one) of the FitzHugh-Nagumo model. The constant current shifts the intersection of the two nullclines, and if it is such that the intersection lies on the part of the $`\dot{x}`$ nullcline with a positive slope, then the stationary point is unstable and there is a stable limit cycle. The limit cycle is born in the supercritical Hopf bifurcation of the stationary solution. The limit cycle is of approximately circular shape only if $`I`$ is quite close to the critical value $`I_0`$, and then the amplitude is of the order of $`\sqrt{II_0}`$. Otherwise it has the shape typical for relaxation oscillators.

Instantaneously coupled identical units:

As the next step, let us fix the parameters $`a,b`$ and $`\gamma `$ such that each of the units displays the excitable behavior, and consider the coupled system but with the instantaneous coupling $`\tau =0`$. Point $`(x_1,y_1,x_2,y_2)=(0,0,0,0)`$ represent a stationary solution, and its local stability is determined by analyzing the corresponding characteristic equation

$$[(a+\lambda )(\gamma +\lambda )+bc\delta (\gamma +\lambda )][(a+\lambda )(\gamma +\lambda )+bc\delta (\gamma +\lambda )]=0,$$
(5)
The sign of the real parts of the four eigenvalues

$`2\lambda _{1,2}`$ $`=`$ $`(a+\gamma c\delta )\pm \sqrt{(a\gamma c\delta )^24b},`$
$`2\lambda _{3,4}`$ $`=`$ $`(a+\gamma +c\delta )\pm \sqrt{(a\gamma +c\delta )^24b},`$ (6)
determines the stability type of the trivial stationary point. If $`(a\gamma )>2\sqrt{b}`$ the point is stable node-node for $`0<c<(a\gamma 2\sqrt{b})/\delta `$, and if $`c`$ is larger the eigenvalue $`\lambda _{1,2}`$ becomes complex, and the point becomes stable focus-node. Otherwise, if $`(a\gamma )<2\sqrt{b}`$ the point is stable focus-focus for $`0<c<(2\sqrt{b}a+\gamma )/\delta `$ and for larger $`c`$ the eigenvalue $`\lambda _{3,4}`$ becomes real and the point is again stable focus-node.

Thus, whatever the stability type of the stationary point in the uncoupled case might be, there is the corresponding value of the coupling constant $`c`$ such that the point becomes focus-node. Then, the complex pair of the eigenvalues $`\lambda _{1,2}`$ correspond to the eigenspace given by $`x_1=x_2`$ and $`y_1=y_2`$. In such a situation the damped oscillations of the two units interfere synchronously, and at some still larger $`c_0`$ given, in both the cases, by

$$c_0=\frac{a+\gamma }{\delta }$$
(7)
where

$$sgn\left(\frac{d\mathrm{Re}\lambda _{1,2}}{dc}\right)_{c=c_0}=sgn\left(\delta /2\right)>0,$$
(8)
the point goes through a direct supercritical Hopf bifurcation. As the result, for small $`ϵ=cc_0>0`$, the stationary solution acquires an unstable two-dimensional manifold, with a stable limit cycle in it. The unstable manifold is in fact a plane given by the equations $`x_1=x_2x`$ and $`y_1=y_2y`$, independently of the form of the coupling function as in (3). Oscillations on the limit cycle are synchronous, with the linear frequency $`\omega =\sqrt{b\gamma ^2}`$, and symmetrical. In this paper, by synchronous we mean coherent in-phase oscillations, and by symmetrical we mean that $`x_1(t)=x_2(t)`$.

The dynamics on the unstable manifold for small $`ϵ`$ is given by the normal form of the Hopf bifurcation

$$\dot{r}=dϵr+\alpha r^3+O(ϵ^2r,ϵr^3,r^5),\dot{\theta }=\omega +eϵ+\beta r^2+O(ϵ^2,ϵr^2,r^4),$$
(9)
where $`\omega =\sqrt{b\gamma ^2}`$, $`d=\delta /2`$, $`e=\gamma \delta /2\omega `$ and $`r`$ and $`\theta `$ are the polar coordinates

$$x=r\mathrm{sin}\theta ,y=r\mathrm{cos}\theta .$$

The parameters $`\alpha `$ and $`\beta `$ depend on the particular form of the coupling function. For example, in the case the coupling function is $`f(x)=\mathrm{tan}^1(x)`$, then

$$\alpha =\frac{3+c_0}{8}+\frac{(a+1)^2\gamma }{4\omega ^2},\beta =\frac{(3+c_0)\gamma }{8\omega }\frac{(a+1)^2(5\gamma ^2+2\omega ^2)}{12\omega ^3}$$

The limit cycle of (9) is a good approximation only for $`ϵ`$ quite small. However numerical analysis shows that the limit cycle remains a global attractor in the full four-dimensional phase space of the system (3) with no time-delay for a large range of $`c>c_0`$ values, where the approximation by the Hopf normal form (9) is no more valid. Thus there is a range of values of the coupling parameter $`c`$ where the system (3) with no time-delay behaves as a system of two coupled identical limit cycle oscillators. Properties of the oscillations on the limit cycle depend on $`c`$. Perturbation analyzes for $`ϵ`$ small, or the numerical analyzes for larger $`c`$, show that oscillations on the limit cycle are synchronous and symmetrical. In figure 1, we illustrate the limit cycles in the coupled excitable systems with no time-delay. The figure illustrates oscillatory dynamics of both units since on the limit cycle $`x_1(t)=x_2(t)`$ and $`y_1(t)=y_2(t)`$. Although the limit cycles deform continuously with $`c`$, the deformation from the small Hopf circle all the way up to the large limit cycle of the shape like for the relaxation oscillators, happens on a small interval of the values of $`c`$, smaller than $`3\%`$ of the interval $`(c_0,c_1)`$.

Further increase of $`c`$, still with $`\tau =0`$, leads to a bifurcation of the stationary solution and of the limit cycle. For $`c>(a\gamma +2\sqrt{b})/\delta `$, there is a pair of real positive and a pair of real-negative eigenvalues at the trivial solution. The limit cycle disappears at some still larger $`c_1`$ when there appear other stable stationary solutions of (3) (with $`\tau =0`$). This value of the coupling constant $`c=c_1`$, when there appears nonzero stable stationary solution, obviously depends on the coupling function.

In conclusion, there are three qualitatively different types of dynamics of the instantaneously coupled excitable systems. For $`0<c<c_0`$ the coupled system behaves as a pair of excitable units, while for $`c_0<c<c_1`$ the system behaves as a pair of identical limit cycle oscillators. For $`c>c_1`$ there appears a nontrivial stable stationary state. However, we shall be interested in the influence of time-delay only when the coupling constant is in the range $`c(0,c_1)`$, i.e. when the instantaneously coupled system behaves either as excitable $`c<c_0`$ or as oscillatory $`c>c_0`$.

Let us now briefly consider the coupled FitzHugh-Nagumo oscillators, just in order to stress the properties which are relevant for comparison of the influence of the time-delay on the dynamics of coupled oscillatory vs. excitable FitzHugh-Nagumo systems. Thus the external current $`I0`$ and in the range such that each of the non-coupled units is an oscillator, either close to the Hopf bifurcation or of the relaxation type. We consider the coupling of the same type like in the case of the coupled excitable units (3) and (4). For convenience, the zero of the coupling function is shifted to coincide with the unstable stationary point of the non-interacting oscillators. The major effect of such coupling is to increase the amplitude of each of the oscillators. The amplitude monotonically increases with the coupling constant $`c`$. Furthermore, for the positive coupling constant smaller then some value the asymptotic dynamics of the oscillators is symmetric. However, the oscillations of the two units on the attractor need not be in-phase for larger values of the coupling constant, contrary to the case with oscillations in the instantaneously coupled excitable systems.

Before we pass onto the analysis of the delayed equations, let us mention that the diffusive coupling, when for in the $`\dot{x}_1`$ and $`\dot{x}_2`$ equations one has $`(x_1x_2)`$ and $`(x_2x_1)`$ respectively, also leads to destabilization of the stationary point and appearance of the stable limit cycle. However, in this case $`x_1(t)x_2(t)`$ and the corresponding oscillations are coherent but with a constant phase lag. On the other hand, the trivial stationary point of the system with reversed diffusive coupling is stable for any positive $`c`$, even with an arbitrary time lag.

Delayed coupling:

Let us now consider the dynamics in a neighbourhood of the stationary point of the delayed system (3). The point $`(x_1,y_1,x_2,y_2)=(0,0,0,0)`$ is also the stationary solution of (3), but its stability depends on $`\tau `$. Linearization of the system and substitution $`x_i(t)=A_ie^{\lambda t}`$, $`y_i(t)=B_ie^{\lambda t}`$, $`x_i(t\tau )=A_ie^{\lambda (t\tau )}`$, results in a system of algebraic equations for the constants $`A_i`$ and $`B_i`$. This system has a nontrivial solution if the following is satisfied:

$$\mathrm{\Delta }_1(\lambda )\mathrm{\Delta }_2(\lambda )=0,$$
(10)
where

$`\mathrm{\Delta }_1(\lambda )`$ $`=`$ $`[\lambda ^2+(a+\gamma )\lambda +a\gamma +bc\delta \lambda \mathrm{exp}(\lambda \tau )c\delta \gamma \mathrm{exp}(\lambda \tau )]`$ (11)
$`\mathrm{\Delta }_2(\lambda )`$ $`=`$ $`[\lambda ^2+(a+\gamma )\lambda +a\gamma +b+c\delta \lambda \mathrm{exp}(\lambda \tau )+c\delta \gamma \mathrm{exp}(\lambda \tau )]`$ (12)
The equation (10) is the characteristic equation of the system (3). Infinite dimensionality of the system is reflected in the transcendental character of (10). However, the spectrum of the linearization of the equations (3) is discrete and can be divided into infinite dimensional hyperbolic and finite dimensional non-hyperbolic parts . As in the finite dimensional case, the stability of the stationary point is typically, i.e. in the hyperbolic case, determined by the signs of the real parts of the roots of (10). Exceptional roots, equal to zero or with zero real part, correspond to the finite dimensional center manifold where the qualitative features of the dynamics, such as local stability, depend on the nonlinear terms.

Let us first answer the question of local stability of the stationary point for all time-lags. We have proved (see Appendix 1) that the stationary point remains locally stable for all time-lags if the coupling constant is below some value $`c^\tau `$, which is smaller than $`c_0`$, given by

$$c^\tau =\sqrt{\frac{a^2\gamma ^22b+2\sqrt{b(2\gamma ^2+2a\gamma +b^2)}}{\delta ^2}}<c_0$$
(13)
The previous expression for $`c^\tau `$ is valid if $`b>\gamma ^2`$ which is always satisfied by initial assumptions about the excitable units. Notice that the interval $`(c^\tau ,c_0)`$ is quite small for the units that display excitable behaviour, and shrinks to zero length as $`(b+\gamma )/a0`$.

Thus there could be two qualitatively different types of local dynamics around the stationary solution of the time-delayed coupled excitable ($`c<c_0`$) units. The stationary solution could be a combination of the stable node or the stable focus for $`c<c^\tau `$ and any $`\tau `$, and for $`c^\tau <c<c_0`$ and sufficiently small time-lags, or it could be an unstable focus for $`c>c^\tau `$ and for some time-lags larger than a critical value. The smallest critical time-lag will be found by studying the bifurcations conditions. In the next section, we shall see that there is also an important global bifurcation due to sufficiently large $`\tau `$ inside the interval $`(0,c_0)`$ which changes the dynamics of the excitations.

Bifurcations due to a non-zero time-lag occur when some of the roots of (10) cross the imaginary axes. Let us first discuss the nonzero pure imaginary roots. Substitution $`\lambda =i\omega `$, where $`\omega `$ is real and positive, into the first factor gives

$`c\delta (\omega ^2+\gamma ^2)\mathrm{sin}(\omega \tau )`$ $`=`$ $`\omega ^3+(b\gamma ^2)\omega `$
$`c\delta (\omega ^2+\gamma ^2)\mathrm{cos}(\omega \tau )`$ $`=`$ $`a\omega ^2+(a\gamma +b)\gamma `$
or into the second factor

$`c\delta (\omega ^2+\gamma ^2)\mathrm{sin}(\omega \tau )`$ $`=`$ $`\omega ^3(b\gamma ^2)\omega `$
$`c\delta (\omega ^2+\gamma ^2)\mathrm{cos}(\omega \tau )`$ $`=`$ $`a\omega ^2(a\gamma +b)\gamma `$
Squaring and adding the previous two pairs of equations results in the same equation

$$\omega ^6+(A+\gamma ^2)\omega ^4+(\gamma ^2A+B)\omega ^2+b\gamma ^2=0$$
(14)
where

$$A=a^2+\gamma ^22bc^2\delta ^2,\mathrm{and}B=(a\gamma +b)^2c^2\delta ^2\gamma ^2.$$
(15)
Since $`\omega ^2\gamma ^2`$ the term $`\omega ^2+\gamma ^2`$ can be factored out from (14) to obtain

$$\omega ^4+A\omega ^2+B=0.$$
(16)
Solutions of (16) give the frequencies $`\omega _\pm `$ of possible non-hyperbolic solutions

$$\omega _\pm =\sqrt{(A\pm \sqrt{A^24B})/2}$$
(17)
The corresponding critical time lag follows from equations (13) and (14). Consider the first set (13). Then, if

$$\mathrm{sin}(\omega \tau )=\frac{\omega _\pm ^3+(b\gamma ^2)\omega _\pm }{c\delta (\omega _\pm ^2+\gamma ^2)}>0$$
(18)
we have

$$\tau _{1,\pm }^j=\frac{1}{\omega _\pm }\left[2j\pi +\mathrm{cos}^1\left(\frac{a\omega _\pm ^2+(a\gamma +b)\gamma }{c\delta (\omega _\pm ^2+\gamma ^2)}\right)\right],j=0,1,2\mathrm{}$$
(19)
and if

$$\mathrm{sin}(\omega \tau )=\frac{\omega _\pm ^3+(b\gamma ^2)\omega _\pm }{c\delta (\omega _\pm ^2+\gamma ^2)}<0$$
(20)
we have

$$\tau _{1,\pm }^j=\frac{1}{\omega _\pm }\left[(2j+2)\pi \mathrm{cos}^1\left(\frac{a\omega _\pm ^2+(a\gamma +b)\gamma }{c\delta (\omega _\pm ^2+\gamma ^2)}\right)\right],j=0,1,2\mathrm{}$$
(21)
The analogous critical time-lags from the second factor of the characteristic equation are given as follows. If

$$\mathrm{sin}(\omega \tau )=\frac{\omega _\pm ^3+(b\gamma ^2)\omega _\pm }{c\delta (\omega _\pm ^2+\gamma ^2)}>0$$
(22)
we have

$$\tau _{2,\pm }^j=\frac{1}{\omega _\pm }\left[2j\pi +\mathrm{cos}^1\left(\frac{a\omega _\pm ^2(a\gamma +b)\gamma }{c\delta (\omega _\pm ^2+\gamma ^2)}\right)\right],j=0,1,2\mathrm{}$$
(23)
and if

$$\mathrm{sin}(\omega \tau )=\frac{\omega _\pm ^3(b\gamma ^2)\omega _\pm }{c\delta (\omega _\pm ^2+\gamma ^2)}<0$$
(24)
we have

$$\tau _{2,\pm }^j=\frac{1}{\omega _\pm }\left[(2j+2)\pi \mathrm{cos}^1\left(\frac{a\omega _\pm ^2(a\gamma +b)\gamma }{c\delta (\omega _\pm ^2+\gamma ^2)}\right)\right],j=0,1,2\mathrm{}$$
(25)
The previous formulas give bifurcation curves in the plane $`(\tau ,c)`$ for fixed values of the parameters $`a,b`$ and $`\gamma `$. We denote any bifurcation value of the time-lag by $`\tau _c`$ and add the subscripts and superscripts to specify a particular branch of $`\tau _c(c)`$. The bifurcations are either subcritical Hopf on the curves $`\tau _{1,}^j`$ and $`\tau _{2,}^j`$ leading to a disappearance of one unstable plane, or supercritical Hopf on $`\tau _{1,+}^j`$ and $`\tau _{2,+}^j`$ resulting in appearance of an unstable plane. The type of the bifurcation can be seen by calculation of the variations of the real parts $`\mathrm{𝐑𝐞}\lambda `$ as the time-lag is changed through the critical values. Again, differentiation of the characteristic equation gives

$$\left(\frac{\mathrm{\Delta }_1}{\lambda }\mathrm{\Delta }_2+\mathrm{\Delta }_1\frac{\mathrm{\Delta }_2}{\lambda }\right)\frac{d\lambda }{d\tau }=\left(\frac{Eq_1}{\tau }\mathrm{\Delta }_2+\mathrm{\Delta }_1\frac{\mathrm{\Delta }_2}{\tau }\right)$$

and

$$sgn\left(\frac{d\mathrm{Re}\lambda }{d\tau }\right)_{\tau =\tau _c}=sgn\left\{\mathrm{Re}\left(\frac{d\lambda }{d\tau }\right)^1\right\}_{\tau =\tau _c}=sgn\left(\frac{2\omega ^2+A}{c^2\delta ^2(\omega ^2+\gamma ^2)}\right)$$

Substitution of $`\omega =\omega _\pm `$ finally gives

$$\left(\frac{d\mathrm{Re}\lambda }{d\tau }\right)_{\tau =\tau _+}>0,\left(\frac{d\mathrm{Re}\lambda }{d\tau }\right)_{\tau =\tau _{}}<0.$$
(26)
Let us now discuss the zero solution of (10). Such solution would imply that: $`c=(a\gamma +b)/\gamma \delta `$. For all examples of the coupling functions that we have considered, like: linear, sigmoid, $`\mathrm{tan}^1`$ or $`\mathrm{tanh}`$, this value of $`c`$ was always larger than $`c_1`$, i.e. there were nonzero stable stationary points of (3), So we disregard such solutions of the characteristic equation (10), and concentrate only the Hopf bifurcations $`\lambda =\pm i\omega `$.

The bifurcation curves $`\tau _c(c)`$, given by (19),(21),(23) and (25), are shown in figures 2, 3 and 4, for the first few $`j=0,1,2`$, and for the parameters $`a,\gamma `$ and $`b`$ fixed to some typical values, and for the coupling function with $`f_0^{}=\delta =1`$. Bracketed letters indicate the number of stable and unstable planes in the considered area of the $`(c,\tau )`$ parameter space. For example $`(u^2,u)`$ means two pairs of unstable eigenvalues of the first factor in (10) and one pair of the unstable eigenvalues of the second factor. Analogously, $`(s,s)`$ means that all eigenvalues have negative real parts, i.e. the stationary solution is stable.

Consider first the coupled excitable units when the coupling is in the range $`c(c^\tau ,c_0)`$. In this range, the condition (20) applies for $`\omega =\omega _+`$, and the condition (18) for $`\omega =\omega _{}`$. Accordingly $`\tau _{1,+}^j`$ branches should be calculated with formula (21), and $`\tau _{2,+}^j`$ using (23). The $``$ branches should be computed with (19) for $`\tau _{1,}^j`$ and (25) for $`\tau _{2,}^j`$. This gives the bifurcation curves for $`c(c^\tau ,c_0)`$ presented in figure 2a. The first unstable $`(c,\tau )`$ domain is between the curves $`\tau _{2,+}^0`$ and $`\tau _{2,}^0`$. The unstable plane is given by $`x_1=x_2`$ and $`y_1=y_2`$. The corresponding bifurcation on $`\tau _{2,+}^0`$ is the inverse Hopf bifurcation which results in the destabilization of the stationary point and a collapse of an unstable limit cycle. The origin of the latter is in a global bifurcation which will be discussed in the next section, together with the unique global attractor that exists in this parameter domain. The next unstable region between $`\tau _{1,+}^0`$ and $`\tau _{1,}^0`$, is bordered by a direct supercritical Hopf bifurcation at $`\tau _{1,+}^0`$ and the sub-critical Hopf bifurcation at $`\tau _{1,}^0`$. The unstable invariant manifold of the stationary point is given by $`x_1=x_2,y_1=y_2`$. The stable limit cycle in it supports coherent in-phase oscillations of the two units. The unstable domains bordered by the different branches start to overlap for sufficiently large time-lags, leading to multi-dimensional unstable manifolds of the stationary point. The global attractors for large time-lags are studied in the next section.

Next we consider the range of coupling $`c(c_0,c_1)`$ (see figures 3 and 4). Then, for a sufficiently small $`\tau >0`$, there is only one pair of roots of (11) in the right half-plane, and all the other roots of (11) and (12) have negative real parts. There is an unstable stationary solution and the stable limit cycle. If $`c>c_0`$ but is smaller then some $`c_s`$, the Hopf bifurcation that happens for the smallest time lag corresponds to $`\tau _{1,}^0`$. The value $`c_s`$ corresponds to the intersection of the branches $`\tau _{1,}^0`$ and $`\tau _{2,+}^0`$. If $`c>c_s`$ the first bifurcation occurs for $`\tau _{2,+}^0`$. In the parameter area below the two curves $`\tau _{1,}^0`$ and $`\tau _{2,+}^0`$, denoted by $`(u;s)`$, the stationary solution has qualitatively the same properties as for $`\tau =0`$ i.e. it is unstable and has the unstable 2D manifold with the stable limit cycle in it. The bifurcation at $`\tau _{1,}^0`$ is inverse subcritical Hopf, which results in a stabilization of the stationary solution and in a creation of an unstable limit cycle.

From the set of frames in figure 4, we see that as $`(b+\gamma )0`$ the value $`c_s`$ approaches $`c_0`$ and the $`(s,s)`$ domain beyond $`c_0`$ shrinks to nothing. In fact, in this singular limit, there are only positive solutions of the equation (16), and the stabilization of the stationary point by the time-delay can not happen.

In order to claim that the parameter domain denoted $`(s,s)`$, where the stationary point is stable, corresponds to the phenomenon of oscillators death, the local stability of the stationary point is not sufficient. We need to analyze the global dynamics of (3), and this depends on the full form of the coupling function.

Let us briefly discuss the modifications of the presented analyzes that would be implied by the substitution of the coupling function of the form like in (3) by the diffusive or a more general coupling $`f(x_1,x_2)`$. The analysis of the linear stability in the delayed case, in particular the formulas for the critical time-lags and eigenvalues, would remain unchanged provided that the parameters $`a`$ and $`\delta `$ are changed as follows

$$a\overline{a}=a+_1f_0\delta \overline{\delta }=_2f_0.$$

In particular, for the diffusive coupling, with $`f(x_1,x_2^\tau )=(x_1x_2^\tau )`$, there will be a bounded $`(s,s)`$ region where the stationary point is stable for some finite, nonzero $`\tau `$ and unstable for $`\tau =0`$, and the other parameters unchanged. For example, if $`a=0.25,b=\gamma =0.02,c=0.132`$ the stationary point is stable for $`\tau =0.85`$ up to $`\tau =24.5`$, and unstable for $`\tau =0`$ up to $`\tau =0.85`$ and above $`\tau =24.5`$. However, the unstable manifolds with the limit cycles for small $`\tau `$ would not be given by $`x_1=x_2,y_1=y_2`$ plane.

Contrary to the case of coupled excitable units, the stationary solution of coupled identical FitzHugh-Nagmo oscillators, with the same type of coupling, remains unstable for any value of the time-lag. Thus, there could be no oscillator death in the case of coupled FitzHugh-Nagmo oscillators with the considered type of coupling. On the other hand, it is known (, ) that different type of coupling, like reverse diffusive, does lead to the amplitude death at least when the oscillators are near the Hopf bifurcation, i.e. $`I>I_0;II_0`$. However, as we have pointed earlier, there is no Hopf bifurcation of the trivial stationary point of the excitable systems with the reverse diffusive coupling so that the stationary state is in this case always stable.

## 3  Global dynamics of the system with two units

We study the global dynamics of (3) by numerical computations of orbits for different typical values of the parameters $`(c,\tau )`$ in each of the domains in the local bifurcation diagram of the stationary point up to moderately large values of $`\tau `$. Our main interest is to determine if there is one or more then one attractors, and, if there are stable oscillatory solutions, what is the dimensionality and properties of the oscillations. As examples of the coupling we used different functions, with the same qualitative conclusions, and for illustration we use $`f(x)=\mathrm{tan}^1(x)`$. On the other hand, diffusive coupling $`f(x_1,x_2^\tau )=(x_1x_2^\tau )`$ implies quite different global dynamics, which we shall indicate at the end of the section.

Before presenting the results let us comment on the initial conditions for the system (3) that we used in calculations. In order to uniquely fix a solution of the delay-differential equations (3) one needs to specify the solution on the interval $`[\tau ,0]`$. In our calculations we always use a physically plausible initial functions on $`[\tau ,0]`$ given by solutions of the equations (3) with $`c=0`$. Thus, before the period $`\tau `$ has elapsed each of the units was evolving independently of the other unit.

Firstly we discuss the global dynamics for the coupling constant $`c<c^\tau `$, i.e. when the trivial stationary solution of the whole system is stable for any time-lag. Intuitively, one would expect that if the time-lag $`\tau `$ is such that $`2\tau `$ is smaller than the minimal time needed for an excitable orbit to approach the stable stationary solution, then the coupling would just induce both units to fire one spike each, with some time-delay, and relax to the stationary solution. However, if $`2\tau `$ is larger than the indicated minimal time, then the excitation of one of the units would arrive just in time to kick the state of the other unit from close to the stationary into the excitable domain, even if the coupling constant is rather small. Thus, the excitable orbit of the coupled system becomes periodic. Nevertheless the equilibrium state remains a stable stationary solution. The system as a whole is bi-stable excitable with periodically spiking excitations. This picture is confirmed by numerical computations. Two typical orbits are illustrated in figure 2b. The periodic orbits are supported on the stable limit cycle. The latter is created in a global fold limit cycle bifurcation, together with an unstable one. As expected, the motion on the limit cycle is coherent and out-of-phase, with the frequency that increases with the time-lag.

Qualitatively different global dynamical pictures can occur for the coupling constant in the range $`(c^\tau ,c_0)`$ and for various $`\tau >0`$. In the domain of the $`(c,\tau )`$ plane bordered by $`\tau _{2,+}^0`$ and $`\tau _{2,}^0`$, there is only one global attractor given by the large stable limit cycle, since the value $`\tau _{2,+}^0`$ is above the critical time-lag of the global fold limit cycle bifurcation discussed in the previous paragraph, and the stationary point is unstable (see figure 2c). The large cycle is not affected qualitatively by the local bifurcation on $`\tau _{2,+}^0`$ or $`\tau _{2,}^0`$, so the dynamics on it is given by coherent out of phase oscillations of the coupled units. If the domain is entered by crossing $`\tau _{2,+}^0`$ the unstable limit cycle collapses on the stationary solution, and if the domain is left through $`\tau _{2,}^0`$ yet another unstable limit cycle is created.

The next unstable domain, with $`c`$ still in $`(c^\tau ,c_0)`$ is bounded by $`\tau _{1,+}^0`$ from below and $`\tau _{1,}^0`$ from above. The supercritical Hopf bifurcation on $`\tau _{1,+}^0`$ results in a creation of a 2D unstable manifold of the stationary point given by $`x_1=x_2`$ and $`y_1=y_2`$. In it, there is a stable limit cycle supporting in-phase coherent oscillations. However, the large limit cycle with out-of-phase oscillations is not affected by the local bifurcation at $`\tau _{1,+}^0`$, so that is this domain the dynamics is bi-stable with two limit cycle attractors. This is illustrated in figure 2d.

Next we consider possible attractors for the coupling constant beyond $`c_0`$, i.e. when instantaneously coupled units beheve as two limit cycle oscillators. Qualitatively different dynamics corresponding to different domains in $`(c,\tau )`$, are illustrated in figure 5. Frame 5b corresponds to a typical values of $`(c,\tau )`$ in the domain $`(u,s)`$. We have not found any global bifurcation that would occur as $`(c,\tau )`$ are varied inside the $`(u,s)`$ region. The dynamics is characterized by the 2D unstable manifold of the stationary solution, given by $`x_1=x_2,y_1=y_2`$ (see fig 5a). There is a globally stable limit cycle inside this manifold. Oscillations of $`x_1`$ and $`x_2`$ on this limit cycle are obviously in phase.

Two frames, 5c and 5d, correspond to the situations in $`(s,s)`$ with one stable stationary solution but with two globally different dynamics. In 5c the system is bi-stable. There is the stable stationary solution and the stable large limit cycle in the plane $`x_1=x_2,y_1=y_2`$. There is also a small unstable limit cycle, that is created in the inverse sub-critical Hopf bifurcation at $`\tau =\tau _{1,}^0`$. This cycle acts as a threshold between the sub-excited damped oscillations and periodic synchronous spiking of both units. As $`\tau `$ is increased, but still for $`(c,\tau )(s,s)`$, the unstable limit cycle approaches the stable one, and they disappear in a fold limit cycle bifurcation, which occurs in the invariant plane $`x_1=x_2,y_1=y_2`$. Thus, there is a parameter region inside $`(s,s)`$ where $`(0,0,0,0)`$ is globally stable, and that corresponds to the death of the identical oscillators. However, let us stress again that the global dynamics for the parameters in the domain $`(s,s)`$ could correspond to either spiking excitability, with sub-threshold dumped oscillations and sup-threshold periodic spiking, or to the death of oscillators. In the latter regime the whole system is excitable with the stationary point as the only attractor.

The global dynamics above the curve $`\tau _{2,+}^0`$ is characterized by one large limit cycle as the global attractor. The oscillations on it are coherent and out-of-phase. The same type of the global attractor occur above the critical line $`\tau _{1,+}^1`$ as illustrated in frames 5e and 5f. The oscillations are further illustrated in figure 6b by plotting the limit cycle as seen in the coordinates $`(x_1x_2,y_1y_2)`$. Convergence to the limit cycle is much slower than in the case of the symmetric oscillators that occurs for smaller $`\tau `$, illustrated in figure 6a. In fact, in all domains up to a quite large value of the time-lag the global attractor is a stable limit cycle (could be imbedded in a multi-dimensional unstable manifold for larger $`\tau `$) which supports asymmetric phase shifted oscillations of $`x_1`$ and $`x_2`$. However, there are domain for larger values of the time-lag, for example for $`\tau =55`$ and any $`c(c_0,c_1)`$, where the global attractor is the symmetric limit cycle, with coherent and in-phase oscillations.

It should be pointed out that, for all larger time-lags up to quite large values, equal to several refractory times of the non-coupled units, the attractor is always a limit cycle. On the limit cycle, all variables oscillate with the same frequency, and could be either symmetric or phase shifted. The two regimes interchange as the time-lag is increased. These are the only two possible stable attracting patterns, despite the large dimensionality of the unstable manifold of the stationary point. It should be pointed out that in the case of identical FitzHugh-Nagumo oscillators with the same coupling various types of quasi-periodic attractors occur for moderate values of the time-lag. However, also in this case, stronger coupling and larger time-lags imply synchronization, either identical or phase shifted, like for the coupled excitable systems. The dynamics for time-lags much larger than the refractory time has not been systematically studied.

We now briefly comment on the dynamics in the case of the diffusive coupling. As stated before, at some $`c`$ and for $`\tau `$ zero or small, the only attractor is the stable limit cycle, with coherent and phase-shifted oscillations of the two units. The time delay can stabilize the trivial stationary point, but the system remains bi-stable with the limit cycle and the stationary point as the attractors, for all values of $`(c,\tau )`$ in $`(s,s)`$ domain. The phenomenon of oscillator death does not occur in the case of the diffusive coupling.

## 4  $`N>2`$ lattice

The goal of this section is to present numerical evidence that for some common types of lattices with $`N>2`$ there are regions in the parameter plane $`(c,\tau )`$ analogous to $`(u,s),(s,s),(u,u)\mathrm{}`$ in figures 2 and 3. We have analyzed examples of systems of identical FitzHugh-Nagumo excitable units arranged in linear or circular lattices, with unidirectional or bi-directional symmetrical coupling by few typical coupling functions. Lattices of the size $`N=10,20`$ and $`N=30`$ have been studied systematically.

The conclusions are illustrated using the following model:

$`\dot{x}_i`$ $`=`$ $`x_i^3+(a+1)x_i^2ax_iy_i+cf(x_{i1}^\tau )+cf(x_{i+1}^\tau ),`$
$`\dot{y}_i`$ $`=`$ $`bx_i\gamma y_i,i=2,\mathrm{}N1,`$ (27)
$`\dot{x}_{1,N}`$ $`=`$ $`x_{1,N}^3+(a+1)x_{1,N}^2ax_{1,N}y_{1,N}+cf(x_{2,N1}^\tau ),`$
$`\dot{y}_{1,N}`$ $`=`$ $`bx_{1,N}\gamma y_{1,N},`$ (28)
where the coupling is given by $`f(x)=\mathrm{tan}^1(x)`$, and $`N=20`$.

Firstly, in the case of instantaneously coupled units, there is the Hopf bifurcation at some $`c=c_0`$. As in $`N=2`$, for the coupling constant below some $`c^\tau <c_0`$ and any time-lag the trivial stationary point is stable. If the time-lag is sufficiently large, there is also the stable large limit cycle. On it, all units oscillate coherently. However, the nearest neighbours oscillate exactly anti-phase, so that two clusters are formed.

The coupling above the threshold $`c_0`$, and for small time lags, leads to the appearance of a globally stable limit cycle representing synchronous oscillations in the plane given by $`x_2=\mathrm{}=x_{N1},y_2=\mathrm{}=y_{N1}`$ and $`x_1=x_N,y_2=y_N`$, ( $`(u,s)`$ region). As expected, the synchronization period could be quite large if the value of the coupling constant is near $`c_0`$, i.e. when each of the units is near the Hopf bifurcation.

Increasing the time lag leads to the inverse Hopf bifurcation. For any $`c`$ in some interval $`(c_0,c_s)`$ we have been able to find intervals of time-lags $`(\tau _c^{},\tau _c^+)`$ that correspond to bi-stability or to death of all $`N`$ oscillators ($`(s,s)`$ region), illustrated in figure 7. The same figures illustrate the attractors in the dynamics of any of the identical neurons. Again, the inverse Hopf and the subsequent fold limit cycle bifurcations, due to increasing time-delay, are responsible for the amplitude death in the systems of the form (27). On the contrary, the stationary point of a lattice like (27) with the same coupling but with Fichugh-Nagumo oscillators is always unstable.

Larger time lags do not change the topological nature of the attractor. It is always a limit cycle, but the synchronization pattern between the coherent oscillations of the units depends on $`\tau `$. Non-symmetric oscillations with equal frequencies appear after long transients, as is illustrated in figure 8a. Dynamics in the transient period can be quite complicated. Properties of the asymptotic synchronization patterns could depend on the geometry of the lattice.

Existence of the presented types of dynamics, and the order of their appearance as $`(c,\tau )`$ are varied, was confirmed in all examples that we have studied. We conjecture that qualitative properties of the dynamics of all small 1D lattices with nearest-neighbor delayed coupling of the form like in the equations (27) between the identical FitzHugh-Nagumo excitable systems are the same, at least for not very large values of the time lag, in the sense that the same types of bifurcations appear and determine the dynamics.

## 5  Summary and discussion

We have studied small lattices of excitable identical units with time-delayed coupling, where each unit is given by the excitable FitzHugh-Nagumo model. The coupling is always between the voltages of the nearest neighbors, but could be of a quite general form. Our primary interest was in the bifurcations and the typical dynamics that occur for time-lags which are not very large on the time-scale set up by the refractory or the inter-spike period. Detailed study, in the case of only two units, of the local stability and bifurcations of the stationary solution suggested, but does not uniquely determine, the possible global bifurcations and dynamics. These are studied numerically.

There are only few possible types of dynamics, at least for time lags as large as several refractory times. For small coupling constants and small time-lags there is only one attractor in the form of the stable stationary solution. The whole system beehives qualitatively as the simple excitable. Relatively small coupling constants $`c<c^\tau `$ and sufficiently large time-lag result in the limit cycle attractor co-existent with the stable stationary solution. The whole system is bi-stable with spiking excitability. The oscillations on the limit cycle are coherent and out-of-phase. For the coupling constant above $`c^\tau `$ the sequence of Hopf bifurcations due to the time-delay of the stationary solution are possible. For $`c(c^\tau ,c_0)`$ and small time lags the stable stationary point is the only attractor, but as time-lag is increased the system could be either bi-stable or there could be only one attractor in the form of the limit cycle. The bi-stability is manifested either in the form of the stable stationary solution and the stable limit cycle, or could be in the form of two stable limit cycles (one in-phase and one out-of-phase). The interval $`c(c^\tau ,c_0)`$ is rather a small part of the $`c`$ values for which there is only one stationary solution. For $`cc_0`$ each of the instantaneously coupled units is near a direct super-critical Hopf bifurcation, but as soon as $`cc_0>0`$ is bigger then some quite small $`ϵ_0`$ the resulting limit cycle has quite large radius and the harmonics become influential, unlike in the case of the Hopf limit cycle. Increasing the time-lag $`\tau `$ could lead to stabilization of the stationary point, via indirect sub-critical Hopf, resulting in a bi-stable dynamics, with a stable stationary point, small unstable limit cycle as a threshold, and a large stable limit cycle. Further increasing $`\tau `$ leads to a fold limit cycle bifurcation, in which the unstable and the stable limit cycles disappear, and the stationary point remains the only attractor. In this, oscillator death regime, the system again displays the simplest form of excitability, like in the case of the weak coupling $`c<c_0`$ and zero or small time lags. Still further increase of the time-lag $`\tau `$ leads to the super-critical Hopf. The oscillations on the limit cycle are coherent but are phase shifted, and the oscillators need not have the same amplitude. Described sequence of bifurcations happens for time lags that are all small, up to $`10\%`$ , with respect to the refractory period of the single isolated unit. Further increase of the time-lag leads to more dimensional unstable manifold of the stationary solution. However, the global attractor is always a simple limit cycle. The asymptotic dynamics is always coherent, and is either in-phase or phase shifted. Unlike the case of coupled FitzHugh-Nagumo oscillators, nothing more complicated than the limit cycle could be the attractor of the coupled excitable FitzHugh-Nagumo systems.

Our analyzes shows that the most common type of excitations of the whole system, in response to an impulse submitted to either of the units, is in the form of coherent out-of-phase oscillations. However, if the transmission is sufficiently strong and for moderately large transmission delays of signals between the units, the compound system would respond by synchronous in-phase oscillations. Furthermore, our results suggest that relatively small but non-zero time-delay together with sufficiently strong interaction could result in a simple excitable behaviour of the compound system. For such values of the parameters the system would operate as a powerful amplifier of a quite small impulse administered to its single unit. Due to the particular model of the excitable system and to the type of coupling that we have studied in detail, the most relevant possible application of our results is in modeling coupled neurons. In fact, relatively recent experiments and analyzes  show that the FitzHugh-Nagumo equations, despite the common opinion, might represent better qualitative model of an excitable neuron than the more detailed Hodgkin-Huxly system. Our results indicate that the fine tuning between the synaptic coupling and delay could lead to the in-phase synchronous operation of a collection of neurons.

Although there is a quite substantial amount of work done on the systems of dynamical units with the delayed coupling, such systems are comparatively much less studied than the corresponding systems with the instantaneous coupling. For the purpose of comparison with our work, we shall try to classify the existing contributions into typical groups.

Firstly, we consider the model and the results presented in ,. In these papers, a network of $`N=2`$ and $`N>2`$ oscillators described by the equations of the normal form of the Hopf bifurcations with delayed inverse diffusive coupling is studied. At zero coupling, and/or for small time lags, all oscillators have small limit cycles just created by the direct Hopf bifurcation. It is shown that the time-delay can lead to the stabilization of the trivial stationary point, even for the identical oscillators, which is interpreted as the amplitude death. Analogously, in our case, the Hopf direct bifurcation is responsible for the appearance of the oscillations when the excitable units are instantaneously coupled, and the time-delay leads to stabilization of the stationary point. The death of oscillations in our case appears after the fold bifurcation of the stable and the unstable limit cycles, which are created in the same plane. The oscillator death occurs only in a domain in the $`(c,\tau )`$ parameter space smaller that the domain of the stability of the stationary point.

Next, we compare our model and the results with those that appear in the studies of the delayed coupled relaxation oscillators, for example in : ,. In these studies, each unit is a relaxation oscillator, and the primary objective of the analysis are the periodic orbits that appear in the delayed coupled system. Singular approximation, or an approximate or numerical construction of the Poincare map, are used to analyze various types of synchronous or asynchronous oscillations. The phenomenon of the oscillator death was not observed ( ). In our case, the noniteracting units are not oscillators and the oscillations are introduced by coupling, via the Hopf bifurcation. The domain of parameters $`(c,\tau )`$ that implies oscillator death shrinks to nothing in the singular limit $`(b+\gamma )0`$. Furthermore, the FitzHugh-Nagumo model is type II excitable, which reflects in the type of bifurcations that might occur in the coupled systems.

Less directly related to our work is the analyzes of the influence of the time delay in the systems of coupled phase oscillators. In fact, if the rate of attraction to the limit cycles of two voltage-coupled neural oscillators is sufficiently strong the dynamics can be described by the coupled phase oscillators. The coupling between the phases mimics the voltage coupling, and is not of the diffusive type. The phenomenon of oscillator death in such instantaneously coupled phase oscillators was studied for example in . The influence of time-delay in coupled phase oscillators was studied for example in  and  (and also in ), where it is shown that the time-delay can not introduce significant changes into the dynamics of a class of such systems , unless the time-lag is of the order of several oscillation periods . Independently of neuronal models, collective behavior of the phase coupled (phase) oscillators with time-delayed coupling, have been studied using the dynamical ,,,), or statistical  methods. In our case, the coupling is between the voltages, could be of a quite general form, and all analyzes and the observed phenomena occur already for quite small time lags.

Finally, the influence of time delay has been studied in the Cohen-Grasberger-Hopfield (CGH) type of neural networks, as early as in 1967 . More recent references are for example , , and for small networks ,, , (see also  and the references there in). In the non-delayed case, the stability of the stationary point in such networks is proven using an energy-Lyapunov function. Using the corresponding Lyapunov functional in the delayed case, it was shown in , that the stationary point remains globally stable for sufficiently small time lags. On the other hand, destabilization of the stationary point occurs via the Hopf bifurcation, as was shown in ,, for the networks with $`N=2`$ and $`N=3`$ units, and multiple time-delayed coupling.

As is seen, the model treated here, and our results, have some features common with few other models. As in the CGH networks, each isolated unit has a globally stable stationary point, and the time-delayed weak coupling does nothing to the dynamics, provided that the time-lag is sufficiently small. If the coupling is strong enough, the system behaves either as a collection of near-Hopf oscillators, or as a collection of relaxation oscillators. Death of oscillators due to time-delay is observed in both types of dynamics, although the phenomenon happens for a smaller range of time-lags if the system behaves as a collection of relaxation oscillators.

Let us finally mention few related questions that we shall study in the future. Firstly, it should be interesting to see if the systems of slightly different units share the same type of dynamics. Secondly, examples of type I excitable (and not oscillatory) systems coupled with time-delay should be analyzed, in order to underline the role of the type of excitability. Finally, the external pulse perturbations, like for example in , ,, could introduce different transitions from excitability to the oscillatory regime, and the consequent changes in the dynamics due to time-delay should be analyzed.

Acknowledgements This work is supported by the Serbian Ministry of Science contract No. 1443.

## 6  Appendix

We start with the characteristic function (10) in the form:

$$\varphi (z)=\left[(z+\gamma )(z+a)+b\right]^2c^2\delta ^2(z+\gamma )^2\mathrm{exp}(2z\tau ),$$

and consider the following expression:

$$f(z)=\frac{\varphi (z)}{P_4(z)}=1\frac{c^2\delta ^2(z+\gamma )^2}{P_4(z)}\mathrm{exp}(2z\tau ),$$

where $`P_4(z)=\left[(z+\gamma )(z+c)+b\right]^2`$ is actually the characteristic function of the single non-coupled unit.

Consider the contour $`C_R`$ in the complex half-plane $`\mathrm{𝐑𝐞}z>0`$, formed by the segment $`[iR,iR]`$ of the imaginary axis and the semi-circle with the radius $`R`$ centered at the origin. As the condition $`4(b/\gamma )>(a1)^2`$, equivalent to the existence of a unique stationary solution, is by assumption always satisfied, the polynomial $`P_4(z)`$ has no zeros in the half-plane $`\mathrm{𝐑𝐞}z>0`$. In that case, the number of poles of $`f(z)`$ is $`P_c=0`$. Using the argument principle we infer the number of zeros $`N_{C_R}`$ of $`f(z)`$. If $`lim_R\mathrm{}N_{C_R}=0`$ then all the roots of the characteristic function $`\varphi (z)`$ satisfy $`\mathrm{𝐑𝐞}z>0`$. Thus, we need to find the conditions on the parameters $`a,b,\gamma `$ and $`c`$, such that the image of the contour $`C_R`$ when $`R\mathrm{}`$ by the function $`f(z)`$ does not encircle the point $`z=0`$. Then the variation of the argument is zero, so that $`lim_R\mathrm{}N_{C_R}=0`$, and consequently the zeros of the characteristic function satisfy $`\mathrm{𝐑𝐞}z>0`$ for any $`\tau `$. This is the essence of the amplitude-phase method ( see for example .)

It is enough to consider the image of the segment $`[iR,iR]`$ by the function

$$\omega _\tau (z)\frac{c^2\delta ^2(z+\gamma )^2}{P_4(z)}\mathrm{exp}(2z\tau ),$$

or, in fact, by just $`\omega _0(z)`$ since $`|\omega _0(iy)|<1`$ if and only if $`|\omega _\tau (iy)|<1`$, and the image of the semi-circle shrinks to a point as $`R\mathrm{}`$.

Since

$$|\omega _0(iy)|=|\frac{c\delta (iy+\gamma )}{(iy+\gamma )(iy+a)+b}|^2=\frac{c^2\delta ^2(\gamma ^2+y^2)}{y^4+(a^2+\gamma ^22b)y^2+(a\gamma +b)^2}$$

we obtain that $`|\omega _0(z)|<1`$ is equivalent with

$$y^4+Ay^2+B>0$$

where $`A`$ and $`B`$ are given by the same formula as in (15), i.e.

$$A=a^2+\gamma ^22bc^2\delta ^2,\mathrm{and}B=(a\gamma +b)^2c^2\delta ^2\gamma ^2.$$

For the parameters such that $`b>\gamma ^2`$ and $`4(b/\gamma )>(a1)^2`$ the above condition is equivalent to

$$c<\sqrt{\frac{a^2\gamma ^22b+2\sqrt{b(2\gamma ^2+2a\gamma +b^2)}}{\delta ^2}}.$$

The right side is the critical value that we denoted $`c^\tau `$ in the main text.

FIGURE CAPTIONS

1)Figure 1 The figure illustrates continuous transition of the limit cycles from near Hopf to that of a relaxation oscillator for the coupled system with $`\tau =0`$. The fixed parameters are: $`a=0.25,b=\gamma =0.02`$. The cycle is created at $`c=0.27`$ and the smallest cycle on the figure is for $`c=0.2702`$, the next to the largest for $`c=0.27048`$ and the largest for $`0.3`$.

2) Figure 2 The figures illustrate typical dynamics below $`c_0=0.27`$ for $`a=0.25,b=\gamma =0.02`$: a) First few branches of the bifurcation curves $`\tau _c(c)`$ given by equations (19),(21),(23) and (24), for the parameters $`a=0.25,b=0.02,\gamma =0.02`$ and $`c<c_0`$; b) Examples of quickly relaxing (1) and periodic excitable (2) orbits (projections on $`(x_1,y_1)`$; (c) Projection of the global attractor limit cycle on $`(x_1,x_2)`$ in the domain $`(s,u)`$; (d) Projection of the two limit cycle attractors on $`(x_1,x_2)`$ in the domain $`(u,u)`$.

3) Figure 3: First few branches of the bifurcation curves $`\tau _c(c)`$ given by equations (19),(21),(23) and (24), for the parameters $`a=0.25,b=0.02,\gamma =0.02`$ and $`c>c_0`$.

4) Figure 4: The same as figure 3, but for few fixed values of $`a,b,\gamma `$: a)$`a=0.25,b=0.005,\gamma =0.005`$; b)$`a=0.25,b=0.003,\gamma =0.003`$; c)$`a=0.25,b=0.0015,\gamma =0.0015`$; d)$`a=0.25,b=0.00075,\gamma =0.00075`$

5)Figure 5: Phase portraits in $`(x_1,y_1)`$ (a,b,c,d,e) plane, or $`(x_1,x_2)`$ plane (f). The initial points if there are different orbits are indicated by numbers. The fixed parameters are $`a=0.25,b=0.02,\gamma =0.02,c=0.3`$, except in (a) where $`c=0`$. The time-lag is: (a),(b) $`\tau =0`$; (c) $`\tau =4`$; (d)$`\tau =6`$; (e),(f) $`\tau =27`$. Dynamics illustrated in (e) and (f) is typical also for other values of $`(c,\tau )`$ above $`\tau _{1,+}^1`$ curve.

6) Figure 6: (a) The asymptotic state is symmetric in any domain below $`\tau _{1,+}^1`$ curve in figure 3. In domains $`(u,s)`$, $`(s,u)`$ or $`(u,u)`$ the symmetric state are the synchronous oscillations, and in $`(s,s)`$ the stable stationary point. (b) The asymptotic synchronous oscillations are not symmetric for $`(c,\tau )`$ above the curve $`\tau _{1,+}^1`$ in figure 3, as is illustrated for a pair $`(c,\tau )(u,u)`$, but become symmetric for $`\tau 55`$ (not illustrated, see the main text).

7) Figure 7: Bi-stability (a) and oscillator death in the lattice (27) with $`a=0.25,b=0.02,\gamma =0.02`$ and (a) $`(c,\tau )=(0.16,4)`$ or (b) $`(c,\tau )=(0.16,6)`$.

8) Figure 8: Asymptotic states of the lattice (27) for $`(c,\tau )=(0.16,15)`$ are coherent oscillations but with a fixed time lag, represented by the projection of the limit cycle on the $`(x_4x_{15},y_4y_{15})`$ plane in the frame (b). For such a small $`c`$ the synchronization period is more then 10 times larger than the characteristic period, as is illustrated in frame (b) with the time dependence of the time-series $`x_4(t)`$ and $`x_{15}(t)`$


format: text

text:
# LOCAL QUASI-EQUILIBRIUM DESCRIPTION OF SLOW RELAXATION SYSTEMS

## I Introduction

Systems whose global evolution results from the existence of a wide variety of time and length scales frequently exhibit slow relaxation dynamics, which is manifested through the aging behaviour of the correlations and the nonexponential decay of the response function angell -stillinger . To characterize the dynamics of these systems is one of the most challenging problems of nonequilibrium statistical mechanics nowadays. Many efforts have been made in this sense by proposing different approaches including spin models, stochastic dynamics and thermodynamical approaches (see, for example, Refs. Anderson -theo ).

Experiments performed in diverse systems as in amorphous polymers struick  and supercooled liquids and glasses angell , israeloff -coll-glasses , seem to indicate that these systems undergo, in general, non-Markovian and non-stationary dynamics leading to the peculiar behaviour of the two-time correlation functions and to the dependence of the relaxation times on the initial time.

In this paper we present a dynamical description of slow relaxation systems based on the formulation of a Fokker-Planck equation for non-Markovian processes. This equation is obtained through the generalization of the Onsager theory, originally proposed to account for the dynamics of fluctuations around equilibrium states onsager  and extended to nonequilibrium stationary states dorfman , to nonequilibrium *aging* states.

An important feature of the dynamics of these systems is that at the time scales of experiments, they still relax and never reach thermal equilibrium with the bath. We have assumed the existence of a local quasi-equilibrium state virasoro  characterized by a non-stationary probability density which can be used to determine average quantities as, for example, the local quasi-equilibrium temperature librozwanzig . The entropy of the system has been expressed in terms of the conditional probability density by means of the Gibbs entropy postulate, and used to derive the generalized Fokker-Planck equation by taking into account probability conservation and the rules of nonequilibrium thermodynamics degroot . For stationary processes, this equation was derived in adelman  and pop-mori , and has subsequently been treated in the literature of stochastic processes in oxtoby -FP-nonmarkovian .

The unified description of slow relaxation systems we propose can be applied to, among others, glassy systems kurchan , granular matter dufty1  and anomalous diffusion problems tokuyama , oliveira .

The paper is distributed as follows. In Sec. II, we will analyze the non-stationary dynamics of slow relaxation systems under the framework of a generalized Onsager theory. Section III is devoted to the formulation of relaxation equations and of the generalized fluctuation-dissipation theorem for non-stationary systems. An illustration of the theory presented is given in Sec. IV by analyzing the dynamics of systems having two well-separated time scales (glassy systems). In Sec. V, we will discuss the Brownian motion in a granular gas. We will compare our results with those obtained by means of kinetic theory and simulations. In the discussion section we will summarize our main results.

## II Fokker-Planck equation for non-stationary dynamics

The non-stationary dynamics of slow relaxation systems can be characterized at the mesoscopic level by means of the conditional probability density $`P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)`$, which satisfies the normalization condition $`P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)𝑑\underset{¯}{\alpha }=1`$. This quantity depends in general on a set of fluctuating variables $`\underset{¯}{\alpha }=(\alpha _1,\mathrm{},\alpha _n)`$, or fields $`\underset{¯}{\alpha }(\stackrel{}{k})`$, hydrofluc , necessary to univocally determine the state of the system during the relaxation and which may represent, for instance, fluctuating positions and velocities, order parameters, reaction coordinates or hydrodynamic fields.

The evolution in time of $`P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)`$, with $`\underset{¯}{\alpha }_0`$ the state of the system when the fluctuation takes place at $`t_0`$ (in glassy systems literature it is referred to as waiting time or time elapsed after quenching), is governed by the continuity equation

$$\frac{}{t}P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)=\frac{}{\underset{¯}{\alpha }}\left[P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)\underset{¯}{v}_{\underset{¯}{\alpha }}(\underset{¯}{\alpha }_0,\underset{¯}{\alpha };t_0,t)\right],$$
(1)
expressing probability conservation during its natural evolution in $`\underset{¯}{\alpha }`$-space. The probability current $`P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)\underset{¯}{v}_{\underset{¯}{\alpha }}(\underset{¯}{\alpha }_0,\underset{¯}{\alpha };t_0,t)`$, with $`\underset{¯}{v}_{\underset{¯}{\alpha }}(\underset{¯}{\alpha }_0,\underset{¯}{\alpha };t_0,t)`$ a stream velocity in $`\underset{¯}{\alpha }`$-space, has to be determined in order to render the dynamical description complete.

This task can be accomplished by assuming that the scheme of nonequilibrium thermodynamics degroot  remains valid in $`\underset{¯}{\alpha }`$-space. Thus, a linear relationship can be established between the diffusion current $`P(t_0|t)\underset{¯}{v}_{\underset{¯}{\alpha }}(t_0,t)`$ and its conjugate force in $`\underset{¯}{\alpha }`$-space, which can be obtained from the entropy production calculated in the framework of Onsager’s theory degroot . To derive the entropy production we will assume that, for sufficiently long times $`t`$, the system reaches a local quasi-equilibrium state, characterized by the entropy per unit mass $`s_{qe}(t)`$ and the probability density $`P_{qe}(\underset{¯}{\alpha },t)`$, related formally to the conditional probability through $`lim_t\mathrm{}P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)=P_{qe}(\underset{¯}{\alpha },t)`$.

In accordance with inertial , the variation of the entropy $`s_{qe}(t)`$ satisfies the Gibbs equation

$$T(t)\delta s_{qe}(t)=\delta e_{qe}(t)\mu _{qe}(t)\delta P(t)𝑑\underset{¯}{\alpha },$$
(2)
where the temperature $`T(t)`$ of the system at local quasi-equilibrium will be assumed to be only a function of time, $`e_{qe}(t)`$ is the mean internal energy, and $`\mu _{qe}(t)`$ the chemical potential at local quasi-equilibrium, which do not depend on the $`\underset{¯}{\alpha }`$’s. For simplicity we have omitted the explicit dependence of the probability density on phase variables. The entropy production $`\sigma `$ corresponding to the diffusion process in $`\underset{¯}{\alpha }`$-space can be computed by using the Gibbs entropy postulate kampenslow

$$s(t)=\frac{k_B}{m}P(t_0|t)\mathrm{ln}\frac{P(t_0|t)}{P_{qe}(t)}d\underset{¯}{\alpha }+s_{qe}(t),$$
(3)
where $`s(t)`$ is the entropy of the system, $`k_B`$ the Boltzmann’s constant and $`m`$ the molecular mass. Notice that, according with mesoscopic nonequilibrium thermodynamics (MNET), the entropy $`s_{qe}(t)`$ appearing in Eq. (3), constitutes a macroscopic quantity whose variations are governed by Eq. (2). The evolution in time of $`s_{qe}(t)`$ can also be described with the MNET scheme by using the expression of Eq. (2), and the corresponding Eq. (3) referred to the equilibrium state (see, for example, inertial  and MNETmazur ). In that case, by taking the variation $`\delta s_{qe}=s_{qe}(\left\{P_{qe}+\delta P\right\})s_{qe}(\left\{P_{qe}\right\})`$ at first order in $`\delta P`$ and combining the resulting expression with Eq. (2), we obtain

$$P_{qe}(t)=P_ee^{\frac{m}{k_BT(t)}[\mu _{qe}e_{qe}]},$$
(4)
which incorporates the topology of the energy landscape through $`e_{qe}`$, and where $`P_e`$ is the equilibrium probability distribution.

By combining the rate of change of $`s_{qe}(t)`$ obtained from Eq. (2) with the time derivative of Eq. (3), after using Eq. (1) and integrating by parts assuming that $`P(t_0|t)\underset{¯}{v}_{\underset{¯}{\alpha }}(t_0,t)`$ vanishes at the boundary, we obtain gradtemp

$$\frac{}{t}s(t)=\frac{1}{T(t)}\frac{}{t}e(t)\frac{1}{T(t)}P(t_0|t)\underset{¯}{v}_{\underset{¯}{\alpha }}(t_0;t)\frac{}{\underset{¯}{\alpha }}\left[\frac{k_BT(t)}{m}\mathrm{ln}\frac{P(t_0|t)}{P_{qe}(t)}+\mu _{qe}(t)\right]𝑑\underset{¯}{\alpha },$$
(5)
where $`\delta e(t)=\delta e_{qe}(t)+\frac{k_BT(t)}{m}P(t_0|t)\delta \mathrm{ln}P_{qe}(t)𝑑\underset{¯}{\alpha }`$ is the variation of the nonequilibrium internal energy $`e(t)`$. If no external thermodynamic forces as, for example, velocity or temperature gradients act on the system, it may be shown by using the Gibbs-Duhem relation that $`\delta e(t)=0`$. As a consequence of this, the entropy production is only due to the underlying diffusion process in $`\underset{¯}{\alpha }`$-space. Moreover, in the last term of Eq. (5) one may identify the nonequilibrium “chemical potential” $`\mu (t_0;t)`$ as

$$\mu (t_0;t)=\frac{k_BT(t)}{m}\mathrm{ln}\frac{P(t_0|t)}{P_{qe}(t)}+\mu _{qe}(t).$$
(6)
Important to notice is the fact that $`\mu (t_0;t)`$ may depend, in general, on interaction potentials or entropic barriers Entropic barriers  through the corresponding dependence of $`P_{qe}(t)`$.

According to the rules of nonequilibrium thermodynamics expressing currents as linear functions of the conjugate thermodynamic forces degroot , from the entropy production (5) the velocity $`\underset{¯}{v}_{\underset{¯}{\alpha }}(t_0;t)`$ can now be expressed in terms of $`\frac{}{\underset{¯}{\alpha }}\mu (t)`$, obtaining

$$\underset{¯}{v}_{\underset{¯}{\alpha }}(t_0;t)=\underset{¯}{\underset{¯}{B}}(t_0;t)\left[\underset{¯}{X}(\underset{¯}{\alpha },t)+\frac{k_BT(t)}{m}\frac{}{\underset{¯}{\alpha }}\mathrm{ln}P(t_0|t)\right],$$
(7)
where we have identified the generalized force FP-nonmarkovian , degroot

$$\underset{¯}{X}(\underset{¯}{\alpha },t)=\frac{}{\underset{¯}{\alpha }}\left[\frac{k_BT(t)}{m}\mathrm{ln}P_{qe}(t)\right],$$
(8)
and introduced the coefficients $`\underset{¯}{\underset{¯}{B}}(t_0;t)`$ which may, in general, depend on the state variables $`\underset{¯}{\alpha }`$. These coefficients, incorporating memory effects through their time dependence, can be determined in terms of the conditional relaxation functions of the corresponding variables $`\underset{¯}{\alpha }`$, degroot . Notice that, in general, the $`B_{ij}`$’s do not satisfy reciprocal relations since for non-stationary systems and non-conservative interactions, the time reversal symmetry and the principle of detailed balance are not necessarily fulfilled. The time dependence of the $`B_{ij}`$’s is the origin of anomalous diffusion, which is characterized by the power law behaviour in time of the correlation function tokuyama , oliveira , and the nonexponential time decay of the relaxation functions angell .

By substituting the expression of the probability current obtained from Eq. (7) into Eq. (1), we obtain the Fokker-Planck equation

$$\frac{P(t_0|t)}{t}=\frac{}{\underset{¯}{\alpha }}\underset{¯}{\underset{¯}{B}}(t_0;t)\left[\underset{¯}{X}(\underset{¯}{\alpha },t)P(t_0|t)+\frac{k_BT(t)}{m}\frac{P(t_0|t)}{\underset{¯}{\alpha }}\right].$$
(9)
Similar equations describing stationary non-Markovian dynamics were obtained in Refs. adelman -FP-nonmarkovian , to study linear and nonlinear transport. From Eq. (9), it follows that the characteristic relaxation times of the system $`B_{ij}^1(t_0;t)`$, are functions of the waiting time $`t_0`$ as a consequence of the non-Markovian nature of the dynamics.

At the time scales considered, the system relaxes towards the time-dependent local quasi-equilibrium state and then thermal equilibrium between system and bath can not be established. The complete description of the dynamics can then be performed by giving the relation between the system and bath temperatures. In order to derive it, the influence of the bath is incorporated through the second law $`\delta s_T=\delta s_{qe}+\delta s_B0`$, with $`s_T`$ the total entropy of the system and $`s_B`$ that of the bath. From this expression, by assuming that the variation $`\delta s_T`$ takes place at constant volume (see Eq. (2)), it can be expressed as $`\delta s_T=[1c]\frac{1}{T(t)}[\delta e_{qe}(t)\delta g_{qe}(t)]`$, where we have taken into account the Gibbs equation (2), used the Gibbs free energy per mass unit $`\delta g_{qe}(t)=\mu _{qe}(t)\delta P(t)𝑑\underset{¯}{\alpha }`$ and defined $`c=\frac{\sigma _B}{\sigma }1`$, with $`\sigma _B`$ the entropy production of the bath. The equality occurs when the process is reversible, i.e., $`\sigma _B=\sigma `$. In a similar way, by using the Gibbs equation of the bath we obtain $`\delta s_T=[1c^1]\frac{1}{T_B(t)}\delta e_B(t)`$, with $`e_B(t)`$ and $`T_B(t)`$ the energy and temperature of the bath. A simple comparison of the expressions for $`\delta s_T`$ leads to the relation

$$T(t)=AT_B(t),$$
(10)
where we have defined $`A=c\left[\frac{\delta g_{qe}(t)\delta e_{qe}(t)}{\delta e_{B\text{ }}}\right]`$. Since $`\delta g_{qe}(t)=s_{qe}\delta T(t)\mu _{qe}\delta \rho `$, we conclude that the factor $`A`$ also depends on the density and the temperature of the system. When the process is reversible, $`A=1`$ since thermal equilibrium between system and bath is attained. In this case, $`P_{qe}(t)`$ reduces to a local equilibrium probability density in which the temperature is that of the bath Entropic barriers , landau .

## III Relaxation equations and the fluctuation- dissipation theorem

Once the Fokker-Planck equation has been derived, we will proceed to establish the relaxation equations for the moments of the probability density. The analysis of the general case can be accomplished by considering that the coefficients are functions of the state variables FP-nonmarkovian . For the sake of simplicity, we will focus on the case in which the coefficients $`B_{ij}`$ are only functions of time.

First, we will consider the two-time correlation function

$$C_{\alpha _i\alpha _{0j}}(t_0,t)\alpha _i\alpha _{0j}P(\underset{¯}{\alpha }_0,t_0)P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)𝑑\underset{¯}{\alpha }_0𝑑\underset{¯}{\alpha },$$
(11)
whose evolution equation may be found by taking the time derivative of Eq. (11), substituting Eq. (9) and integrating by parts, obtaining

$$\frac{d}{dt}C_{\alpha _i\alpha _{0j}}(t_0,t)=B_{il}(t_0;t)X_l(t)\alpha _{0j}(t_0),$$
(12)
where we have defined $`X_i(t)\alpha _{0j}(t_0)X_i(\underset{¯}{\alpha },t)\alpha _{0j}P(\underset{¯}{\alpha }_0,t_0)P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)𝑑\underset{¯}{\alpha }_0𝑑\underset{¯}{\alpha }`$. This equation implicitly provides an operational definition of the transport coefficients $`B_{ij}`$, showing aging effects through their dependence on the waiting time. The relation between the dissipation of the energy and the response of the system to an external time-dependent force $`\underset{¯}{ϵ}(t)=\underset{¯}{ϵ}\mathrm{\Theta }(tt_0)`$, with $`\mathrm{\Theta }(x)`$ being the step function, may be determined by means of the evolution equation for the first moment of the probability density, $`\alpha _i(t)=\alpha _iP(\underset{¯}{\alpha }_0,t_0)P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)𝑑\underset{¯}{\alpha }_0𝑑\underset{¯}{\alpha }`$. In the presence of the external force, the system is described by means of the corresponding Fokker-Planck equation

$`{\displaystyle \frac{P(t_0|t)}{t}}={\displaystyle \frac{}{\underset{¯}{\alpha }}}\underset{¯}{\underset{¯}{B}}(t_0;t)\left[\underset{¯}{X}(\underset{¯}{\alpha },t)P(t_0|t)+{\displaystyle \frac{k_BT_B(t)}{m}}A{\displaystyle \frac{P(t_0|t)}{\underset{¯}{\alpha }}}\right]`$
$`{\displaystyle \frac{}{\underset{¯}{\alpha }}}\left[\underset{¯}{ϵ}(t)P(t_0|t)\right],`$ (13)
where the external force has introduced the additional drift $`\underset{¯}{ϵ}(t)`$. According to linear response theory kubo , the response $`R_{\alpha _i\alpha _j}(t_0,t)`$ of the system is then given through the functional derivative

$$R_{\alpha _i\alpha _j}(t_0,t)=\frac{\delta \alpha _i(t)}{\delta ϵ_j(t_0)}.$$
(14)
Notice that, due to the time dependence of the coefficients $`B_{ij}(t_0;t)`$, the response function will contain modifications to the Markovian response.

The response function can be calculated from the equation for the first moment $`\alpha _i(t)`$ which follows from the Fokker-Planck equation (13). We obtain

$$R_{\alpha _i\alpha _{0j}}(t_0,t)=\frac{m}{k_B}\frac{1}{T(t)}\frac{}{t_0}C_{\alpha _i\alpha _{0j}}(t_0,t),$$
(15)
which constitutes the formulation of the fluctuation-dissipation theorem. This expression indicates that the fluctuation-dissipation theorem remains valid when the system is in local quasi-equilibrium; the temperature entering that expression is that of the system at local quasi-equilibrium. If one uses Eqs. (10) and (15), it can also be expressed as

$$R_{\alpha _i\alpha _{0j}}(t_0,t)=\frac{m}{k_B}\frac{A^1}{T_B(t)}\frac{}{t_0}C_{\alpha _i\alpha _{0j}}(t_0,t).$$
(16)
The result (16) can in general be interpreted as a violation of the fluctuation-dissipation theorem kurchan9  and kurchan . When analyzing its equivalent form (15) it expresses validity of the fluctuation-dissipation theorem in terms of the local quasi-equilibrium temperature.

The quantity $`A`$ has been derived in Section II from thermodynamic grounds. An expression in terms of correlation functions can be obtained from the equal-time correlation function $`\mathrm{\Phi }_{\alpha \alpha }(t)\frac{1}{2}\underset{¯}{\alpha }\underset{¯}{\alpha }P(\underset{¯}{\alpha }_0,t_0)P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)𝑑\underset{¯}{\alpha }_0𝑑\underset{¯}{\alpha }`$. After substituting Eq. (9) and integrating by parts, one obtains the evolution equation

$$\frac{d}{dt}\mathrm{\Phi }_{\alpha \alpha }(t)=B(t_0;t)\underset{¯}{X}(t)\underset{¯}{\alpha }(t)+\frac{k_BT(t)}{m}B(t_0;t),$$
(17)
where for simplicity, we have considered $`\underset{¯}{\underset{¯}{B}}(t_0;t)=B(t_0;t)\underset{¯}{\underset{¯}{1}}`$, with $`\underset{¯}{\underset{¯}{1}}`$ the unit tensor. For times $`t\beta ^1(t_0;t)`$, the time derivative of Eq. (17) may be neglected and the system relaxes to the reference state, leading to the expression

$$A\frac{m}{k_BT_B(t)}\underset{¯}{X}(t)\underset{¯}{\alpha }(t)_{qe},$$
(18)
where $`\mathrm{}_{qe}`$ represents the average in respect to the local quasi-equilibrium probability density $`P_{qe}(t)`$. This equation indicates that $`A`$ may in general be expressed by means of the moments of the probability density. According to Eqs. (10) and (18), it follows that the temperature $`T(t)`$ entering the Gibbs equation (2) results from the contribution of the energies related to the corresponding degrees of freedom $`\underset{¯}{\alpha }`$ by

$$T(t)\frac{m}{k_B}\underset{¯}{X}(t)\underset{¯}{\alpha }(t)_{qe}.$$
(19)
An important fact to be noticed from this expression is that the temperature is not robust since it implicitly depends on the number of relevant variables or observables used to describe the system, as has been evidenced in origin . Given that the probability density of the reference state depends on the initial condition, it follows that $`T(t)`$ is also a function of $`t_0`$. In the case when the force is linear, Eq. (19) reduces to $`T(t)\frac{m}{k_B}\alpha ^2(t)_{qe}`$, which coincides with the corresponding result obtained when local thermodynamic equilibrium between system and bath exists. In that case, when $`\underset{¯}{\alpha }=\stackrel{}{u}`$ with $`\stackrel{}{u}`$ the velocity of a particle, Eq. (19) becomes the kinetic definition of the temperature.

## IV Slow dynamics in glassy systems

As a first application of the theory presented, we will analyze the dynamics of systems characterized by two sets of variables: fast $`\beta `$-variables and slow $`\alpha `$-variables. The $`\beta `$-variables relax to a local equilibrium state characterized by the temperature of the bath $`T_B`$, whereas the $`\alpha `$-variables relax to a quasi-equilibrium state characterized by the temperature $`T(t)`$.

The probability density is now a function of the state vector $`\underset{¯}{\mathrm{\Gamma }}=(\underset{¯}{\alpha },\underset{¯}{\beta })`$. Since $`\beta `$-variables equilibrate with the bath, we may assume that the probability density characterizing the reference state can be factorized in the form $`P_{qe}(t)=P_{le}(\underset{¯}{\beta })P_{qe}(\underset{¯}{\alpha },t)`$, with $`P_{le}(\underset{¯}{\beta })`$ and $`P_{qe}(\underset{¯}{\alpha },t)`$ stationary and non-stationary distributions. Consequently, $`s_{qe}(t)=s_{le}^\beta +s_{qe}^\alpha (t)`$, where $`T_B\delta s_{le}^\beta `$ and $`T(t)\delta s_{qe}^\alpha `$ enter Gibbs equations similar to Eq. (2). Consistent with the analysis of Sec. II, the probability densities of the reference states are given by

$$P_{le}(\underset{¯}{\beta })=e^{\frac{m}{k_BT_B}[\mu _{le}^\beta e_{le}(\underset{¯}{\beta })]}\text{and}P_{qe}(\underset{¯}{\alpha },t)=P_ee^{\frac{m}{k_BT(t)}[\mu _{qe}^\alpha e_{qe}(\underset{¯}{\alpha })]}.$$
(20)
By taking the time derivative of Eq. (3) and using Eq. (20) and the corresponding continuity equation, the entropy production can be written as

$$\sigma =P(t_0|t)\underset{¯}{v}(\underset{¯}{\mathrm{\Gamma }},t)\frac{}{\underset{¯}{\mathrm{\Gamma }}}\mathrm{\Theta }(t_0;t)𝑑\underset{¯}{\mathrm{\Gamma }},$$
(21)
where we have defined the nonequilibrium Massieu function

$$\mathrm{\Theta }(t_0;t)\frac{k_B}{m}\mathrm{ln}P(t_0|t)+\frac{e_{le}(\underset{¯}{\beta })}{T_B}+\frac{e_{qe}(\underset{¯}{\alpha })}{T(t)},$$
(22)
which is the appropriate function in the entropy representation callen , N-part .

Since $`\underset{¯}{v}(\underset{¯}{\mathrm{\Gamma }},t)(\underset{¯}{v}_\beta ,\underset{¯}{v}_\alpha )`$ and $`\frac{}{\underset{¯}{\mathrm{\Gamma }}}(\frac{}{\underset{¯}{\beta }},\frac{}{\underset{¯}{\alpha }})`$, from the entropy production (21) we can obtain the expressions for $`\underset{¯}{v}_\beta `$ and $`\underset{¯}{v}_\alpha `$ (see Eq. (7)). In general, the rules of nonequilibrium thermodynamics introduce couplings between the velocities $`\underset{¯}{v}_\beta `$ and $`\underset{¯}{v}_\alpha `$ and the corresponding forces through the mobility $`\underset{¯}{\underset{¯}{B}}`$ which is in general a function of $`t_0`$, $`t`$, $`\underset{¯}{\alpha }`$ and $`\underset{¯}{\beta }`$.

In this case, since no couplings appear due to the separation of time scales, the mobility matrix $`\underset{¯}{\underset{¯}{B}}`$ can be written as the direct sum of two submatrices $`\underset{¯}{\underset{¯}{B}}^\alpha `$ and $`\underset{¯}{\underset{¯}{B}}^\beta `$. For simplicity, we will assume that the $`B_{ij}^\beta `$’s are constant since the $`\beta `$-variables decay to a local equilibrium state degroot . The $`B_{ij}^\alpha `$’s will be functions of two times, $`B_{ij}^\alpha =B_{ij}^\alpha (t_0,t)`$.

Now, by taking into account that $`\beta `$-variables thermalize at $`T_B`$ and $`\alpha `$-variables at $`T(t)`$, tartaglia -theo , after substituting the expressions for $`\underset{¯}{v}_\beta `$ and $`\underset{¯}{v}_\alpha `$ into the continuity equation one obtains the generalized Fokker-Planck equation

$`{\displaystyle \frac{P(t_0|t)}{t}}={\displaystyle \frac{}{\underset{¯}{\alpha }}}\underset{¯}{\underset{¯}{B}}^\alpha (t_0;t)\left[\underset{¯}{X}^\alpha (\underset{¯}{\alpha })P(t_0|t)+{\displaystyle \frac{k_BT(t)}{m}}{\displaystyle \frac{P(t_0|t)}{\underset{¯}{\alpha }}}\right]+`$
$`{\displaystyle \frac{}{\underset{¯}{\beta }}}\underset{¯}{\underset{¯}{B}}^\beta \left[\underset{¯}{X}^\beta (\underset{¯}{\beta })P(t_0|t)+{\displaystyle \frac{k_BT_B}{m}}{\displaystyle \frac{P(t_0|t)}{\underset{¯}{\beta }}}\right],`$ (23)
where, according with Eqs. (8) and (20), $`\underset{¯}{X}^\beta (\underset{¯}{\beta })=\frac{e_{le}(\underset{¯}{\beta })}{\underset{¯}{\beta }}`$ and $`\underset{¯}{X}^\alpha (\underset{¯}{\alpha })=\frac{e_{qe}(\underset{¯}{\alpha })}{\underset{¯}{\alpha }}`$.

One may obtain the Fokker-Planck equation for the slow $`\alpha `$-variables by multiplying Eq. (23) by $`P(\underset{¯}{\mathrm{\Gamma }}_0,t_0)`$ and integrating over $`\underset{¯}{\beta }_0`$ and $`\underset{¯}{\beta }`$ taking into account that the current $`P(t_0|t)\underset{¯}{v}_\beta `$ vanishes at the boundaries. Our result is Eq. (9). In similar form, the evolution equation for the $`\beta `$-variables can be derived by integrating over $`\underset{¯}{\alpha }_0`$ and $`\underset{¯}{\alpha }`$. The resulting Fokker-Planck equation is

$$\frac{P(\underset{¯}{\beta },t)}{t}=\frac{}{\underset{¯}{\beta }}\underset{¯}{\underset{¯}{B}}^\beta \left[\underset{¯}{X}^\beta (\underset{¯}{\beta })P(\underset{¯}{\beta },t)+\frac{k_BT_B}{m}\frac{}{\underset{¯}{\beta }}P(\underset{¯}{\beta },t)\right],$$
(24)
where $`P(\underset{¯}{\beta },t)`$ is the single-time distribution function.

The computation of the moments of the probability density following the procedure indicated in Sec. III, allows one to find the statistical expression of the temperatures

$$T_B\frac{m}{k_B}\underset{¯}{X}^\beta (\underset{¯}{\beta })\underset{¯}{\beta }_{le},$$
(25)
and

$$T(t)\frac{m}{k_B}\underset{¯}{X}^\alpha (\underset{¯}{\alpha },t)\underset{¯}{\alpha }_{qe}.$$
(26)
Eq. (25) implies that the fast degrees of freedom satisfy the equipartition law and its energy is a quadratic function in $`\beta `$’s, whereas (26) is a definition of the temperature of the system in terms of the correlations at the quasi-equilibrium state.

¿From Eq. (23), we may compute the response of the system to an external perturbation and derive the fluctuation-dissipation theorem related to both sets of variables. Following the procedure indicated in Sec. III, the response functions of the fast variables satisfy the relation

$$R_{\beta _i\beta _{0j}}(t_0,t)=\frac{m}{k_B}\frac{1}{T_B}\frac{}{t_0}C_{\beta _i\beta _{0j}}(t_0,t),$$
(27)
where $`R_{\beta _i\beta _{0j}}(t_0,t)=\frac{\delta \beta _i(t)}{\delta ϵ_j(t_0)}`$. The corresponding expression for the response of the $`\alpha `$-variables is given in Eq. (15). From Eq. (27) it follows that the fast degrees of freedom satisfy the fluctuation-dissipation theorem in terms of the temperature of the bath.

## V Brownian motion in a granular gas

As a second application of our formalism, we will discuss in this section the Brownian motion of a grain in a granular gas, which has been previously addressed in dufty1  by means of kinetic theory for systems with inelastic collisions. The phase-space vector is now $`\underset{¯}{\alpha }=(\stackrel{}{u},\stackrel{}{r})`$, where $`\stackrel{}{u}`$ and $`\stackrel{}{r}`$ are the velocity and position of the Brownian grain. The corresponding continuity equation for $`P(t_0|t)`$ is given by

$$\frac{P(t_0|t)}{t}+\stackrel{}{u}P(t_0|t)=\frac{}{\stackrel{}{u}}\left[P(t_0|t)\stackrel{}{v}_\stackrel{}{u}(t_0;t)\right],$$
(28)
We will assume that the bath is in a non-stationary homogeneous cooling state characterized by the kinetic temperature $`T_B(t)`$, dufty1 . The kinetic temperature is given by $`T(t)=\frac{m}{k_B}\frac{u^2}{2}P(\stackrel{}{u}_0,t_0)P(\stackrel{}{u}_0,t_0|\stackrel{}{u},t)𝑑\stackrel{}{u}_0𝑑\stackrel{}{u}`$, and has to be understood as the temperature of a dilute gas of Brownian particles suspended in the granular gas. According to Eq. (4), in this case the local quasi-equilibrium probability density can be taken as

$$P_{qe}(t)=e^{\frac{m}{k_BT(t)}[\mu _{qe}\frac{1}{2}\stackrel{}{u}^2]},$$
(29)
where $`\mu _{qe}`$ normalizes $`P_{qe}(t)`$.

Once the $`P_{qe}(t)`$ has been specified, the Fokker-Planck equation for the Brownian gas may be derived by using Eqs. (3), (28), (29). Following the steps indicated in Sec. II, we finally obtain the generalized Fokker-Planck equation

$$\frac{P(t_0|t)}{t}+\stackrel{}{u}P=\frac{}{\stackrel{}{u}}B(t_0,t)\left[\stackrel{}{u}P(t_0|t)+\frac{k_BT_B(t)}{m}A\frac{P(t_0|t)}{\stackrel{}{u}}\right],$$
(30)
where we have assumed locality in $`\stackrel{}{u}`$-space for which the velocity $`\stackrel{}{v}_\stackrel{}{u}(t_0;t)`$ couples only to the gradient of the chemical potential at the same value of $`\stackrel{}{u}`$ in accordance with Eq. (5). By considering the action of an external force field $`\stackrel{}{ϵ}(t)`$ which can be an impulsive one (tapping force), we may derive by following the steps indicated in Sec. III, the fluctuation-dissipation relation

$$R_{\stackrel{}{u}\stackrel{}{u}_0}(t_0,t)=\frac{m}{k_B}\frac{A^1}{T_B(t)}\frac{}{t_0}C_{\stackrel{}{u}\stackrel{}{u}_0}(t_0,t).$$
(31)
Recent computer simulations of diffusion in dilute granular gases confirm the validity of Eq. (31), loreto . The corresponding integrated version of the generalized fluctuation-dissipation theorem has been also obtained in dense granular systems in Nature-kurchan . Now, by comparing Eq. (30) with the corresponding one found in dufty1  by following kinetic theory methods, we may identify

$$A=\frac{\left(1+\lambda \right)}{2};B(t_0,t)=\frac{1}{2}\left(1+\lambda \right)\gamma _e(t_0,t),$$
(32)
where the restitution coefficient $`\lambda (1)`$ of the Brownian grains has been assumed constant, and $`\gamma _e(t_0,t)`$ is the friction coefficient for the case when elastic collisions take place, dufty1 .

¿From the Fokker-Planck equation (30), it follows that the evolution equation of the second moment of the velocity of the Brownian grains can be expressed in the form

$$\frac{d}{dt}T(t)=2B(t_0;t)\left[T(t)AT_B(t)\right].$$
(33)
which is a particular case of Eq. (17). Eq. (33) establishes that the relation between system and bath temperatures at the local quasi-equilibrium state given by Eq. (10), is valid for times $`tB^1(t_0;t)`$.

The relation between $`T(t)`$ and $`T_B(t)`$ can be generalized in order to describe shorter time scales by incorporating the thermal relaxation of the heat bath through the evolution equation $`\frac{d}{dt}T_B(t)=\zeta (t)T_B(t)`$, with $`\zeta (t)`$ the cooling rate. Assuming that $`0<\zeta (t)<2B(t_0;t)`$, the time evolution of $`\frac{T(t)}{T_B(t)}`$ is then given by dufty1

$$\frac{d}{dt}\frac{T(t)}{T_B(t)}+2\left(1\xi \right)B(t_0;t)\frac{T(t)}{T_B(t)}=2B(t_0;t)A,$$
(34)
where we have defined $`\xi \frac{\zeta (t)}{2B(t_0;t)}`$. From Eq. (34), we then obtain

$$T(t)=\frac{A}{\left(1\xi \right)}T_B(t),$$
(35)
which gives the relation between the system and bath temperatures for times $`t\left[\left(1\xi \right)B(t_0;t)\right]^1`$.

By contracting the description, we enter the diffusion regime described by the density $`\rho (\stackrel{}{r}^{},t_0|\stackrel{}{r},t)=mP(t_0|t)𝑑\stackrel{}{u}`$, whose evolution is governed by a Smoluchowski equation. This equation may be derived by approximating the hierarchy of evolution equations of the first three moments of the probability density $`P(t_0|t)`$, nosotros . One obtains,

$$\frac{}{t}\rho (t_0|t)=D_{eff}(t_0,t)^2\rho (t_0|t),$$
(36)
where $`D_{eff}(t_0,t)`$ is the effective diffusion coefficient given by

$$D_{eff}(t_0,t)\frac{k_BT_B(t)}{mB(t_0,t)}\frac{A}{\left(1\xi \right)^2},$$
(37)
which has been obtained by taking into account Eqs. (34) and (35). Eq. (37) was derived in dufty1  by means of kinetic theory. In particular, when the cooling rate of the bath satisfies $`\zeta (t)2B`$, the expression of the diffusion coefficient (37) reduces to $`D_{eff}(t_0,t)\frac{k_BT_B}{mB(t_0,t)}A`$. In the general case, the dynamics of the system may be described by means of the corresponding Fokker-Planck equation in which the probability density of the reference state has the form $`P_{qe}(t)=e^{\frac{m}{k_BT_B(t)}A^1\left(1\xi \right)[\mu _{qe}\frac{1}{2}\stackrel{}{u}^2]}`$, and the coefficients $`A`$ and $`B(t_0,t)`$ satisfy more general relations than Eq. (32).

The Fokker-Planck and Smoluchowski equations derived previously may also describe diffusion in a colloidal system. In that case, the diffusion coefficient (37) accounts for the violation of the Stokes-Einstein relation reported in dense granular systems Nature-kurchan  and supercooled colloidal fluids coll-glasses .

## VI Discussion

In this paper, we have proposed an extension of Onsager’s theory, originally formulated to describe *nonequilibrium aged* systems, in order to analyze slow relaxation dynamics. Our theory is based on the assumption of the existence of a local quasi-equilibrium state in phase space characterized by a non-stationary probability density. This hypothesis entails the existence of a time-dependent temperature of the system $`T(t)`$ reflecting lack of thermal equilibrium between system and bath, which enables us to formulate a Gibbs equation. The relation between the local quasi-equilibrium and bath temperatures was given in terms of a scale function obtained by taking into account the influence of the bath.

The entropy production derived from the Gibbs entropy postulate leads to the explicit expression of the probability current in phase space and, consequently, to the generalized Fokker-Planck equation (9). Memory effects were incorporated in the description through the time dependence of the coefficients entering into that equation. They are responsible for the dependence of the characteristic relaxation times on the initial (waiting) time. Our formalism could straightforwardly be extended to the case of a description in terms of fields $`\underset{¯}{\alpha }(\stackrel{}{k})`$, hydrofluc , supercooled .

We have shown that the fluctuation-dissipation relation remains valid in terms of the local quasi-equilibrium temperature $`T(t)`$ of the system whose value depends on the time scale considered. The fact that the local quasi-equilibrium temperature depends on time and on the number of variables relevant to characterize the state, clarifies why two systems with slow relaxation dynamics put in thermal contact do not necessarily satisfy a “zero law”.

As a first application, we have analyzed the slow dynamics of glasses. We have shown that the dynamics of the system can be described by a generalized Fokker-Planck equation governing the time evolution of both the fast and slow degrees of freedom. By using this equation, we have verified that the fluctuation-dissipation theorem (27) is satisfied by the fast modes. For the slow modes the theorem (15) is valid in terms of the local quasi-equilibrium temperature.

We have also studied the problem of Brownian diffusion in a dilute granular gas reproducing the Fokker-Planck and Smoluchowski equations derived previously by means of kinetic theory of granular gases. These equations may also be used to describe the dynamics of supercooled colloidal suspensions leading to a generalization of the Stokes-Einstein relation.

The theory presented offers a unified description of systems undergoing slow relaxation and constitutes a generalization of the Onsager theory describing the evolution of fluctuations around nonequilibrium *aging* states.

## VII Acknowledgments

We want to acknowledge Dr. D. Reguera for valious comments. I.S.H. acknowledges UNAM-DGAPA for economic support. This work was partially supported by DGICYT of the Spanish Government under Grant No. PB2002-01267.