text:
# Emergence of chaotic behaviour in linearly stable systems

## 1 Introduction

Unpredictable evolution in dynamical systems is due to the propagation of information. For instance, the sensitivity of trajectories to infinitesimal perturbations implies that any arbitrarily small inaccuracy in the determination of the initial conditions is exponentially amplified in time, with an average rate associated to the positive component of the Lyapunov spectrum. The integral of this component, the so-called Kolmogorov-Sinai entropy , measures the production rate of information that flows from the less significant digits of the dynamical variables to the more significant ones. In particular, the existence of at least one positive Lyapunov exponent is a sufficient condition for identifying chaotic dynamics; conversely, a fully negative spectrum is suggestive of a periodic evolution.

This approach to unpredicatable evolution based on linear stability analysis has been developed in the context of finite (low-dimensional) systems. Its extension to spatially extended dynamical systems is based on the implcit assumption that they can be viewed as a collection of almost independent, finite dimensional subsystems. The existence of a limit Lyapunov spectrum  provides strong support to this hypothesis and typical chaoticity indicators, like entropies and generalized dimensions, can be turned into their corresponding densities . In fact, a primary interest in the study of space-time chaos is the identification of thermodynamic-like properties.

However, Lyapunov instability is not the only source of unpredictability in such systems. Actually, information can also flow through the boundaries and be transmitted in space by nonlinear mechanisms of front propagation. The so-called chaotic rules of Deterministic Cellular Automata (DCA)  are typical examples of unpredictable evolution in the absence of linear instabilities. There, the discreteness of the state variable prevents the very existence of infinitesimale perturbation, while, on the other hand, even isolated “state-flips” may propagate through the lattice with finite velocity, giving rise to an irregular dynamics. In fact, any DCA rule defined over a lattice of $`L`$ cells is bounded to exhibit a periodic behaviour, since the number of possible states is finite ($`b^L`$ if $`b`$ is the number of possible states in each given site). What makes a “chaotic” DCA rule different from an ordered one is the exponential growth of the recurrence time of the typical configurations. Accordingly, the unpredictable behaviour is dynamically persistent only in the infinite-size limit.

A very similar unpredictable behaviour, denoted as stable chaos, has been observed also in coupled map lattice (CML) models  of the type

$$x_i(t+1)=f(\frac{\epsilon }{2}x_{i1}(t)+(1\epsilon )x_i(t)+\frac{\epsilon }{2}x_{i+1}(t)),$$
(1)
where $`x_i(t)`$ is the continuous state variable at time $`t`$ on the site $`i`$ of a 1d lattice; the function $`f`$ is a linearly stable map of the interval  into itself and $`\epsilon `$ is the strength of the diffusive coupling. It is worth pointing out that this spatial coupling cannot produce any linear instability mechanism and the whole spectrum of Lyapunov exponents is found to be negative. Accordingly, the CML dynamics must eventually approach a periodic stable attractor. This notwithstanding, if $`f`$ is equipped with a sufficiently strong nonlinearity (e.g. a discontinuity or a region with rapidly varying slope), one can find a region in parameter space, where the “transient” evolution towards the periodic attractor grows exponentially with the system size $`L`$, analogously to chaotic DCA rules. Despite there is no rigorous proof of this statement, many independent numerical studies confirm such a scenario (see, e.g., ). Striking features of this “transient” regime are its stationarity and apparent ergodicity: for instance, space and time correlation functions decay exponentially like in usual chaotic phases, while ensemble avarages coincide with time averages (provided that a size-independent pre-transient is discarded). Moreover, the maximum Lyapunov exponent is found to approach a stationary negative value, before quite suddenly turning to the value corresponding to the eventual attractor. Because of its exponential growth, the “transient” represent the truly relevant regime in the thermodynamic limit, while the periodic attractor(s) have no practical significance.

A careful inspection indicates that the mechanism of information production in stable chaos is a flow from the outer (left and right) parts of the chain, like in chaotic DCA rules . Accordingly, the unpredictability of stable chaos relies on a genuine nonlinear propagation mechanism,

As we recall in section II, stable chaotic evolution can be detected by measuring damage spreading, i.e. the average velocity of a front propagating into an unperturbed region. This indicator can be viewed as a sort of generalization of the standard Lyapunov exponent and the use of a proper metric, attributing an increasingly smaller weight to the farther sites would make the correspondence more transparent. However, it is worth noting that the front velocity does not allow defining an analogous of a negative Lyapunov exponent, since in the ordered phase a perturbation does not only regress but also decreases everywhere in size leading to the disappearance of the front itself.

Numerical studies of stable chaos have contributed to shed some light on the relationship between this dynamical regime and the appearance of many interesting complex phenomena, such as nonequilibrium phase transitions, spiral chaos, and the propagation of rough interfaces . However, little is rigorously known about the underlying mechanisms. The aim of this paper is precisely to make some progress in this direction by investigating the conditions under which finite perturbations can propagate instead of die out. In analogy to various analytical techniques introduced to estimate the maximum Lyapunov exponent in standard chaotic CMLs , here we assume that the evolution in phase-space generates a truly random pattern, characterized by short range correlations and thereby introduce a suitable stochastic model describing the evolution in the difference space. Besides diffusion, the model dynamics allows for a random alternancy of a contraction and an expansion process the probability of which depends on the perturbation size. First of all, in Sec, II we verify that, with an appropriate choice of the parameter values, the model is able to reproduce not only qualitatively but also quantitatively the main features of stable chaos. Afterwards, we further simplify the rule determining the expansion process to reduce the subsequent technicalities, while keeping untouched the key ingredients of the model.

In order to understand how a finite perturbation can be sustained even in the presence of an average contraction rate, in Sec. III we first discuss the uncoupled, i.e. 0-dimensional, case, where, by definition, propagation of perturbations is absent. The negative Lyapunov exponent obviously implies that perturbations are eventually absorbed, so that stable chaos cannot exist in this framework. Nonetheless, the presence of a non-uniform contraction process yields nontrivial properties of the dynamics. They are exemplified by the difference existing between the standard multifractal spectrum (associated to infinitesimal perturbations) and the spectrum defined in this paper to describe the evolution of finite perturbations. In fact, a noteworthy result of Sec. III is that it is possible to define a finite-size multifractal spectrum independently of the initial amplitude of the perturbation. Finally, in Section III we comment about the connection with the finite-size Lyapunov exponents recently proposed by some authors as a tool to characterize dynamical unpredictability, beyond standard linear stability analysis .

The methods and concepts introduced in Section III are applied in Section IV to the study of the spatially extended version of the stochastic model. For the sake of space, we limit our discussion to the case of “democratic” coupling, but it is clear that stable chaos arises in a broad region of strong coupling. It is precisely the diffusive coupling to be responsible for the sustainment of stable chaos. Over small scales, contraction is more effective than the sporadic amplification: in such circumstances, diffusion just levels the damping process. Conversely, at larger scales, diffusion proves to be an effective mechamism to propagate locally generated amplifications. The net effect is that, in suitable parameter regions, perturbations self-sustain.

By introducing a factorization hypothesis of the spatial degrees of freedom, we obtain a good estimate of the probability distribution associated with the spatially extended dynamics. In particular, we find a critical value of the contraction rate separating a dynamical regime, where any perturbation eventually vanishes, from a truly chaotic phase, where the average value of the perturbation remains finite in the infinite-time limit. A mean field argument provides a suggestive description of this scenario: the diffusive coupling induces a shift of the finite-size Lyapunov exponent that grows from the negative quenched average (corresponding to the uncoupled limit) towards the annealed average of the expansion rate. Depending on the parameter values, this latter average may be strictly larger than 0, thus implying that the zero-amplitude regime is unstable.

The transition exibited by the stochastic model is reminiscent of the “fuzzy” transition region found in a deterministic CML . However, we cannot push the analogy to a quantitative level, since the stochastic CML dynamics is self-generated and thus it is increasingly regular while approaching the transition region (this is the main reason for the fuzziness observed in Ref. ), while here the stochastic properties of the contraction/expansion process is fixed a priori. A tight analogy, instead, exists with the stochastic synchronization induced by an additive noise . Indeed, in this context, the external noise does not change its robust stochastic features when passing from the synchronized to non-synchronized regime. It is precisely this analogy which suggests that the transition described in this paper should belong to the universality class of directed percolation . The numerical simulations described in Sec. V do confirm such an expectation, but subtle problems still prevent us from still drawing definite conclusions. Such open problems and the possible future perspectives are summarized in Sec. V.

## 2 Generalities

The typical functions used to investigate stable chaos in CML are piecewise linear maps of the unit interval of the form (see also Fig. (1) for a pictorial representation)

$`f(x)=\{\begin{array}{cc}c_1x\hfill & \hfill 0x<1/c_1\\ 1(x1/c_1)(1c_2)/\mathrm{\Delta }\hfill & \hfill 1/c_1x<1/c_1+\mathrm{\Delta }\\ c_2+c_3[x(1/c_1+\mathrm{\Delta })]\hfill & \hfill 1/c_1+\mathrm{\Delta }x1\end{array}.`$ (5)
The most studied case in the literature  corresponds to the limit $`\mathrm{\Delta }0`$, where $`f(x)`$ reduces to a discontinuous function. As we shall comment along this paper, a strong nonlinear component of map $`f(x)`$, rather than a true discontinuity, is sufficient to yield stable chaos. This is why we prefer to consider here the more general case (5). The parameter range of interest for the present study is when all initial conditions (except for a set of zero Lebesgue measure) converge to the same periodic orbit. For instance, for $`c_1=2.7`$, $`c_2=0.07`$, $`c_3=0.1`$ and $`\mathrm{\Delta }=0.01`$, a stable period-3 orbit exists with Lyapunov exponent $`\lambda =0.316..`$. As already mentioned in the Introduction, there is a range of values of the diffusive coupling $`\epsilon `$ defined in (1), where the CML dynamics exhibits a “chaotic” evolution, despite the largest Lyapunov exponent is negative. Damage spreading analysis provides a first hint about this mechanism, that is responsible for the sustainment of irregular behaviour. More precisely, while standard chaos amounts to a flow of information from the less to the more significant digits, stable chaos is generated by a flow of information from the outer (left and right) to the inner parts of an infinite chain. Unfortunately, while the former flow can be “easily” studied thanks to the linearity of the process (in fact, over sufficiently small scales, any smooth function can be linearized), the same argument does not apply to the dynamics on the left and right edges, that is equally nonlinear at any spatial position.

This is the main reason for the difficulty in deriving the necessary and sufficient conditions for the propagation of perturbations. Moreover, while the dynamics of an infinitesimal perturbation can be studied by neglecting propagation phenomena, the opposite is not possible. One cannot study propagation without properly accounting for the local contraction/expansion mechanisms.

The damage spreading analysis is performed by studying the dynamical variable $`u_i(t)=|x_i(t)y_i(t)|`$, i.e. the absolute value of the difference between two test trajectories, $`x_i(t)`$ and $`y_i(t)`$, that are initially set equal to one another on the right of some lattice site, say $`i=0`$, while they are assumed to be totally independent on the left. With such a setting, the damage spreading analysis amounts to studying the propagation of a front separating the region in space where $`u_i(t)𝒪(1)`$ (the tail) from the region where the two trajectories converge to each other, $`u_i(t)0`$ (the forefront).

The first conceptual problem that we have to face is not just the propagation of the front, but its self-sustainment in spite of the local average contraction rate. In order to shed some light on this crucial point, we have simplified the model by assuming that the dynamics is indeed irregular, thereby determining whether this assumption is consistent with the sustainment of an $`𝒪(1)`$ perturbation. This is analogous to the consistency approach developed for the description of standard chaos, where complete randomness of the multipliers of the evolution operator in the tangent space is assumed in order to estimate the maximum Lyapunov exponent .

Accordingly, we introduce a suitable stochastic model to describe the evolution in the “difference” space spanned by $`u_i(t)`$. Let us start from the simple case $`\mathrm{\Delta }=0`$. With reference to the CML dynamics, the dynamical rule is composed of two steps. The first one corresponds to the application of the standard discrete diffusive operator

$$\stackrel{~}{u}_i(t)=(1\epsilon )u_i(t)+\frac{\epsilon }{2}u_{i+1}(t)+\frac{\epsilon }{2}u_{i1}(t),$$
(6)
where $`0\epsilon 1`$ is the coupling parameter. The second step contains the stochastic component of the evolution rule,

$`\begin{array}{c}u_i(t+1)=\{\begin{array}{cc}r_i(\stackrel{~}{u}_i,t),\hfill & \text{w.p.}p=\mathrm{min}[1,b\stackrel{~}{u}_i(t)]\hfill \\ a\stackrel{~}{u}_i(t),\hfill & \text{w.p.}1p\hfill \end{array},\end{array}`$ (10)
where, “w.p.” is the shorthand notation for “with probability”, while $`r_i(\stackrel{~}{u}_i,t)`$ is a random number distributed in the unit interval according to some probability distribution that depends on $`\stackrel{~}{u}_i`$.

This is the non-trivial part of the stochastic model, defined to mimick the evolution of perturbation in the CML dynamics. The first line describes the instability mechanism associated with the discontinuity while the second line describes the contraction of $`u_i(t)`$ by a constant factor 111For the sake of simplicity we assume that the contraction rate is a constant, as in . $`0<a<1`$ . In the CML model, the instability mechanism arises whenever the test trajectories, $`x_i(t)`$ and $`y_i(t)`$, lie on different sides of the discontinuity of the map. In this case, the value taken by the difference variable $`u_i(t+1)`$ is not uniquely determined by the value of $`u_i(t)`$, since it depends also on $`x_i(t)`$. In particular, for small values of $`u_i(t)`$, the instability mechanism occurs quite rarely (it is unlike that the discontinuity is placed across two nearby trajectories) and $`u_i(t+1)`$ is amplified by a big factor; conversely, for large values of $`u_i(t)`$, the discontinuity plays a role much more frequently, but it is less effective. The numerical analysis of the CML model shows that the probability $`p`$ of the instability mechanism grows linearly for small values of the difference variable $`u_i(t)`$ and approaches $`1`$ when $`u_i(t)`$ does the same. Accordingly, in the stochastic model (10), we have decided to schematize this dependence with the simple law $`p=\mathrm{min}[1,b\stackrel{~}{u}_i(t)]`$.

Such a stochastic model can be straightforwardly generalized to cover the case corresponding to a CML dynamics with a nonzero $`\mathrm{\Delta }`$. The main difference is that, whenever $`\stackrel{~}{u}`$ is smaller than $`\mathrm{\Delta }`$, $`u`$ is expanded by a fixed factor $`\mathrm{\Delta }^1`$ with a constant probability $`b\mathrm{\Delta }`$, while the previous stochastic rule still applies for $`\stackrel{~}{u}>\mathrm{\Delta }`$. This amounts to assuming that only perturbations larger than $`\mathrm{\Delta }`$ perceive the steep branch of the map as an effective discontinuity. Accordingly, Eq. (10) is replaced by

$`\begin{array}{cc}u_i(t+1)=\{\begin{array}{cc}r_i(\stackrel{~}{u}_i,t),\hfill & \text{w.p.}p=\mathrm{min}[1,b\stackrel{~}{u}_i(t)]\hfill \\ a\stackrel{~}{u}_i(t),\hfill & \text{w.p.}1p\hfill \end{array},\hfill & \text{if }\stackrel{~}{u}_i(t)>\mathrm{\Delta }\text{,}\\ u_i(t+1)=\{\begin{array}{cc}\stackrel{~}{u}_i(t)/\mathrm{\Delta },\hfill & \text{w.p.}p=b\mathrm{\Delta }\hfill \\ a\stackrel{~}{u}_i(t),\hfill & \text{w.p.}1p\hfill \end{array}\hfill & \text{if }\stackrel{~}{u}_i(t)\mathrm{\Delta }\text{.}\end{array}`$ (17)
One can easily check that these formulae reduce to Eq. (10) in the limit $`\mathrm{\Delta }0`$.

The maximum Lyapunov exponent of the CML dynamics corresponds to the time average of the expansion rates of infinitesimal perturbations. In the context of the stochastic model, it naturally corresponds to the quantity

$$\lambda _0=\underset{t\mathrm{}}{lim}\underset{u(0)0}{lim}\frac{1}{t}\underset{\tau =0}{\overset{t1}{}}\mathrm{ln}\frac{u_i(\tau +1)}{u_i(\tau )}=\underset{t\mathrm{}}{lim}\underset{u(0)0}{lim}\frac{1}{t}\mathrm{ln}\frac{u_i(t)}{u_i(0)},$$
(18)
that we still denote as the “maximum Lyapunov exponent”.

An analytical estimate of $`\lambda _0`$ can be obtained by a mean-field argument, according to which the probability of applying the expansion factor $`1/\mathrm{\Delta }`$ is $`b\mathrm{\Delta }`$, while the probability of appling the contraction factor $`a`$ is $`1b\mathrm{\Delta }`$. One easily obtains

$$\lambda _0b\mathrm{\Delta }\mathrm{ln}\left(\frac{1}{\mathrm{\Delta }}\right)+(1b\mathrm{\Delta })\mathrm{ln}a=\mathrm{ln}ab\mathrm{\Delta }\mathrm{ln}(a\mathrm{\Delta }).$$
(19)
Numerical simulations indicate that the true value of $`\lambda _0`$ is generally slightly larger than this mean-field estimate.

Moreover, for $`\mathrm{\Delta }=0`$, $`\lambda _0\mathrm{ln}a<0`$, while for increasing values of $`\mathrm{\Delta }`$ it may become positive, indicating that a standard choatic regime is attained222This is not entirely correct, since, as discussed in Ref. , the standard chaotic regime occurs when a stronger constraint is met: the linear propagation velocity coincides with the nonlinear one.. Here, we are interested in studying only the parameter region where $`\lambda _0`$ remains negative.

### 2.1 Comparison between stochastic and deterministic models

The reliability of the stochastic models (10) and (17) has been tested by numerical simulations which show that both rules (10) and (17) together with the diffusive coupling (6), exhibit the same qualitative features of the CML dynamics (1), (5). In particular, when the contraction is relatively strong, any initial perturbation $`u_i(0)`$ is quickly absorbed to the fixed point $`u_i=0i`$: this regime corresponds to the “ordered” or, equivalently, to the “synchronized” phase. On the other hand, for weaker contractions, almost any initial condition evolves towards an irregular spatial structure like the one shown in Fig. 2. In this case, the space average

$$\overline{u}(t)=\frac{1}{L}\underset{i=1}{\overset{L}{}}u_i(t)$$
(20)
remains finite and $`\overline{u}(t)=u^{}`$, independently of $`t`$ and $`L`$ for sufficiently large sizes. To be more precise, perturbations of the order $`𝒪(1)`$ survive only for a finite time $`\tau `$ also in the chaotic phase of the stochastic model. Nevertheless, in analogy to the CML model, $`\tau `$ grows exponentially with the lattice size $`L`$. Furthermore, ensemble averages indicate that this regime is asymptotically stationary in time and the comparison with time averages (obviously performed over times much shorter than $`\tau `$) indicate that ergodicity holds as well. It is therefore meaningful to define the single-site probability distribution $`Q(t,v)`$ of finding a perturbation $`u_i`$ in between $`v`$ and $`v+dv`$ at time $`t`$, and its stationary limit $`Q(v)`$, attained for large enough times. It is also worth introducing the space and ensemble averages $`m(t)=\overline{u}(t)`$ (from here on, unless otherwise stated, $``$ denotes an ensemble average) and the stationary limit, $`m(t)\stackrel{~}{m}`$.

In the numerical investigations we have assumed no-flux boundary conditions as they preserve the $`u=0`$ fixed point and are reasonably harmless in the chaotic phase. The $`\stackrel{~}{u}`$-dependent probability distribution of the random variable $`r_i(\stackrel{~}{u}_i,t)`$ has been reconstructed from the CML dynamics. Moreover, in order to avoid boundary effects in the damage spreading analysis, the lattice size has always been chosen in such a way that the perturbation front never reaches the lattice edges during each simulation.

In order to investigate damage spreading phenomena, the initial conditions have been fixed by imposing $`u_i(0)=0`$ for $`i>0`$ and randomly choosing $`u_i(0)`$ with a uniform probability distribution in the interval $`[0,1]`$ for $`i0`$. The statistical fluctuations showed by $`u_i(t)`$ have been smoothed out by performing ensemble averages (over different realizations of the stochastic process) of the spatial configurations at equal times. In the chaotic phase, the initial “kink”-like structure persists: the front connecting the perturbed with the synchronized region moves with a fluctuating velocity. The ensemble averages corresponding to three different times are reported in Fig. 3b, where they are compared with the results directly obtained from the CML model (for the sake of space we limit ourselves to consider the case $`\mathrm{\Delta }=0`$a), where ensemble averages are performed over different initial conditions. Several observations are in order. First, we notice that, in spite of the simplifications introduced in the stochastic model (besides the lack of space-time correlations, we have indeed assumed a constant contraction rate as if the two branches of the local map (5) had the same slope), there is a reasonable agreement with the CML data. In particular, it can be seen that the average height $`u^{}`$ is definitely smaller than 1 in both cases: the reason can be traced back to the combined effect of the contraction mechanism with the diffusive process. Somehow larger differences can be observed in the behaviour of the variance $`\sigma _u^2=(u_i(t))^2u_i(t)^2`$ but they can be attributed more to the approximation in the description of the perturbed region, rather than to peculiarities of the propagation.

In fact, if we look at the width $`w(t)`$ of the average front, defined as the distance between the rightmost sites where $`u_i(t)`$ is larger than $`0.8u^{}`$ and, respectively, $`0.2u^{}`$, we observe a nearly square-root growth in both the CML and the stochastic model (see Fig. 4, where a nearly quantitative agreement is also observed). Moreover, we have verified that, keeping the contraction parameter $`a`$ fixed, and increasing $`\mathrm{\Delta }`$ up to small values (e.g. $`\mathrm{\Delta }0.01`$), $`v_F`$ slightly increases in both models.

On the basis of the results reported in Fig. 4, one can summarize our observations of the front dynamics by effectively assuming a stepwise shape for the profile at any time and approximating its motion with a diffusive process with drift. More precisely, we can write

$`u_i(t)`$ $``$ $`{\displaystyle \frac{u^{}}{\sqrt{2\pi Dt}}}{\displaystyle _0^+\mathrm{}}\mathrm{\Theta }[x(iv_Ft)]\mathrm{exp}\left({\displaystyle \frac{x^2}{2Dt}}\right)𝑑x=`$ (21)
$`=`$ $`{\displaystyle \frac{u^{}}{2}}\left[1\text{erf}\left({\displaystyle \frac{iv_Ft}{\sqrt{2Dt}}}\right)\right],`$
where $`v_F`$ is the average front velocity, while $`D`$ is the “diffusion” coefficient accounting for the square-root growth in time for the standard deviation of the front position $`x`$ (assumed to be continuous, for the sake of simplicity).

### 2.2 Further simplifications of the stochastic models

As we have verified that the qualitative features of the front dynamics do not depend on the shape of the probability distribution of $`r_i(\stackrel{~}{u}_i,t)`$, we have decided to simplify the stochastic models by assuming a $`\delta `$-like distribution, i.e. that $`r_i(\stackrel{~}{u}_i,t)=1`$. Moreover, in order to get rid of unnecessary technical complications we shall assume that $`b=a`$. Notice also that, having chosen $`r_i(\stackrel{~}{u}_i,t)=1`$, we are obliged to assume $`b<1`$, in order to avoid the appearence of a fictitious fixed point $`u^{}=1`$. From now on we study the following model

$`\stackrel{~}{u}_i(t)=(1\epsilon )u_i(t)+{\displaystyle \frac{\epsilon }{2}}u_{i+1}(t)+{\displaystyle \frac{\epsilon }{2}}u_{i1}(t),`$
$`\begin{array}{cc}u_i(t+1)=\{\begin{array}{cc}1,\hfill & \text{w.p.}p=a\stackrel{~}{u}_i(t)\hfill \\ a\stackrel{~}{u}_i(t),\hfill & \text{w.p.}1p\hfill \end{array},\hfill & \text{if }\stackrel{~}{u}_i(t)>\mathrm{\Delta }\text{,}\\ u_i(t+1)=\{\begin{array}{cc}\stackrel{~}{u}_i(t)/\mathrm{\Delta },\hfill & \text{w.p.}p=a\mathrm{\Delta }\hfill \\ a\stackrel{~}{u}_i(t),\hfill & \text{w.p.}1p\hfill \end{array}\hfill & \text{if }\stackrel{~}{u}_i(t)\mathrm{\Delta }\text{,}\end{array}`$ (28)
that we call Continuous Stochastic Model (CSM). In the limit $`\mathrm{\Delta }0`$ it will be named Discontinuous Stochastic Model (DSM).

The only parameter that we are going to consider in the following sections is the contraction parameter $`a`$: this is sufficient to identify and characterize the relevant transition from the chaotic/unsynchronized phase to the ordered/synchronized one.

## 3 The zero dimensional case

The most important problem of linearly stable chaos is to understand how finite perturbations can propagate in spite of the average local contraction. To clarify this point, in this section we consider the zero-coupling limit, i.e. the 0-dimensional case. We shall see that even if “chaotic” motion cannot be sustained whenever the Lyapunov exponent is negative, the lack of a uniform contraction induces anyhow non-trivial properties.

Our first observation concerns the probability $`P(t,u(0))`$ for a finite-amplitude perturbation $`u(0)`$ to be never amplified by the instability mechanism over a time $`t`$. In the DSM, such a probability can be easily factorized as

$`P(t,u(0))={\displaystyle \underset{n=1}{\overset{t}{}}}(1u(0)a^n)=\mathrm{exp}\left[{\displaystyle \underset{n=1}{\overset{t}{}}}\mathrm{ln}(1u(0)a^n)\right]=`$
$`=\mathrm{exp}\left[{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{u(0)^k}{k}}{\displaystyle \underset{n=1}{\overset{t}{}}}a^n\right]=\mathrm{exp}\left[{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{u(0)^k}{k}}{\displaystyle \frac{a^k(1a^{kt})}{1a^k}}\right]`$
$`\stackrel{t\mathrm{}}{}\mathrm{exp}\left[{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(u(0)a)^k}{k(1a^k)}}\right]:=\stackrel{~}{P}(u(0)).`$ (29)
One can see that, for $`t\mathrm{}`$, $`P(t,u(0))`$ approaches a finite value $`\stackrel{~}{P}(u(0))`$ that is both strictly larger than 0 and smaller than 1 for any value of $`u(0)`$ (if $`0<a<1`$). The same conclusion can be drawn also for the CSM, although the algebra is more complicate in that case. The inequality $`\stackrel{~}{P}(u(0))>0`$ indicates that the occasional amplifications are not so strong as to prevent the eventual absorption of the perturbation (this is consistent with our goal to deal with linearly stable processes). On the other hand, the inequality $`\stackrel{~}{P}(u(0))<1`$ indicates that the amplification process cannot be neglected. Notice that, since this holds true independently of $`u(0)`$ (although $`\stackrel{~}{P}(u(0))1`$ for $`u(0)0`$), there is a difference with the convergence to a stable fixed point in a topologically chaotic map (think, e.g., of the logistic map in one of the stability windows that follow the first period doubling) since, in that case, there would be a threshold (corresponding to the border of the basin of attraction), below which a monotonous contraction would start. A closer similarity exists with the so-called strange nonchaotic attractors, as they are characterized by a nonmonotonous contraction even arbitrarily close to the attractor .

The multifractal theory  provides the most appropriate framework to characterize this system. By following this approach, devised with reference to infinitesimal perturbations, we introduce the exponential growth rate of a finite initial perturbation $`u(0)`$,

$$\mathrm{\Lambda }(t,u(0))=\frac{1}{t}\mathrm{ln}\left(\frac{u(t)}{u(0)}\right).$$
(30)
In the limit $`u(0)0`$, $`\mathrm{\Lambda }(t,u(0))\lambda (t)`$, the equivalent in this context of the usual finite-time Lyapunov exponent . The proper indicator to look at is the probability distribution $`𝒫(\mathrm{\Lambda },t,u(0))`$ to find a growth rate between $`\lambda `$ and $`\mathrm{\Lambda }+d\mathrm{\Lambda }`$ at time $`t`$ starting from $`u(0)`$. More precisely, we introduce the finite-size multifractal spectrum

$$H(\mathrm{\Lambda })=\underset{t\mathrm{}}{lim}\left[\frac{1}{t}\mathrm{ln}𝒫(\mathrm{\Lambda },t,u(0))\right].$$
(31)
As it will become clear later, $`H(\mathrm{\Lambda })`$ is independent of $`u(0)`$ (as long as $`u(0)>0`$) . Nevertheless, we shall show that it differs from the standard multifractal spectrum $`h(\lambda )`$, obtained by taking the limit $`u(0)0`$ before the infinite-time limit

$$h(\lambda )=\underset{t\mathrm{}}{lim}\underset{u(0)0}{lim}\left[\frac{1}{t}\mathrm{ln}𝒫(\mathrm{\Lambda },t,u(0))\right].$$
(32)
In other words, the order of the two limits is crucial for understanding the difference between the behaviour of infinitesimal and finite perturbations. In the following, we shall compare the multifractal distribution of $`\lambda (t)`$ (linear analysis) with that of $`\lambda (t,u(0))`$ (nonlinear analysis) in both stochastic models.

### 3.1 Linear analysis

We start from the CSM, as the DSM is nothing but a limit case of the former one. The linear approximation amounts to assuming $`u_i(t)<\mathrm{\Delta }i,t`$. In this case, standard combinatorial analysis implies

$`h_\mathrm{\Delta }(\lambda )`$ $`=`$ $`{\displaystyle \frac{\mathrm{ln}(\mathrm{\Delta })+\lambda }{\mathrm{ln}(a\mathrm{\Delta })}}\left[\mathrm{ln}(1a\mathrm{\Delta })\mathrm{ln}\left({\displaystyle \frac{\mathrm{ln}(\mathrm{\Delta })+\lambda }{\mathrm{ln}(a\mathrm{\Delta })}}\right)\right]`$ (33)
$`+{\displaystyle \frac{\mathrm{ln}(a)\lambda }{\mathrm{ln}(a\mathrm{\Delta })}}\mathrm{ln}\left[{\displaystyle \frac{a\mathrm{\Delta }\mathrm{ln}(a\mathrm{\Delta })}{\mathrm{ln}(a)\lambda }}\right],`$
where we have made explicit the dependence on the parameter $`\mathrm{\Delta }`$. In the limit $`\mathrm{\Delta }0`$, the above expression reduces to

$$h_0(\lambda )=\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}𝒫[\lambda (t)]\stackrel{\mathrm{\Delta }0}{}\mathrm{ln}(a)\lambda ,\lambda >\mathrm{ln}a.$$
(34)
On the other hand, performing the limit $`\mathrm{\Delta }0`$ before the $`t\mathrm{}`$ limit would yield the trivial result $`h_0(\lambda )=0`$ with the support of $`h_0`$ restricted to the point $`\lambda =\mathrm{ln}a`$. Thus, the non-commutativity of the limits $`\mathrm{\Delta }0`$ and $`t\mathrm{}`$ reveals that the discontinuous case is a singular limit of the CSM. In other words, the multifractal spectrum of the discontinuous model depends on the way it is defined. We prefer to adopt a “physical” point of view, i.e. to consider the discontinuity as the limit of a negligible $`\mathrm{\Delta }`$, which corresponds to taking first the $`t\mathrm{}`$ limit.

### 3.2 Nonlinear analysis

We now consider a finite perturbation $`u(0)`$ in the simple context $`\mathrm{\Delta }=0`$. From Eq. (30) we see that $`\mathrm{\Lambda }`$ is uniquely determined from the knowledge of $`u(t)`$ and $`u(0)`$ (and, obviously, the time $`t`$). Accordingly, the knowledge of $`𝒫(\mathrm{\Lambda },t,u(0))`$ is fully equivalent to that of the single-site probability distribution $`Q(t,u)`$, together with the initial condition

$$Q(0,u)=\delta (uu(0)).$$
(35)
It proves useful to introduce the notation $`u=a^n`$ and $`u(0)=a^{n_0}`$, where $`n`$ and $`n_0`$ are real variables $`0`$. Eq. (30) can be rewritten as

$$\lambda =\lambda (t,n,n_0)=\frac{1}{t}\mathrm{ln}\left(\frac{a^n}{a^{n_0}}\right)=\frac{nn_0}{t}\mathrm{ln}a.$$
(36)
Once $`Q(t,n)`$ is known, Eq. (36) allows one reconstructing the corresponding probability distribution $`𝒫(\mathrm{\Lambda },t,u(0))`$. Two possibilities are in order, either the system has never been “kicked”, i.e. reset to 1 by the instability mechanism, in which case the initial value has been contracted $`t`$-times by a factor $`a`$, or it has received at least one kick, loosing memory of the initial condition. In the former case, occurring with probability $`P(t,u(0))`$, $`\mathrm{\Lambda }(t,u(0))=\mathrm{ln}a`$ or, equivalently, $`n=n_0+t`$. In the latter case, $`\mathrm{\Lambda }(t,u(0))>\mathrm{ln}a`$ and the accessible values of $`n`$ are restricted to positive integer numbers strictly smaller than $`n_0+t`$ (since that the maximum possible contraction factor in $`t`$ time steps is $`a^t`$). Accordingly, $`n`$ can be interpreted as the elapsed time since the last kick, and the probability distribution $`Q(t,u)`$ can be factorized as the product of the probability $`G(tn,n_0)`$ of receiving a kick at time $`tn`$ for an initial perturbation $`u(0)=a^{n_0}`$ by the probability $`P(n,1)`$ of not being kicked anymore for the remaining $`n`$ time steps,

$$Q(t,n)=G(tn,n_0)P(n,1)$$
(37)
(in the remaining part of this section and in App. A, with no ambiguity arising, we denote with $`Q`$ also the probability density of the logarithmic variable $`n`$). Moreover, we impose the condition $`P(0,u(0)):=1`$ to extend the validity of Eq. (37) to the case $`n=0`$.

The probability $`G(t,n_0)`$ can be recursively expressed as the probability of receiving the very first kick at time $`t`$ plus the probability of receiving the second last kick at any previous time, i.e.

$$G(t,n_0)=a^{t+n_0}\underset{k=1}{\overset{t1}{}}(1a^{k+n_0})+\underset{k=1}{\overset{t1}{}}G(k,n_0)a^{tk}P(t1k,1),$$
(38)
with $`G(1,n_0):=a^{n_0+1}`$. Eq. (38) can be numerically iterated in the large time limit to obtain the expression of the multifractal distribution (31). The analysis performed in App. A shows that $`H(\mathrm{\Lambda })`$, defined as in (31), is a segment of straight line restricted to the open interval $`(\mathrm{ln}a,0)`$ of negative values. The slope of this straight line depends on the contraction parameter $`a`$ but is independent of $`u(0)`$. In fact, in the appendix we show that the slope can be obtained by solving an eigenvalue problem, where $`u(0)`$ enters to specify the initial condition but not the operator itself. An approximate analytic solution is also determined, which confirms the numerical observation that the slope increases monotonously from -1 (for $`a0`$) to 0 (for $`a1`$). Curves 2 and 4 in Fig. 5 correspond to the nonlinear and linear analysis of the discontinuous case for $`a=0.7`$, respectively. The latter curve lies well below the former one, indicating that the linear analysis leads to an underestimation of the fluctuations. This is a general fact holding for all values of $`a`$ in the meaningful range $`[0,1]`$. The difference must be attributed to the sporadic amplifications due to the discontinuity: it is remarkable that the finite-size spectrum is independent of the initial condition.

It is now important to test whether the difference between linear and nonlinear curves persists also when the discontinuity is removed, i.e. in the CSM. For the sake of simplicity, we fix $`\mathrm{\Delta }`$ equal to $`a^{\overline{n}}`$ and suppose that both $`\overline{n}`$ and $`n_0`$ are non-negative integers, so that $`u(t)`$ is defined on a discrete subset of the unit interval, $`\{a^n\}_{n=0,1,\mathrm{}}`$. This assumption does not affect the main conclusions while it allows one writing a simple recursive equation for the probability distribution $`Q(t,n)`$:

$`Q(t+1,0)`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{\overline{n}}{}}}Q(t,k)a^{k+1},`$
$`Q(t+1,n)`$ $`=`$ $`\{\begin{array}{cc}(1a^n)Q(t,n1)+a^{\overline{n}+1}Q(t,n+\overline{n})\hfill & 0<n\overline{n}\\ (1a^{\overline{n}+1})Q(t,n1)+a^{\overline{n}+1}Q(t,n+\overline{n})\hfill & n>\overline{n}\end{array}`$ (41)
Equations (41), with the initial conditions $`Q(0,n)=\delta _{n,n_0}`$ and boundary conditions $`Q(t,n)=0`$ $`t>n+n_0`$, can be numerically iterated to obtain the finite-size multifractal distribution in the continuous case.

Since the development of analytical techniques to determine $`H_\mathrm{\Delta }`$ is by far more complex than in the previous case, we have limited ourselves to determine the multifractal spectrum numerically. The linear and nonlinear spectra are reported in Fig. 5 for $`\mathrm{\Delta }=0.0097`$ (curves 3 and 1, respectively). Let us first notice that the finite-size spectrum is again independent of the initial condition $`u(0)`$ and this makes $`H_\mathrm{\Delta }`$ a well defined quantity: the effect of $`u(0)`$ is just to affect the convergence to the asymptotic spectrum. Furthermore, we can see that the nonlinear curve lies well above the standard multifractal spectrum, indicating that it is not the discontinuity to be responsible for the difference between the behaviour of infinitesimal and finite-size perturbations. Moreover, the overall closeness of curves 1 and 2 reveals that the removal of the discontinuity does not introduce significant differences in the finite-size spectrum.

Finally, let us closely compare $`H_\mathrm{\Delta }(\mathrm{\Lambda })`$ with $`h_\mathrm{\Delta }(\lambda )`$: the inset in Fig. 5 reveals that the two coincide for $`\mathrm{\Lambda },\lambda <\lambda _c`$ (see the dashed line), while above some critical value $`\lambda _c`$ the finite-size spectrum continues as a straight line, while $`h_\mathrm{\Delta }`$ decreases faster. In the language of thermodynamic formalism, the linear behaviour of $`H_\mathrm{\Delta }`$ is suggestive of a phase transition, from small values of $`\mathrm{\Lambda }`$, that are correctly described by the linear analysis, to large values of $`\mathrm{\Lambda }`$, where the finite character of the perturbation cannot be neglected. This point would certainly require a more detailed analysis to provide a more solid background to the above arguments, but we avoid this as it would drive us too far from what is the main goal of the present paper.

The analysis performed in this section has allowed us introducing well defined observables to deal with finite-size perturbations. We cannot, however, avoid commenting on an alternative class of tools that have been devised to deal with this problem, although mainly in the context of truly chaotic systems. Since in realistic physical conditions perturbations are always finite, it is very tempting to introduce a growth rate to characterize perturbations of different sizes. However, serious conceptual problems are immediately encountered if one tries to define truly finite-size Lyapunov exponents. On the one hand, the finiteness of the size induces a dependence on the norm utilized, on the other hand it requires that the Lyapunov exponent has to be defined for a finite-time resolution, since perturbations change size over time. This latter implication is rather crucial in that so defined Lyapunov exponents are not self-averaging quantities (see Ref.  for a detailed discussion of the problem).

In spite of such limitations, finite-size Lyapunov exponents may carry useful information, although one has to be careful in interpreting them. Let us, for instance, look at

$$\mathrm{\Lambda }_0(u)=\mathrm{ln}\frac{u^{}}{u}$$
(42)
where $`u^{}`$ is the first iterate of $`u`$ in the DSM and the average is performed over all possible realizations of the stochastic process. Simple algebra yields

$$\mathrm{\Lambda }_0(u)=\mathrm{ln}aau\mathrm{ln}(au).$$
(43)
In the limit $`u0`$, $`\mathrm{\Lambda }_0(u)`$ reduces correctly to the true Lyapunov exponent $`\mathrm{log}a`$. If $`a`$ is not too small ($`a>1/e`$, where $`e`$ is the Neper number), $`\mathrm{\Lambda }_0(u)`$ has a maximum value for an intermediate value of $`u`$ and, more interestingly, the maximum value of $`\mathrm{\Lambda }`$ is larger than 0, if $`a>\mathrm{exp}(1/e)=0.692..`$. Therefore, this result is suggestive of a phase transition from a regime (small $`a`$) where perturbations of all sizes decrease, to a regime where sufficiently large pertubations expand (and can, in principle, self-sustain). However, we already know that this conclusion is incorrect: the reason is precisely that $`\mathrm{\Lambda }_0(u)`$ is an average quantity, and fluctuations must be taken into account. In fact, a different scenario arises, if we define the contraction rate by taking the logarithm of the average expansion factor. In this case, we obtain

$$\mathrm{\Lambda }_1(u)=\mathrm{ln}a+\mathrm{ln}(2au),$$
(44)
an expression that can be larger than zero even at $`u=0`$, where it actually attains its maximum value! In practice, one has to be very careful in drawing meaningful conclusions from any expression of the finite-size Lyapunov exponent. In particular, since $`u`$ can never be larger than 1, while it can become arbitrarily small, the typical negative contraction rate operating at small $`u`$-values eventually wins, making always $`u=0`$ the only stable fixed point.

## 4 The spatially extended case

The past study of stable chaos has revealed that no qualitative difference exists, as long as strictly finite chains are considered, with respect to standard stable systems. After an exponentially long transient, a periodic behaviour is always attained. As it was remarked in Section II, a similar scenario is indeed exhibited by our stochastic models.

Eq. (29) shows that in the absence of any coupling, there is a finite probability $`\stackrel{~}{P}(u(0))`$ for an initial perturbation $`u(0)`$ to be contracted for an infinite amount of time steps and thus of being effectively absorbed. It is worth defining here

$$u_M(t)=\underset{i}{\mathrm{max}}u_i(t).$$
(45)
In particular, this allows us applying the above reasoning to the coupled case, with $`u_M(0)`$ playing the role of $`u(0)`$; the probability $`\stackrel{~}{P}`$ that an arbitrary, spatially extended, perturbation is contracted forever on every site $`i`$ satisfies the inequality

$$\stackrel{~}{P}(\{u_i\}_{i=1,\mathrm{},L}\left[\underset{n=1}{\overset{\mathrm{}}{}}(1u_M(0)a^n)\right]^L=\left[\stackrel{~}{P}(u_M(0))\right]^L,$$
(46)
where $`\stackrel{~}{P}(u_M(0))^L`$ is small but finite for every finite lattice length $`L`$. Thus, the average time $`\tau `$ needed for any perturbation to die out (i.e.to be contracted below some arbitrarily small thresold value) is finite and does not grow faster than exponentially with the system size,

$$\tau \left[\stackrel{~}{P}(u_M(0))\right]^L.$$
(47)
In the active, i.e. chaotic, phase where perturbations propagate with a finite velocity, $`\tau `$ grows exactly exponentially (with the system size), since the only way for a perturbation to die out is to be contracted in all sites, while in the inactive phase, the latter is only a sufficent condition and not a necessary one, and perturbations die out on significantly shorter time scales.

At variance with the 0-dimensional case, direct numerical simulations of the DSM in a 1-d lattice reveal the existence of a regime where finite perturbations self-sustain. This is clearly shown in Fig. 6, where we have reported the space-time and ensenble average value $`\stackrel{~}{m}`$ of the perturbation for different values of the contraction rate. Above $`a=0.60(5)`$, nonzero amplitudes are actually observed. It is, therefore, crucial to understand the reason why the spatial interactions can stabilize finite perturbations in spite of the diffusive nature of the coupling.

The natural extension of what we have learnt in zero-dimension consists in looking at the joint probability distribution $`R(t,u_1,u_2,\mathrm{},u_n,\mathrm{})`$ over the whole lattice. If the system is sufficiently above the transition to stable chaos, it seems reasonable, in a first approximation, to neglect spatial correlations. This is certainly incorrect for those sites that are close to propagating fronts but the fraction of such lattice sites is definitely negligible. Accordingly, we approximate the joint probability distribution as a product of single-site probabilities $`Q(t,u)`$,

$$R(t,u_1,u_2,\mathrm{},u_n,\mathrm{})Q(t,u_1)Q(t,u_2)\mathrm{}Q(t,u_n)\mathrm{}.$$
(48)
Within this approximation, the single-site probability distribution corresponding to the stochastic dynamics (28) satisfies the following Frobenius-Perron equation,

$`Q(t+1,u)`$ $`=`$ $`g(u){\displaystyle _0^{\mathrm{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}dv_iQ(t,v_i)\delta \left(u{\displaystyle \frac{a}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}v_i\right)`$
$`+a\mathrm{\Delta }{\displaystyle _0^{\mathrm{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}dv_iQ(t,v_i)\delta \left(u{\displaystyle \frac{1}{N\mathrm{\Delta }}}{\displaystyle \underset{i=1}{\overset{N}{}}}v_i\right)`$
$$Q(t+1,1)=_{a\mathrm{\Delta }}^{\mathrm{}}𝑑uu_0^{\mathrm{}}\underset{i=1}{\overset{N}{}}dv_iQ(t,v_i)\delta \left(u\frac{a}{N}\underset{i=1}{\overset{N}{}}v_i\right),$$
(49)
where $`\delta `$ is the Dirac’s distribution, $`v_i`$ is the amplitude of the perturbation in the $`i`$th neighbouring site, $`N`$ is the number of democratically coupled sites (for later convenience we leave $`N`$ unspecified - notice that $`N=3`$ in a 1-d lattice with nearest-neighbour coupling) and

$$g(u)=\{\begin{array}{cc}1u\hfill & \hfill \mathrm{\Delta }<u<1\\ 1a\mathrm{\Delta }\hfill & \hfill 0u\mathrm{\Delta }\end{array}.$$

It is easy to verify that the support of the single-site probability distribution $`Q(t,u)`$ remains confined to the unit interval, provided that this holds true for the initial condition as well. Furthermore, due to the factorization hypothesis, the space and ensemble average $`m(t)`$ defined in Section II coincides with the simplest ensemble average, i.e. the mean value of $`Q(t,u)`$.

In Fig. 7 we have plotted the single site probability distributions obtained by directly iterating the stochastic model and the approximate Frobenius-Perron equation (49) (solid and dashed line, respectively). The reasonable overlap confirms the validity of the factorization hypothesis (at least away from the critical region). The mean value, equal to $`0.773`$ in former case, compares with $`0.770`$ in the latter one, while the variances are respectively equal to $`0.063`$ and $`0.061`$. Such small differences are due to the different behavior of the probability distributions for small values of $`u`$.

For $`N=1`$ (no coupling), Eq. (49) corresponds to the 0-dimensional dynamics discussed in the previous section, and the evolution equation is exact. It reduces to the linear equation,

$`Q(t+1,u)`$ $`=`$ $`g(u){\displaystyle \frac{1}{a}}Q(t,{\displaystyle \frac{1}{a}}u)+a\mathrm{\Delta }^2Q(t,a\mathrm{\Delta }u)`$
$`Q(t+1,1)`$ $`=`$ $`a{\displaystyle _\mathrm{\Delta }^{\mathrm{}}}𝑑uuQ(t,u),`$ (50)
The only fixed point of this equation is $`Q(t,u)=\delta (u)`$, i.e. the absorbing state. The multifractal spectrum discussed in Sec. III is nothing but a sophisticate characterization of the convergence towards such a fixed point.

On the opposite side of the 0-dimensional limit, there is the mean-field approximation that corresponds to the limit $`N\mathrm{}`$. In this limit, statistical fluctuations vanish and the dynamics reduces to the evolution of the mean value $`m(t)`$, that reads

$$m(t+1)=\{\begin{array}{cc}m(t)[2aa^2m(t)]\hfill & \hfill \text{if}m(t)>\mathrm{\Delta }\\ m(t)[2aa^2\mathrm{\Delta }]\hfill & \hfill \text{if}m(t)\mathrm{\Delta }.\end{array}$$
(51)
For $`a<a_c=(1\sqrt{1\mathrm{\Delta }})/\mathrm{\Delta }`$, Eq. (51) displays the stable fixed point $`m_1=0`$. This is the same regime found in 0-dimension and corresponds to the eventual absorption of any initial finite difference. Increasing $`a`$ above $`a_c`$, the system undergoes a bifurcation: $`m_1`$ becomes unstable and a second (stable) fixed point $`m_2=(2a1)/a^2`$ appears. In the discontinuous limit ($`\mathrm{\Delta }0`$) $`a_c1/2`$. In Fig. 6 the predictions of the mean-field approach are compared with the results of direct simulations: we see that, in spite of the approximations, there are no severe differences and the critical point is underestimated by approximately 17%.

It is rather instructive to notice that the predictions of the mean-field analysis do coincide with the finite-size Lyapunov exponent $`\mathrm{\Lambda }_1`$ (for the sake of simplicity, we limit ourselves to consider the DSM). Since the mean field approximation reduces the CSM dynamics to the evolution of a single variable $`m(t)`$, the comparison can be performed by interpreting $`\mathrm{ln}[m(t+1)/m(t)]`$ as $`\mathrm{\Lambda }_1`$ and recalling that Eq. (44) has been derived for the DSM only. Even though Eq. (44) can be obtained by Eq. (51), there is an important difference between the consequences of the two results. As we discussed in Section III, fluctuations keep $`u^{}=0`$ stable for every value of $`a`$, while in the case of the mean field analysis, the lack of any fluctuation, due to the formally infinite number $`N`$ of neighbours, implies that the fixed point $`u=0`$ is truly unstable when $`a>1/2`$.

The absence of fluctuations in the mean field limit implies that all definitions of the finite-size Lyapunov exponent are equivalent. Therefore, we observe the same scenario previously observed for the standard maximum Lyapunov exponent: the diffusive coupling shifts the Lyapunov exponent from the average value of the logarithm of the multiplier (the so-called quenched average holding for the single map) towards the logarithm of the average multiplier (annealed average, predicted by the mean field analysis). The important consequence of this shift is that, in the present context, it can change the stability of the $`u=0`$ solution leading to the onset of the chaotic phase.

Anyway, one should not forget that the mean-field anlaysis provides an approximate solution. For $`N`$ finite and strictly larger than 1, Eq. (49) defines a non trivial evolution operator $`𝒬`$ in a functional space. In practice, one can expand the evolution equation into an infinite set of equations for, e.g. the momenta $`M_k`$ of $`Q`$. An approximate solution can thus be found by either suitably truncating the hierarchy of equations or introducing a closure Ansatz. In App. B we parametrize the probability distribution as the sum of a $`\delta `$-distribution and a Gaussian. This allows us deriving three evolution equations for the DSM.

In both cases that we have investigated ($`N=2`$, 3), we find a scenario similar to the one predicted by the mean field analysis. There exists a critical value $`a_c`$ (equal to 0.548 for $`N=2`$ and to 0.536 for $`N=3`$) below which the dynamics is characterized by the stable fixed point $`A=1`$, $`v=0`$, $`V=0`$ (corresponding to the absorbing state $`Q(t,v)=\delta (v)`$) and above which the previous solution becomes unstable, giving rise to a stable nontrivial solution. The dependence of $`m`$ on $`a`$ reported in Fig. 6 (see the dashed curve in Fig. 6) indicates that the critical value predicted by this analytic approach improves the mean-field estimate, but the growth of $`m`$ above threshold is not as good as one would like. In fact, there is a qualitative difference with the mean-field approach: a further bifurcation (at $`a=\overline{a}=0.855`$ for $`N=2`$ and $`\overline{a}=0.869`$ for $`N=3`$), where the new solution destabilizes too. Such a bifurcation and the slow growth of $`\stackrel{~}{m}`$ with $`a`$ are both consequences of a defect of the approximation: the support of the Gaussian extends out of the unit interval. This unphysical property becomes increasingly important as soon as the average amplitude of $`u`$ is comparable with 1.

Indeed, a better agreement with the direct simulations is obtained by iterating numerically the Frobenius-Perron equation (see Fig. 6). Neither simulations performed with $`\mathrm{\Delta }=0`$ nor with $`\mathrm{\Delta }=0.01`$ reveal the second bifurcation found with the Gaussian approximation, confirming that it is an artifact of the approximation. In the continuous case, the bifurcation occurs at the critical value $`a_c=0.591\mathrm{}`$ (for $`N=2`$) and $`a_c=0.568\mathrm{}`$ (for $`N=3`$), to be compared with the mean field prediction $`a_c=0.501\mathrm{}`$. In the discontinuous case we find $`a_c=0.585\mathrm{}`$ (for $`N=2`$) and $`a_c=0.567\mathrm{}`$ (for $`N=3`$), to be compared with the mean field prediction $`a_c=0.5`$.

From the data reported in Fig. 6 for $`N=3`$ we see that the factorization hyptohesis reproduces fairly well the behaviour of the full stochastic model everywhere except for the transition region. This is not unexpected as it is well known that the correlation length diverges in the critical region.

## 5 Open problems and conclusions

In this paper we have shown that a simple stochastic model, specifically designed to simulate a different response to finite and infinitesimal perturbations, is able to capture the key features of irregular behaviour in linearly stable systems. In particular, we have seen that replacing the sequence of jumps generated by the CML dynamics with a genuine stochastic process allows for a faithful reconstruction of the front propagation. The main theoretical advantage of the stochastic model is the disentanglement between the generation of a pseudo-random pattern and the evolution of perturbations. In reality the two issues are interlaced: their separation has allowed us clarifying under which conditions (finite) perturbations can be effectively sustained throughout an infinite lattice. In particular, a full consistency exists between the CML and the stochastic model in the chaotic regime, since we can state that the amplification of finite perturbations contributes to sustain an irregular regime.

On the other hand, the transition to the ordered phase observed in the CSM/DSM models does not reproduce the analogous behaviour displayed by the CML. In this latter case, it was observed that the critical region is not point-like, but rather extended to a what has been called “fuzzy region”, where ordered and chaotic dynamics alternate in a quite irregular manner . The reason for the difference is that in the CML model, the absence of local chaos makes the sequence of multipliers increasingly less random in the transition region. In the stochastic model, instead, the randomness of multipliers is always assumed a priori. In spite of such a difference, it is nevertheless instructive to notice that a transition persists without modifying the stochasticity in the real space.

A more precise analogy for the transition investigated in the previous section is provided by the correspondence with the problem of synchronization in the presence of external noise. In fact, in this latter context, the noise represents the (unvaoidable) source of stachasticity in the synchronous as well as in the asynchronous regime. Since a recent numerical study of the synchronization transition in linearly stable system has suggested that it belongs to the universality class of Directed Percolation (DP) , it is tempting to verify whether the same holds true in our stochastic models.

As already mentioned, even though the front velocity $`v_F`$ is a good order parameter to characterize the transition, it is quite difficult to obtain a reliable estimate of the critical value of the control parameter $`a_c`$ from the vanishing of $`v_F`$. In fact, finite-size and transient effects combined with the existence of wild fluctuations prevent a careful analysis. A more efficient method amounts to measuring the dependence of the so-called absorption time $`\tau `$ on the system size $`L`$. This is defined as the time required for the space averaged perturbation $`\overline{u}(t)`$ to become smaller than some very small, but finite threshold $`\mathrm{\Gamma }`$. Fluctuations of $`\tau `$ can be efficiently reduced by averarging over a sufficiently large ensemble of initial conditions. In the active phase, $`\tau `$ is expected to diverge exponentially with $`L`$ (stable chaotic regime), while in the absorbing phase, it should depend at most logarithmically on $`L`$. Only at the critical point, $`\tau `$ exhibits a power law dependence

$$\tau (L,a_c)L^z,$$
(52)
where $`z=\nu _l/\nu _{}`$ is the so called dynamical exponent. We have performed numerical simulations for both the CSM (with $`\mathrm{\Delta }=0.01`$) and DSM, averaging over 3000 realizations of the stochastic process and over randomly sampled initial conditions. In the CSM we find $`a_c=0.6055\mathrm{}`$, with $`z=1.58\pm 0.02`$; in the DSM we obtain $`a_c=0.6065\mathrm{}`$ and $`z=1.56\pm 0.06`$. The errors have been estimated as the maximum deviation from linearity in the log-log plot that has been used for extracting the scaling law (52) (see, for instance, Fig. 8, where $`\tau `$ has been plotted versus $`L`$ for $`\mathrm{\Delta }=0.01`$ and different values of $`a`$). These results agree with the most accurate numerical estimates of the DP value, $`z=1.5807`$ .

We also measured the critical exponent $`\delta `$ associated with the temporal decay of the density of active sites $`\rho (t)`$, i.e. those sites where $`u_i(t)>\mathrm{\Gamma }`$: at the critical point $`\rho (t)`$ is characterized by the scaling law

$$\rho (t,a_c)t^\delta ,$$
(53)
where $`\delta =\beta /\nu _l`$ (as usual $`\nu _l`$, $`\nu _{}`$ and $`\beta `$ are the critical exponents respectively associated with space and time correlation lengths and with the order parameter). By averaging over 3000 realizations and choosing a sufficiently large value of $`L`$ to get rid of finite-size corrections we have found $`\delta =0.150\pm 0.01`$ for the CSM and $`0.155\pm 0.005`$ for the DSM, to be compared with the DP value $`\delta =0.1595`$ .

Altogether, our simulations support the hypothesis that the transition belongs to the same universality class as DP. This may look as an almost trivial result, since local spreading is the only mechanism for the propagation of perturbations (or, in different languages, active sites, infections). However, the whole problem is definitely more subtle, as an absorbing state cannot be identified so clearly. In fact, we have already seen that the determination of the absorption time requires to fix a somehow arbitrary threshold $`\mathrm{\Gamma }`$, and the same is true for the computation of the active sites at a given time. Even though we have found that our results are independent of the choice of $`\mathrm{\Gamma }`$ (provided that it is small enough), this may appear as a numerical trick: there is nothing like a true threshold, since no matter how small is a perturbation, there is always some finite probability that it gives rise to a burst: this is contrary to the existence of a truly absorbing state. On the other hand, the argument presented in the beginning of the previous section to convince the reader that in a finite chain any perturbation eventually dies out, confirms the existence of an absorbing state: since the smaller is $`u_M`$ the more likely is that the perturbation keeps being absorbed, any perturbation $`u_i(t)`$ has a finite probability to enter an infinite “contraction loop” in which every site $`i`$ is monotonously contracted for any time larger than $`t`$.

It would be nice to put our qualitative arguments on a more rigorous basis, by definining a suitable finite-size Lyapunov exponent that is negative below some threshold to indicate that $`u^{}=0`$ is a truly absorbing state. However, it is not clear whether this could be accomplished, since we know that in the chaotic phase, $`u^{}=0`$ should be at the same time “macroscopically” unstable since perturbations eventually drive the system towards the only stable state and “microscopically” stable to mean that small enough perturbations have to be absorbed. In the future, we hope to be able to clarify whether it is possible to define an indicator that contains both messages.

RL and AP wish to thank R. Kapral for early discussions about a meaningful definition of the stochastic model. A. Pikovsky, V. Ahlers and A. Torcini are acknowledged for useful exchanges of ideas about the characterization of the transition. A profitable discussion on the 0-dimensional model has been carried on with Y. Elskens. Part of this work has been completed thanks to the financial support of the NATO contract CRG.973054. We also thank I.S.I. in Torino where part of this work was performed.

## Appendix A Multifractal distribution

In this appendix we report the analytical calculation of the (constant) slope of the multifractal distribution (31) for the zero dimensional DSM. We shall prove that $`H^{}(\mathrm{\Lambda },u(0))`$ (the prime denotes derivative with respect to $`\mathrm{\Lambda }`$) is independent of $`\mathrm{\Lambda }`$ and increases monotonously with the contraction rate $`a`$ from -1 (for $`a0`$) to 0 (for $`a1`$).
First of all, note that the condition $`0nn_0+t`$, together with Eq. (36), implies that the support of $`H(\mathrm{\Lambda })`$ is confined in the interval $`[\mathrm{ln}a,0]`$, since the ratio $`\rho =n/t`$ can take values between 0 and 1 and in the $`t\mathrm{}`$ limit $`n_0/t`$ vanishes.

We are now interested in the $`\mathrm{\Lambda }>\mathrm{ln}a`$ case (i.e. the one in which the system received at least one kick), where, as it was stated in Section III, the actual size of the perturbation $`u`$ at time $`t`$ is unambiguously determined by the time $`n<t`$ elapsed since the last kick (see Eq. (36)) and thus, for any finite $`t`$, both $`u`$ and $`\mathrm{\Lambda }`$ can assume only a discrete set of values labelled by $`n`$.

Let us now denote the discrete $`\mathrm{\Lambda }`$-derivative of a generic function $`f(\mathrm{\Lambda })`$ as

$$D_\mathrm{\Lambda }(f)=\frac{f(\mathrm{\Lambda }+\mathrm{\Delta }\mathrm{\Lambda })g(\mathrm{\Lambda })}{\mathrm{\Delta }\mathrm{\Lambda }}.$$
(54)
From Eq. (31), we can approximate the derivative of the multifractal distribution $`H`$ as

$$H^{}(\mathrm{\Lambda },u(0))=\underset{t\mathrm{}}{lim}\frac{D_\mathrm{\Lambda }(\mathrm{ln}𝒫(\mathrm{\Lambda },t,u(0))}{t},$$
(55)
where $`\mathrm{\Delta }\mathrm{\Lambda }`$ is naturally fixed by Eq. (36) and the discrete character of $`n=\mathrm{ln}u/\mathrm{ln}a`$,

$$\mathrm{\Delta }\mathrm{\Lambda }=\mathrm{\Lambda }_{\tau +1}(t,u(0))\mathrm{\Lambda }_\tau (t,u(0))=\frac{\mathrm{ln}a}{t}.$$
(56)
As, for $`t\mathrm{}`$, $`\mathrm{\Delta }\mathrm{\Lambda }`$ goes to $`0`$, Eq. (55) becomes asymptotically exact.

Morover, Eq. (36) allows one switching to the $`Q`$ rapresentation of probabilities writing

$$H^{}(\mathrm{\Lambda },u(0))=\frac{1}{\mathrm{ln}a}\underset{t\mathrm{}}{lim}\mathrm{ln}\left[\frac{𝒫(\mathrm{\Lambda }+\mathrm{\Delta }\mathrm{\Lambda },t,u(0))}{𝒫(\mathrm{\Lambda },t,u(0))}\right]=\frac{1}{\mathrm{ln}a}\underset{t\mathrm{}}{lim}\mathrm{ln}\left[\frac{Q(t,n+1)}{Q(t,n)}\right]|_{n=n_0+\frac{t\mathrm{\Lambda }}{\mathrm{ln}a}}$$
(57)
(remember also that $`u(0)=a^{n_0}`$). Making use of Eqs. (37) and (29), we obtain

$$H^{}(\mathrm{\Lambda },u(0))=\frac{1}{\mathrm{ln}a}\left[\mathrm{ln}\left(\frac{G(tn1,n_0)}{G(tn,n_0)}\right)+\mathrm{ln}\left(1a^{n+1}\right)\right]|_{n=n_0+\frac{t\mathrm{\Lambda }}{\mathrm{ln}a}}.$$
(58)
In the limit $`t\mathrm{}`$, the last term in the r.h.s. vanishes, provided $`n`$ diverges too (i.e. $`\mathrm{\Lambda }<0`$). We are thus left with the following equality, which holds true for any $`\mathrm{\Lambda }`$ in the open interval $`(\mathrm{ln},0)`$,

$$H^{}(\mathrm{\Lambda },u(0))=\frac{\mathrm{ln}\eta (a,n_0)}{\mathrm{ln}a},$$
(59)
where

$$\eta (a,n_0)=\underset{t\mathrm{}}{lim}\frac{G(t+1,n_0)}{G(t,n_0)},$$
(60)
the $`n`$-dependence has been eliminated by shifting the time origin and the dependence on the parameter $`a`$ has been made explicit. From this equation, we see that the problem of determining the slope of the multifractal spectrum is equivalent to an eigenvalue problem. In fact, we can formally write Eq. (38) as

$$\stackrel{}{v}(t,n_0)=M\stackrel{}{v}(t1,n_0),$$
(61)
where $`\stackrel{}{v}(t,n_0)`$ is the infinite-dimensional vector

$$\stackrel{}{v}(t,n_0)(G(t,n_0),G(t1,n_0),G(t2,n_0),\mathrm{},)$$
(62)
and $`M`$ is a linear infinite-dimensional operator. Therefore, $`\eta (a)`$ is nothing but the maximum eigenvalue of the operator $`M`$.

From the explicit expression of the recursive equation for $`G(t,n_0)`$ (see Eq. (38))

$$G(t,n_0)=G(t1,n_0)P(0,1)a+G(t2,n_0)P(1,1)a^2+\mathrm{}+a^{n_0}\underset{k=1}{\overset{t1}{}}(1a^{k+n_0})a^t,$$
(63)
we see that, since both $`G(t,n_0)`$ and $`P(t,1)`$ are not larger than 1 (they are probabilities), the sum in the r.h.s. is bounded from above by the sum of the first $`t`$ powers of $`a`$. In the limit $`t\mathrm{}`$, the function $`G(t,n_0)`$ converges exponentially fast, so that we are able to truncate Eq. (63) to order $`k`$ with an arbitrary precision. This truncation makes the problem numerically solvable, as the operator $`M`$ can be approximated by the finite dimensional matrix $`M_k`$ (here and below, the subscript $`k`$ stands for $`k`$th order approximation) that we report here below,

$$M_k=\left(\begin{array}{c}\alpha _1\\ 1\\ 0\\ 0\\ \mathrm{}\\ 0\end{array}\begin{array}{c}\alpha _2\\ 0\\ 1\\ 0\\ \mathrm{}\\ 0\end{array}\begin{array}{c}\alpha _3\\ 0\\ 0\\ 1\\ \mathrm{}\\ 0\end{array}\begin{array}{c}\mathrm{}\\ \mathrm{}\\ \mathrm{}\\ \mathrm{}\\ \mathrm{}\\ \mathrm{}\end{array}\begin{array}{c}\alpha _k\\ 0\\ 0\\ 0\\ \mathrm{}\\ 1\end{array}\right),$$
(64)
where $`\alpha _j=a^jP(j1,1)`$. Numerical estimates of $`\eta _k(a)`$ indicate that $`k=10`$ suffices to attain a good convergence in the whole range of $`a`$ values between 0 and 1. For $`k=2`$ the eigenvalue $`\eta _2(a)`$ can be computed analytically yielding,

$$H_2^{}(\mathrm{\Lambda },u(0))=1\frac{1}{\mathrm{ln}a}\mathrm{ln}\frac{1+\sqrt{54a}}{2},$$
(65)
which goes monotonously from -1 ($`a0`$) to $`0^{}`$ ($`a1`$).

## Appendix B An approximate solution of the 1-dimensional problem

The Frobenius-Perron operator defined by Eq. (49) is a functional equation whose solution requires projecting it onto a finite dimensional space either via some truncation or a suitable closure hypothesis. Looking at a typical shape of the stationary probability distribution in Fig. 7, we can see that the multi-peaked structure slows down the convergence of an expansion in either moments or cumulants. Accordingly, we have preferred to approximate the probability distribution as the sum of a Dirac’s $`\delta `$ distribution centered in 1 (the reinjection point) plus a Gaussian distribution, centered around a point to be determined self-consistently,

$$Q(t,v)=A_t\delta (v1)+\frac{1A_t}{\sqrt{2\pi V_t}}\mathrm{exp}\left[\frac{(v\overline{v}_t)^2}{2V_t}\right].$$
(66)
We can see that $`Q(t,v)`$ is parametrized by three quantities: $`A_t`$, the probability of the $`\delta `$ component, $`V_t`$ the variance of the Gaussian, and $`\overline{v}_t`$, its average value. The mean value $`m(t)`$ of $`Q(t,u)`$ is therefore equal to

$$m(t)=1A_t+A_t\overline{v}_t.$$
(67)
It should be noted that our approximation is formally “unphysical” since the support of any Gaussian function is not restricted to the unit interval, but we expect our Ansatz to be reasonably correct as long as the probability to be out of the unit interval is small enough.

Entering the above definition of $`Q(t,v)`$ into Eq. (49) and computing separately the new weight of the $`\delta `$ component, the average and the variance of the Gaussian, we obtain three evolution equations. In the simplest nontrivial case, $`N=2`$, we have

$`A_{t+1}`$ $`=`$ $`a\left[A_t(1\overline{v}_t)+\overline{v}_t\right]`$
$`\overline{v}_{t+1}`$ $`=`$ $`v_t^2\left[2A_t\right]{\displaystyle \frac{a^2}{2}}+\overline{v}_ta\left[1aA_t\right]a\left[1+{\displaystyle \frac{a}{2}}(V_t2+A_t){\displaystyle \frac{1a}{1A_t}}\right]`$ (68)
$`V_{t+1}`$ $`=`$ $`v_t^3\left[A_t\left(a{\displaystyle \frac{3}{4}}a^2\right)(1a)^2\right]a+\overline{v}_t^2\left[A_t(2a){\displaystyle \frac{3}{4}}a^2+a^2a+{\displaystyle \frac{1a}{1A_t}}\right]+`$
$`\overline{v}_t\left[A_t(1V_t){\displaystyle \frac{3}{4}}a^3{\displaystyle \frac{3}{2}}a^3V_t+a^2(V_t2)+2a{\displaystyle \frac{2a(1a)}{1A_t}}\right]`$
$`A_t[a(V_t1)+2]{\displaystyle \frac{a^2}{4}}+a^3+{\displaystyle \frac{a^2}{2}}(V_t2)+{\displaystyle \frac{a^2a^3}{1A_t}}.`$
For $`N=3`$, the case that we have investigated in detail as it corresponds to a 1-d lattice of democratically coupled maps, the equations read

$`A_{t+1}`$ $`=`$ $`a\left[A_t(1\overline{v}_t)+\overline{v}_t\right]`$
$`\overline{v}_{t+1}`$ $`=`$ $`\overline{v}_t^2[3A_t]{\displaystyle \frac{a^2}{3}}+\overline{v}_t(a{\displaystyle \frac{4}{3}}a^2A_t)a\left[1+{\displaystyle \frac{a}{3}}(V_t32A_t){\displaystyle \frac{1a}{1A_t}}\right]`$
$`V_{t+1}`$ $`=`$ $`\overline{v}_t^3\left[{\displaystyle \frac{2}{9}}A_t^2a^2+{\displaystyle \frac{2}{9}}A_t(5a2)a(1a)^2\right]a+`$
$`\overline{v}_t^2\left[{\displaystyle \frac{2}{3}}A_t^2a^3{\displaystyle \frac{2}{3}}A_t(2a3)a^2+a^2a+{\displaystyle \frac{1a}{1A_t}}\right]`$
$`{\displaystyle \frac{2\overline{v}_t}{3}}\left[A_t^2a^3A_t(3+2a+V_ta)2a\left(a^21+{\displaystyle \frac{V_t}{2}}(23a)a3{\displaystyle \frac{1a}{1A_t}}\right)\right]+`$
$`A_t^2\left({\displaystyle \frac{2}{9}}a^3+a^2a\right)A_t\left[{\displaystyle \frac{2}{3}}a^3(V_t+{\displaystyle \frac{4}{3}}){\displaystyle \frac{1}{3}}a^2+a\right]+`$
$`{\displaystyle \frac{8}{9}}a^3+V_t{\displaystyle \frac{a^2}{3}}a+{\displaystyle \frac{2(aa^3)}{1A_t}}.`$
In spite of the great simplifications involved in the derivation of both sets of equations, it is still impossible to obtain an analytic expression even for the critical point $`a_c`$. The results of the numerical solution of Eq. (B) and (B) are discussed in the text (see Sec. IV).


format: text

text:
# Abstract

### Abstract

We discuss the $`\phi ^4`$ and $`\phi ^6`$ theory defined in a flat $`D`$-dimensional space-time. We assume that the system is in equilibrium with a thermal bath at temperature $`\beta ^1`$. To obtain non-perturbative result, the $`1/N`$ expansion is used. The method of the composite operator (CJT) for summing a large set of Feynman graphs, is developed for the finite temperature system. The ressumed effective potential and the analysis of the $`D=3`$ and $`D=4`$ cases are given.

## 1 Introduction

The conventional perturbation theory in the coupling constant or in $`\mathrm{}`$ i.e., the loop expansion can only be used for the study of small quantum corrections to classical results. When discussing quantum mechanical effects to any given order in such an expansion, one is not usually able to justify the neglect of yet higher order. In other words, for theories with a large $`N`$ dimensional internal symmetry group, there exist another perturbation scheme, the $`1/N`$ expansion, which circumvents this criticism. Each term in the $`1/N`$ expansion contains an infinite subset of terms of the loop expansion. The $`1/N`$ expansion has the nice property that the leading-order quantum corrections are of the same order as the classical quantities. Consequently, the leading order which adequately characterizes the theory in the large $`N`$ limit preserves much of the nonlinear structure of the full theory. In the next section we derive the effective action to leading order in $`1/N`$ in $`D`$-dimensional space-time and consequently the effective potential. It is known that, in $`D>4`$, such theories with $`\phi ^4`$ interaction are in fact free field theory, while in $`D<4`$ they have a non-trivial continuum limit as an interacting field theory. For $`D=3`$ it has been shown that, in the large $`N`$ limit, the $`\phi ^6`$ theory has a UV fixed point and therefore must have a second IR fixed point . At least for large $`N`$ the $`(\phi ^6)_{D=3}`$ theory is known to be qualitatively different from $`(\phi ^4)_{D=4}`$ theory.

In a previous work  by use of the composite operator formalism and, we re-examinate the behavior at finite temperature of the $`O(N)`$-symmetric $`\lambda \phi ^4`$ model in a generic D-dimensional Euclidean space. In the cases $`D=3`$ and $`D=4`$, an analysis of the thermal behavior of the renormalized squared mass and coupling constant are done for all temperatures. It results that the thermal renormalized squared mass is positive and increases monotonically with the temperature. It is interesting to stress that the behavior of the thermal coupling constant is quite different in odd or even dimensional space. In $`D=3`$, the thermal coupling constant decreases up to a minimum value different from zero and then grows up monotonically as the temperature increases. In the case $`D=4`$, it is found that the thermal renormalized coupling constant tends in the high temperature limit to a constant asymptotic value. Also for general D-dimensional Euclidean space, we are able to obtain a formula for the critical temperature of the second order phase transition. This formula agrees with previous known results at $`D=3`$ and $`D=4`$ .

It is well known that the introduction of the $`\phi ^6`$ term generated a rich phase diagram, with the possibility of second order, first order phase transitions or even both transitions occurring simultaneously. This situation defines the tricritical phenomenon. Some systems such antiferromagnets in the presence of a strong external field or the $`He^3He^4`$ mixture exhibits such behavior. In a previous paper the massive $`(\phi ^6)_{D=3}`$ model was analyzed at finite temperature at the two-loops approximation. We demonstrate the existence of the tricritical point . A natural extension of this paper was done in ref. . In this paper we proved the existence of the tricritical point using a non-pertubative approach. This was done using the CJT formalism i.e. the composite operator formalism .

Here we continuous to study the composite operator method. Still studying the $`\phi ^6`$ theory in the large N expansion, the effective potential at finite temperature is calculated. The organization of the paper is the following. In section II we derive the effective potential using the composite operator (CJT) formalism. In section III the thermal effective potential is found for a D-dimensional generic space. Conclusions are given in section IV. In this paper we use $`\mathrm{}=c=k_b=1`$.

## 2 The effective potential (The CJT formalism)

We are interested here in the most general renormalizable scalar field model $`\lambda \phi ^4+\eta \phi ^6`$ possessing an internal symmetry $`O(N)`$, in a generic $`D`$-dimensional space-time. Of course, for $`D=4`$ this theory is non-renormalizable. In this case to ensure renomalisability we must make $`\eta =0`$. Let us define at the beginning the field in a generic D-dimensional space-time. For simplicity we will call this theory a $`\phi ^6`$ model.

Using the method of composite operator developed by Cornwall, Jackiw and Tomboulis , Townsend derived the effective potential of $`\phi ^6`$ theory in the $`1/N`$ expansion for $`D=3`$ at zero temperature . This author proved that $`1/N`$ expansion is consistent for $`\phi ^6`$ to leading order.

The Lagrangian density of the $`O(N)`$ symmetry $`\phi ^6`$ theory is :

$$(\phi )=\frac{1}{2}(_\mu \phi )^2\frac{1}{2}m_0^2\phi ^2\frac{\lambda _0}{4!N}\phi ^4\frac{\eta _0}{6!N^2}\phi ^6,$$
(1)
where the quantum field is an $`N`$-component vector in the $`N`$-dimensional internal symmetry space. For definiteness, we work at zero-temperature; however, the finite temperature generalizations can be easily obtained . We are interested in the effective action $`\mathrm{\Gamma }(\varphi )`$ which governs the behavior of the expectation values $`\phi _a(x)`$ of the quantum field where $`\varphi `$ is given by

$$\varphi (x)\frac{\delta W(J)}{\delta J(x)}=<0|\phi (x)|0>,$$
(2)
where $`W(J)`$ is the generating functional for connected Green’s functions.

$`\mathrm{\Gamma }(\varphi )`$ can be shown to be the sum of one-particle irreducible (1PI) Feynman graphs with a factor $`\varphi _a(x)`$ on the external line. We make use of the formalism of composite operator which reduces the problem to summing two particle irreducible (2PI) Feynman graphs by defining a generalized effective action $`\mathrm{\Gamma }(\varphi ,G)`$ which is a functional not only of $`\varphi _a(x)`$, but also of the expectation values $`G_{ab}(x,y)`$ of the time ordered product of quantum fields $`<0|T(\phi (x)\phi (y))|0>`$, i.e.

$$\mathrm{\Gamma }(\varphi ,G)=I(\varphi )+\frac{i}{2}TrLnG^1+\frac{i}{2}TrD^1(\varphi )G+\mathrm{\Gamma }_2(\varphi ,G)+\mathrm{},$$
(3)
where $`I(\varphi )=𝑑x^D(\varphi )`$, $`G`$ and $`D`$ are matrices in both the functional and the internal space whose elements are $`G_{ab}(x,y)`$, $`D_{ab}(\varphi ;x,y)`$ respectively and $`D`$ is defined by

$$iD^1=\frac{\delta ^2I(\varphi )}{\delta \varphi (x)\delta \varphi (y)}.$$
(4)
The quantity $`\mathrm{\Gamma }_2(\varphi ,G)`$ is computed as follows. In the classical action $`I(\phi )`$ we have to shift the field $`\phi `$ by $`\varphi `$. The new action $`I(\phi +\varphi )`$ possesses terms cubic and higher in $`\phi `$. This define an interaction part $`I_{int}(\phi ,\varphi )`$ where the vertices depend on $`\varphi `$. $`\mathrm{\Gamma }_2(\varphi ,G)`$ is given by sum of all (2PI) vacuum graphs in a theory with vertices determined by $`I_{int}(\phi ,\varphi )`$ and the propagators set equal to $`G(x,y)`$. The trace and logarithm in eq.(3) are functional. After these procedures the interaction Lagrangian density becomes

$`_{int}(\phi ,\varphi )`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\lambda _0\varphi _a}{3N}}+{\displaystyle \frac{\eta _0\varphi ^2\varphi _a}{30N^2}}\right)\phi _a\phi ^2\left({\displaystyle \frac{8\eta _0\varphi _a\varphi _b\varphi _c}{6N^2}}\right)\phi _a\phi _b\phi _c{\displaystyle \frac{1}{4!N}}\left(\lambda _0+{\displaystyle \frac{\eta _0\varphi ^2}{10N}}\right)\phi ^4`$ (5)
$`\left({\displaystyle \frac{12\eta _0\varphi _a\varphi _b}{6!N^2}}\right)\phi _a\phi _b\phi ^2{\displaystyle \frac{1}{5!}}\left({\displaystyle \frac{\eta _0\varphi _a}{N^2}}\phi _a\phi ^4\right){\displaystyle \frac{\eta _0}{6!N^2}}\phi ^6.`$
The effective action $`\mathrm{\Gamma }(\varphi )`$ is found by solving for $`G_{ab}(x,y)`$ the equation

$$\frac{\delta \mathrm{\Gamma }(\varphi ,G)}{\delta G_{ab}(x,y)}=0,$$
(6)
and substituting the solution in the generalized effective action $`\mathrm{\Gamma }(\varphi ,G)`$.

The vertices in the above equation contain factors of $`1/N`$ or $`1/N^2`$, but a $`\phi `$ loop gives a factor of N provided the $`O(N)`$ isospin flows around it alone and not into another part of the graph. We usually call such loops bubbles. Then at leading order in $`1/N`$, the vacuum graphs are bubble trees with two or three bubbles at each vertex. The (2PI) graphs are shown in figure.(1).

It is straightforward to obtain

$$\mathrm{\Gamma }_2(\varphi ,G)=\frac{1}{4!N}d^Dx\left(\lambda _0+\frac{\eta _0\varphi ^2}{10N}\right)[G_{aa}(x,x)]^2\frac{\eta _0}{6!N^2}d^Dx[G_{aa}(x,x)]^3.$$
(7)
Therefore eq.(6) becomes

$`{\displaystyle \frac{\delta \mathrm{\Gamma }(\varphi ,G)}{\delta G_{ab}(x,y)}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(G^1)_{ab}(x,y)+{\displaystyle \frac{i}{2}}D^1(\varphi ){\displaystyle \frac{1}{12N}}\left(\lambda _0+{\displaystyle \frac{\eta _0\varphi ^2}{10N}}\right)[\delta _{ab}G_{cc}(x,x)]\delta ^D(xy)`$ (8)
$`{\displaystyle \frac{3\eta _0}{6!N}}\delta _{ab}[G_{cc}(x,x)]^2\delta ^D(xy)=0.`$
Rewriting this equation, we obtain the gap equation

$`(G^1)_{ab}(x,y)`$ $`=`$ $`D_{ab}^1(\varphi ;x,y)+{\displaystyle \frac{i}{6N}}\left(\lambda _0+{\displaystyle \frac{\eta _0\varphi ^2}{10N}}\right)[\delta _{ab}G_{cc}(x,x)]\delta ^D(xy)+`$ (9)
$`{\displaystyle \frac{i\eta _0}{5!N^2}}\delta _{ab}[G_{cc}(x,x)]^2\delta ^D(xy).`$
Hence

$$\frac{i}{2}TrD^1G=\frac{1}{12N}d^Dx\left(\lambda _0+\frac{\eta _0\varphi ^2}{10N}\right)[G_{aa}(x,x)]^2+\frac{3\eta _0}{6!N^2}d^Dx[G_{aa}(x,x)]^3+cte.$$
(10)
Using eqs.(9) and eq.(10) in eq.(7) we find the effective action

$`\mathrm{\Gamma }(\varphi )`$ $`=`$ $`I(\varphi )+{\displaystyle \frac{i}{2}}Tr[LnG^1]+{\displaystyle \frac{1}{4!N}}{\displaystyle d^Dx\left(\lambda _0+\frac{\eta _0\varphi ^2}{10N}\right)[G_{aa}(x,x)]^2}+`$ (11)
$`{\displaystyle \frac{2\eta _0}{6!N^2}}{\displaystyle d^Dx[G_{aa}(x,x)]^3},`$
where $`G_{ab}`$ is given implicitly by eq.(9). The trace in (11) are both the functional and the internal space. The last two terms on the r.h.s of eq.(11) are the leading contribution to the effective action in the $`1/N`$ expansion. As usual we may simplify the situation by separating $`G_{ab}`$ into transverse and longitudinal components, so

$$G_{ab}=(\delta _{ab}\frac{\varphi _a\varphi _b}{\varphi ^2})g+\frac{\varphi _a\varphi _b}{\varphi ^2}\stackrel{}{g},$$
(12)
in this form we can invert $`G_{ab}`$,

$$(G)_{ab}^1=(\delta _{ab}\frac{\varphi _a\varphi _b}{\varphi ^2})g^1+\frac{\varphi _a\varphi _b}{\varphi ^2}\stackrel{}{g}^1.$$
(13)
Now we can take the trace with respect to the indices of the internal space,

$$G_{aa}=Ng+O(1),(G)_{aa}^1=Ng^1+O(1).$$
(14)
From this equation at leading order in $`1/N`$, $`G_{ab}`$ is diagonal in $`a,b`$. Substituting eq.(14) into eq.(11) and eq.(9) and keeping only the leading order one finds that the daisy and superdaisy resummed effective potential for the $`\phi ^6`$ theory is given by:

$`\mathrm{\Gamma }(\varphi )`$ $`=`$ $`I(\varphi )+{\displaystyle \frac{iN}{2}}tr(\mathrm{ln}g^1)+{\displaystyle \frac{N}{4!}}{\displaystyle d^Dx\left(\lambda _0+\frac{\eta _0\varphi ^2}{10N}\right)g^2(x,x)}+`$ (15)
$`{\displaystyle \frac{2N\eta _0}{6!}}{\displaystyle d^Dxg^3(x,x)}+O(1),`$
where the trace is only in the functional space, and the gap equation becomes

$$g^1(x,y)=i\left[\mathrm{}+m_0^2+\frac{\lambda _0}{6}(\frac{\varphi ^2}{N}+g(x,x))+\frac{\eta _0}{5!}(\frac{\varphi ^2}{N}+g(x,x))^2\right]\delta ^D(xy)+O(\frac{1}{N}).$$
(16)
It is important to point out that this calculation was done by Townsend . We interested to generalize these results assuming that the system is in equilibrium with a thermal bath a temperature $`T=\beta ^1`$. Since we are studying the equilibrium situation it is convenience to use the Matsubara formalism. Consequently it is convenient to continue all momenta to Euclidean values $`(p_0=ip_4)`$ and take the following Ansatz for $`g(x,y)`$,

$$g(x,y)=\frac{d^Dp}{(2\pi )^D}\frac{\mathrm{exp}^{i(xy)p}}{p^2+M^2(\varphi )}.$$
(17)
Substituting eq.(17) in eq.(16) we get the expression for the gap equation:

$$M^2(\varphi )=m_0^2+\frac{\lambda _0}{6}\left(\frac{\varphi ^2}{N}+F(\varphi )\right)+\frac{\eta _0}{5!}\left(\frac{\varphi ^2}{N}+F(\varphi )\right)^2,$$
(18)
where $`F(\varphi )`$ is given by

$$F(\varphi )=\frac{d^Dp}{(2\pi )^D}\frac{1}{p^2+M^2(\varphi )},$$
(19)
and the effective potential in the D-dimensional Euclidean space can be expressed as

$$V(\varphi )=V_0(\varphi )+\frac{N}{2}\frac{d^Dp}{(2\pi )^D}\mathrm{ln}\left[p^2+M^2(\varphi )\right]\frac{N}{4!}(\lambda _0+\frac{\eta _0\varphi ^2}{10N})F(\varphi )^2\frac{2N\eta _0F(\varphi )^3}{6!},$$
(20)
where $`V_0(\varphi )`$ is the classical potential. In the next section using the Matsubara formalism we present the effective potential of the model at finite temperature.

## 3 The effective potential for $`\phi ^6`$ theory at finite temperature

Let us suppose that our system is in equilibrium with a thermal bath. To study the temperature effects in quantum field theory we will use the imaginary time Green function approach . In this formalism the Euclidean-time $`\tau `$ is restricted to the interval $`0\tau \beta `$, and the bosonic filed satisfies periodic boundary conditions in Euclidean-time. This is equivalence to replace the continuous four momenta $`k_4`$ by discrete $`\omega _n`$ and the integration by a summation ($`\beta =\frac{1}{T}`$):

$`k_4`$ $``$ $`\omega _n={\displaystyle \frac{2\pi n}{\beta }},n=0,\pm 1,\pm 2,\mathrm{}`$
$`{\displaystyle \frac{d^Dk}{(2\pi )^D}}`$ $``$ $`{\displaystyle \underset{n}{}}{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{d^{D1}k}{(2\pi )^{D1}}}.`$ (21)
It is important to stress that all the Feynman rule are the same as the temperature case, except, as we stressed that momentum space integrals over the zeroth component is replace by summ over discret summs. The effective potential at finite temperature can be write as:

$`V_\beta (\varphi )`$ $`=`$ $`V_0(\varphi )+{\displaystyle \frac{N}{2\beta }}{\displaystyle \underset{n}{\overset{\mathrm{}}{}}}{\displaystyle \frac{d^{D1}p}{(2\pi )^{D1}}\mathrm{ln}\left[\omega _n+p^2+M_\beta ^2(\varphi )\right]}`$ (22)
$`{\displaystyle \frac{N}{4!}}(\lambda _0+{\displaystyle \frac{\eta _0\varphi ^2}{10N}})F_\beta (\varphi )^2{\displaystyle \frac{2N\eta F_\beta (\varphi )^3}{6!}}.`$
where $`F_\beta (\varphi )`$ is a finite temperature generalization of $`F(\varphi )`$, where,

$$F_\beta (\varphi )=\frac{1}{\beta }\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\frac{d^{D1}p}{(2\pi )^{D1}}\frac{1}{\omega _n^2+p^2+M_\beta ^2(\varphi )}.$$
(23)
The gap equation for this theory at finite temperature is given by,

$$M_\beta ^2(\varphi )=m_0^2+\frac{\lambda _0}{6}\left(\frac{\varphi ^2}{N}+F_\beta (\varphi )\right)+\frac{\eta _0}{5!}\left(\frac{\varphi ^2}{N}+F_\beta (\varphi )\right)^2.$$
(24)
In order to regularize $`F_\beta (\varphi )`$ given by eq.(23), we use a mixing between dimensional regularization and analytic regularization. For this purpose we define the expression $`I_\beta (D,s,m)`$ as :

$$I_\beta (D,s,m)=\frac{1}{\beta }\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\frac{d^{D1}k}{(2\pi )^{D1}}\frac{1}{(\omega _n^2+k^2+m^2)^s}.$$
(25)
Using the analytic extension of the inhomogeneous Epstein zeta function it is possible to obtain the analytic extension of $`I_\beta (D,s,m)`$;

$$I_\beta (D,s,m)=\frac{m^{D2s}}{(2\pi ^{1/2})^D\mathrm{\Gamma }(s)}\left[\mathrm{\Gamma }(s\frac{D}{2})+4\underset{n=1}{\overset{\mathrm{}}{}}\left(\frac{2}{mn\beta }\right)^{D/2s}K_{D/2s}(mn\beta )\right]$$
(26)
where $`K_\mu (z)`$ is the modified Bessel function of the third kind. Fortunately for $`D=3`$ the analytic extension of the function $`I_\beta (D,s=1,m=M_\beta (\varphi ))=F_\beta (\varphi )`$ is finite and can be expressed in a closed form  (note that in $`D=3`$ we have no pole, at least in this approximation), and in particular as

$$F_\beta (\varphi )=I_\beta (3,1,M_\beta (\varphi ))=\frac{M_\beta (\varphi )}{4\pi }\left(1+\frac{2\mathrm{ln}(1e^{M_\beta (\varphi )\beta })}{M_\beta (\varphi )\beta }\right).$$
(27)
This result is not a peculiarity of this method of regularization, because this happens also in dimensional regularization at zero temperature, that is, in odd dimensions integrals which are divergent by naive power counting may be to regulated to finite value with no poles occurring, for example in $`D=3`$. In order to regularize the second term of eq.(22), we use the following method: We define,

$$LF_\beta (\varphi )=\frac{1}{\beta }\underset{n=1}{\overset{\mathrm{}}{}}\frac{d^{D1}p}{(2\pi )^{D1}}\mathrm{ln}\left[\omega _n+p^2+M_\beta ^2(\varphi )\right]$$
(28)
then,

$$\frac{LF_\beta (\varphi )}{M_\beta }=(2M_\beta )\frac{1}{\beta }\underset{n=1}{\overset{\mathrm{}}{}}\frac{d^{D1}p}{(2\pi )^{D1}}\frac{1}{\omega _n+p^2+M_\beta ^2(\varphi )}$$
(29)
and from the equation (23), we have that,

$$\frac{LF_\beta (\varphi )}{M_\beta }=(2M_\beta )F_\beta (\varphi ),$$
(30)
in this way the effective potential could be regularized. For $`D=3`$, $`F_\beta (\varphi )`$ is finite and is given by eq. (31) and integrating the eq.(30), we obtain:

$$LF_\beta (\varphi )_R=\frac{M_\beta (\varphi )^3}{6\pi }\frac{M_\beta (\varphi )Li_2(e^{M_\beta (\varphi )\beta })}{\pi \beta ^2}\frac{Li_3(e^{M_\beta (\varphi )\beta })}{\pi \beta ^3}.$$
(31)
The definition of general polylogarithm function $`Li_n(z)`$ can be found in ref. .

The daisy and super daisy resummed effective potential at finite temperature for $`D=3`$ is given by:

$$V_\beta (\varphi )=V_0(\varphi )+\frac{N}{2}LF_\beta (\varphi )_R\frac{N}{4!}(\lambda _0+\frac{\eta _0\varphi ^2}{10N})(F_\beta (\varphi )_R)^2\frac{2N\eta (F_\beta (\varphi )_R)^3}{6!}.$$
(32)
and the gap equation (see eq.(24)):

$`M_\beta ^2(\varphi )`$ $`=`$ $`m_0^2+{\displaystyle \frac{\lambda _0}{6}}\left({\displaystyle \frac{\varphi ^2}{N}}{\displaystyle \frac{M_\beta (\varphi )}{4\pi }}\left[1+{\displaystyle \frac{2\mathrm{ln}(1e^{M_\beta (\varphi )\beta })}{M_\beta (\varphi )\beta }}\right]\right)`$ (33)
$`+{\displaystyle \frac{\eta _0}{5!}}\left({\displaystyle \frac{\varphi ^2}{N}}{\displaystyle \frac{M_\beta (\varphi )}{4\pi }}\left[1+{\displaystyle \frac{2\mathrm{ln}(1e^{M_\beta (\varphi )\beta })}{M_\beta (\varphi )\beta }}\right]\right)^2.`$
For the case $`D=4`$, and for $`\eta =0`$, where the theory is just renormalizable, the effective potential can be obtain in the same way, that is:

$$V_\beta (\varphi )=V_0(\varphi )+\frac{N}{2}(LF_\beta (\varphi )_R)_{D=4}\frac{\lambda _0N}{4!}(F_\beta (\varphi )_R)_{D=4}^2,$$
(34)
where $`V_0`$ is the classical potential,

$$V_0(\varphi )=\frac{1}{2}m^2\phi ^2+\frac{\lambda }{4!N}\phi ^4,$$
(35)
and $`(F_\beta (\varphi )_R)_{D=4}`$ is given by:

$$(F_\beta (\varphi )_R)_{D=4}=\frac{\lambda M_\beta ^2(\varphi )}{2\pi ^2}_1^{\mathrm{}}\frac{(p^21)^{\frac{1}{2}}}{e^{M_\beta (\varphi )\beta p}1}𝑑p.$$
(36)
and in the limit of high temperature, we could write the above equation as:

$$(F_\beta (\varphi )_R)_{D=4}=\frac{1}{12\beta ^2}\frac{M_\beta (\varphi )}{4\pi \beta }\frac{M_\beta ^2(\varphi )}{8\pi ^2}\mathrm{ln}(M_\beta (\varphi )\beta );$$
(37)
and,

$$(LF_\beta (\varphi )_R)_{D=4}=\frac{M_\beta (\varphi )^2}{12\beta ^2}\frac{M_\beta (\varphi )^3}{6\pi \beta }\frac{M_\beta (\varphi )^4\mathrm{ln}(M_\beta (\varphi )\beta )}{16\pi ^2}+\frac{M_\beta (\varphi )^4}{64\pi ^2},$$
(38)
and the gap equation for $`D=4`$ is given by:

$$M_\beta ^2(\varphi )=\stackrel{~}{m}^2(\varphi )+\frac{\lambda }{6}\left(\frac{1}{12\beta ^2}\frac{M_\beta (\varphi )}{4\pi \beta }\frac{M_\beta ^2(\varphi )}{8\pi ^2}\mathrm{ln}(M_\beta (\varphi )\beta )\right),$$
(39)
where $`\stackrel{~}{m}^2(\varphi )=m^2+\frac{\lambda }{6}\frac{\varphi ^2}{N}`$, and $`m^2`$, $`\lambda `$ are the renormalized mass and coupling constant at zero temperature respectively. We note that, from the gap equation in eq.(39), we find that for the coupling constant $`\lambda 1`$, the condition $`M_\beta ^2(\varphi )/T^21`$ is consistent with $`\stackrel{~}{m}^2(\varphi )/T^2`$, which is exactly the required condition for the high temperature expansion .

## 4 Conclusions

In this paper we have performed an analysis of the daisy and super daisy effective potential for the theory $`\phi ^4`$ and $`\phi ^6`$ in $`D`$-dimensional Euclidean space at finite temperature. The form of effective potential have been found explicitly using resummation method in the leading order $`1/N`$ approximation (Hartree-Fock approximation). We have seen how dimensional regularization and analytic regularization can be used to compute the effective potential at finite temperature. In odd dimensional theory when power counting indicates that the diverges should occur, dimensional regularization and analytic does not give rise to a pole.

## 5 Acknowlegements

This paper was supported by Conselho Nacional de Desenvolvimento Cientifico e Tecnologico do Brazil (CNPq).