text:
# Relaxation and Noise in Chaotic Systems

## I Introduction

The purpose of the lectures that are summarized in this review is to describe the behavior of ensembles of chaotic systems. As in the case of statistical mechanics, the dynamics of ensembles turns out to be simpler then the one of individual systems. For chaotic systems  the long time asymptotic behavior is rigorously known to exhibit exponential decay of correlations only for a class of idealized systems (hyperbolic systems). Physical model systems do not belong to this class. In what follows nonrigorous methods, that enable the exploration of the long time behavior of chaotic systems, are presented and their application is demonstrated for several systems. The review is pedagogical and descriptive in nature, and is intended for an overview of the subject. The reader should consult the references for the details and the precise statements.

The dynamics of the systems considered in this review is of Hamiltonian nature. The dynamics of the continuous systems is determined by the Hamilton equations:

$`\dot{p}=`$ $`{\displaystyle \frac{}{q}}`$ (1)
$`\dot{q}=`$ $`{\displaystyle \frac{}{p}}`$ (2)
where $``$ is the Hamiltonian, while $`q`$ and $`p`$ are the position and momentum respectively. The phase space density $`\rho (𝐱)`$, where $`𝐱=(q,p)`$ are the phase space points, satisfies the Liouville theorem,

$$\frac{d\rho }{dt}=0$$
(3)
or

$$\frac{\rho }{t}=\{,\rho \},$$
(4)
where $`\{\mathrm{}\}`$ are the Poisson brackets. The main issue that will be discussed is the way $`\rho `$ spreads in phase space for chaotic systems. For classical chaotic systems finer and finer structures develop for longer and longer times. These structures are reflected in $`\rho `$. Coarse graining, with some fixed scale of resolution, results in the truncation of the evolution of these very fine structures. On the finite (but arbitrarily small) resolution scale the equilibrium uniform phase space density is approached. In this review the asymptotic (in time) relaxation to this density will be discussed and in the end the coarse graining scale will be taken to $`0`$. The results differ from the ones obtained without coarse graining, since the limits of infinite time and vanishing coarse graining do not commute. The approach, where a finite resolution scale is used and then the limit where this scale tends to zero is taken, is relevant for experimental realizations where the idealized classical description on the finest scales is destroyed resulting of the coupling to the environment.

Maps are transformations of phase space in discrete time, denoted by $`n`$. A map $`𝐅`$ is a transformation

$$𝐱_{n+1}=𝐅(𝐱_n).$$
(5)
We study maps since they are easier to handle analytically and numerically and they reproduce the most interesting results found for systems evolving continuously in time.

Maps can be derived from Hamiltonians of the form:

$$=\frac{p^2}{2}+V(q)\underset{n}{}\delta (tn).$$
(6)
For such maps the phase space area is preserved. Maps can also be defined with no reference to a Hamiltonian. The dynamics of phase space densities will be explored here for several area preserving maps: a. Kicked Rotor (Standard Map), b. Kicked Top, c. Arnold Cat Map and d. Baker Map.

a. Kicked Rotor (Standard Map)
A planar rotor that is periodically kicked is modeled by the Hamiltonian

$$=\frac{J^2}{2}+K\mathrm{cos}\theta \underset{n}{}\delta (tn),$$
(7)
where $`\theta `$ is the coordinate and $`J`$ is the conjugate momentum. The Hamilton equations (1) are:

$`\dot{\theta }`$ $`=`$ $`J`$ (8)
$`\dot{J}`$ $`=`$ $`K\mathrm{sin}\theta {\displaystyle \underset{n}{}}\delta (tn).`$ (9)
Integration with respect to time results in the Standard Map

$`\theta _{n+1}`$ $`=`$ $`\theta _n+J_n`$ (10)
$`J_{n+1}`$ $`=`$ $`J_n+K\mathrm{sin}\theta _{n+1},`$ (11)
where $`\theta _n`$ and $`J_n`$ are the angle and angular momentum just after the $`n`$-th kick. The equations (10) define a map of the form (5)

$$(\theta _{n+1},J_{n+1})=𝐅(\theta _n,J_n).$$
(12)
It is easily checked that the map is area preserving. It becomes more chaotic as the stochasticity parameter $`K`$ increases. The phase space is plotted in Fig. 1.

A variant of this map provides a good description of driven laser cooled atoms  and of beams deflected by dielectrics with a modulated index of refraction .

b. Kicked Top
A large spin can be described by a classical vector

$$𝐉=j(\mathrm{sin}\theta \mathrm{cos}\phi ,\mathrm{sin}\theta \mathrm{sin}\phi ,\mathrm{cos}\theta )$$
(13)
where $`\theta `$ and $`\phi `$ are the polar angles. The kicked top map is defined by the transformation

$$𝐅=R_z(\tau \mathrm{cos}\theta )R_z(\beta _z)R_y(\beta _y),$$
(14)
applied to $`𝐉`$ where $`R_i(\beta )`$ is the rotation around the axis $`i`$ by the angle $`\beta `$. The nonlinearity results of the dependence of a rotation on the angle $`\theta `$. The chaoticity of the map increases with $`\tau `$ that is the stochasticity parameter. The canonical phase space variables for this map are $`q=\phi `$ and $`p=\mathrm{cos}\theta `$. The phase space is plotted in Fig. 2.

c. Arnold Cat Map
It is defined in the $`[0,1]\times [0,1]`$ square of the $`(x,y)`$ phase plane by

$`x_{n+1}=`$ $`x_n+y_n`$ $`\text{mod}1`$ (15)
$`y_{n+1}=`$ $`x_n+2y_n`$ $`\text{mod}1.`$ (16)
The evolution is demonstrated in .

d. Baker map
It is defined in the $`[0,1]\times [0,1]`$ square of the $`(x,y)`$ phase plane by

$$(x_{n+1},y_{n+1})=𝐅(x_n,y_n)=\{\begin{array}{cc}(2x_n,y_n/2)\hfill & \text{ for }0x<\frac{1}{2}\hfill \\ (2x_n1,(y_n+1)/2)\hfill & \text{ for }\frac{1}{2}x<1\hfill \end{array}$$
(17)
This transformation is demonstrated in Fig. 3.

The systems (a) and (b) model physical problems. Their phase space is mixed, namely the dynamics in some parts is chaotic and in some parts it is regular. Systems (c) and (d) are very idealized. For these systems the motion is chaotic in the entire phase space.

A system is chaotic if points, that are initially close, spread exponentially in phase space. For chaotic systems the separation between points grows with the number of steps as

$$ϵ(n)=ϵ(0)\mathrm{\Lambda }^n,$$
(18)
in the limit $`n\mathrm{}`$ and $`ϵ(0)0`$. The Lyapunov exponent is

$$\lambda =\mathrm{ln}\mathrm{\Lambda }$$
(19)
and $`\mathrm{\Lambda }`$ is sometimes called the Lyapunov number. A system that evolves in a bounded region in phase space is called chaotic if its Lyapunov exponent $`\lambda `$ according to the definition (18,19) is positive. (More generally the Lyapunov numbers are the eigenvalues of the tangent map. The Lyapunov number defined by (18) is equal to the magnitude of the largest eigenvalue.) One can study the local expansion and contraction in the various directions. If for a map, for any point of phase space, in each direction, either expansion or contraction takes place, or in other words if none of the real parts of the local Lyapunov exponents vanishes, the system is called hyperbolic.

Exponential spread of trajectories in phase space does not imply exponential decay of correlations in time. Such correlations may result of sticking to small structures in phase space, as is the case for the stadium billiard and for the Sinai billiard as well as for chaotic trajectories of mixed systems, such as the kicked rotor and the kicked top. For hyperbolic systems, such as the baker map, typically correlations in time decay exponentially.

Chaotic motion is reminiscent of a random walk. Therefore it is instructive to review the dynamics of the probability density of random walkers. Assume a random walk on a one dimensional lattice, with equal probability of $`1/2`$ to make a step to the left or to the right. The path of each walker is complicated but the evolution of a probability density of random walkers is simple. It is simple in particular in the continuum approximation (the limit of vanishing lattice spacing) where it is described by the diffusion equation

$$\frac{}{t}\rho (x,t)=D\frac{^2}{x^2}\rho (x,t),$$
(20)
$`D`$ is the diffusion coefficient and $`\rho (x,t)`$ is the density of random walkers at the point $`x`$ at time $`t`$. Let us assume that the random walk is confined to an interval of length $`s`$, for example $`s/2<x<s/2`$. The current at the ends of the interval vanishes, resulting in the boundary conditions

$$\frac{\rho }{x}|{}_{x=\pm s/2}{}^{}=0$$
(21)
on (20). This problem is similar to the standard quantum mechanical problem of a particle in an infinite square well, but with the boundary conditions (21) and with no $`i\mathrm{}`$ on the left hand side of (20). Therefore the density of states can be expanded as:

$$\rho (x,t)=\underset{k=0}{\overset{\mathrm{}}{}}a_ke^{\gamma _kt}u_k(x),$$
(22)
where

$$u_k=\{\begin{array}{cc}\sqrt{\frac{2}{s}}\left[\mathrm{cos}\frac{\pi x}{s}k+\left(\sqrt{\frac{1}{2}}1\right)\delta _{k0}\right]\hfill & k\text{even}\hfill \\ \sqrt{\frac{2}{s}}\mathrm{sin}\frac{\pi x}{s}k\hfill & k\text{odd}\hfill \end{array}$$
(23)
form an orthonormal basis and

$$\gamma _k=\left(\frac{\pi k}{s}\right)^2D.$$
(24)
The units of $`a_k`$ are of $`1/\sqrt{\text{[length]}}`$. The expansion coefficients $`a_k`$ are determined by the initial density $`\rho (x,t=0)`$. The equilibrium density is

$$\underset{t\mathrm{}}{lim}\rho (x,t)=a_0\sqrt{\frac{1}{s}},$$
(25)
and the asymptotic approach to this density is

$$\rho (x,t)\frac{a_0}{\sqrt{s}}a_1e^{\gamma _1t}u_1(x)$$
(26)
with the relaxation rate $`\gamma _1=\left(\frac{\pi }{s}\right)^2D`$. The density-density correlation function is

$$C(t)=_{s/2}^{s/2}𝑑x\rho (x,0)\rho (x,t)a_0^2=\underset{k=1}{\overset{\mathrm{}}{}}|a_k|^2e^{\gamma _kt}|a_1|^2e^{\gamma _1t}.$$
(27)
Hence also the correlations decay with the rate $`\gamma _1`$.

The example of the probability density of random walkers demonstrates that although the trajectories of specific random walkers are complicated their probability density $`\rho `$ follows simple dynamics. We turn now to study the phase space probability density $`\rho (𝐱)`$, that for Hamiltonian dynamics satisfies (4). The introduction of such probability densities, that are sufficiently smooth, is actually a coarse graining over some fine scale. It is an averaging over phase space of the same type that was performed for random walkers. For any given area preserving map $`𝐅`$ it is instructive to introduce the one step evolution operator $`\widehat{U}`$,

$$\rho _{n+1}(𝐱)=\widehat{U}\rho _n(𝐱),$$
(28)
where $`\rho _n(𝐱)`$ is the phase space density after $`n`$ steps of the map. If the operator is defined on a space of sufficiently smooth functions it is called the Frobenius-Perron operator. If $`𝐅`$ is area preserving and invertible, namely if $`𝐅^1`$ exists, then $`\widehat{U}`$ is unitary. All eigenvalues of $`\widehat{U}`$ are on the unit circle in the complex plane, they take the form $`e^{i2\pi \alpha }`$, with real $`\alpha `$. The corresponding “eigenfunctions” are not square integrable. An example of such a function will be presented now . Assume the map $`𝐅`$ has a periodic orbit of period $`n`$, that is:

$`𝐱_{j+1}`$ $`=`$ $`𝐅(𝐱_j),j=1,2,\mathrm{}n`$ (29)
$`𝐱_1`$ $`=`$ $`𝐅(𝐱_n).`$ (30)
The function

$$\psi _\alpha (𝐱)=\underset{j=1}{\overset{n}{}}\delta (𝐱𝐱_j)e^{i2\pi \alpha j}$$
(31)
with $`\alpha =l/n`$, where $`l`$ is an integer, is an “eigenfunction” of $`\widehat{U}`$ with the eigenvalue $`e^{i2\pi \alpha }`$. Here $`\psi _\alpha `$ is called “eigenfunction” if

$$\widehat{U}\psi _\alpha =e^{i2\pi \alpha }\psi _\alpha $$
(32)
and no requirement is made that it belongs to the space where $`\widehat{U}`$ is defined. For example if $`\widehat{U}`$ is defined on the space of square integrable functions, $`\psi _\alpha `$ of (31) does not belong to this space and is not a function in the usual sense (it is a distribution).

The operator $`\widehat{U}`$ can be represented by an infinite dimensional matrix. In physics applications it is natural to use finite dimensional approximations. Let us assume that a basis that is ordered by increased resolution is used (the basis states may be for example sines and cosines or orthogonal polynomials). For the truncation of $`\widehat{U}`$, namely a finite dimensional approximation, the matrix is not unitary, and its eigenvalues are inside the unit circle in the complex plane. They may vary with the dimension of the matrix. The natural question is: Do these eigenvalues approach the unit circle in the limit of an infinite dimensional matrix? It will be demonstrated in what follows that this is not the case for chaotic systems and values inside the unit circle are approached in this limit, in spite of the fact that the infinite dimensional matrix representing $`\widehat{U}`$ is unitary. A heuristic justification for this behavior was proposed by F. Haake . Chaotic systems, in contrast with regular ones, exhibit phase space structures on all scales. These structures are revealed during the evolution. The operator $`\widehat{U}`$ couples the fine scales via its matrix elements that couple to states with high resolution. As a result of the truncation, probability that was originally transferred to the fine scales is lost, resulting in nonunitarity of $`\widehat{U}`$. Convergence of the eigenvalues to values inside the unit circle in the complex plane, in the limit of infinite dimension of the matrix, results of the asymptotic self similarity of the chaotic dynamics. For regular motion, on the other hand, as the dimension of the matrix is increased the eigenvalues approach the unit circle  (see Sec. IV).

It will be shown that for chaotic systems the finite dimensional approximations of $`\widehat{U}`$, in the limit of large dimension, describe the decay of correlations. The phase space density-density correlation function, in analogy with (27) is

$`C(n)`$ $`=`$ $`{\displaystyle 𝑑𝐱\rho (𝐱,0)\rho (𝐱,n)}\rho _{\mathrm{}}^2\mathrm{\Omega }`$ (33)
$`=`$ $`{\displaystyle 𝑑𝐱\rho (𝐱,0)\widehat{U}^n\rho (𝐱,0)}\rho _{\mathrm{}}^2\mathrm{\Omega },`$ (34)
where $`\mathrm{\Omega }`$ is the volume of the chaotic component in phase space. The equilibrium density is $`\rho _{\mathrm{}}=lim_n\mathrm{}\rho (𝐱,n)`$ that is independent of position in phase space. A more general correlation function is:

$$C^{(A,B)}(n)=𝑑𝐱A(𝐱)\widehat{U}^nB(𝐱),$$
(35)
where we assumed

$$\underset{n\mathrm{}}{lim}\widehat{U}^nA(𝐱)=\underset{n\mathrm{}}{lim}\widehat{U}^nB(𝐱)=0.$$
(36)
This can always be obtained by the subtraction of the asymptotic value. For simplicity it will be assumed that both $`A`$ and $`B`$ are real.

It is useful to study the Laplace transforms of the correlation functions. For this purpose we introduce the resolvent

$$\widehat{R}(z)=\frac{1}{z}\underset{j=0}{\overset{\mathrm{}}{}}\widehat{U}^jz^j=\frac{1}{z\widehat{U}}$$
(37)
where $`z`$ is a complex number. Since $`\widehat{U}`$ is unitary the sum is convergent for $`|z|>1`$. This is analogous to the usual definition used in quantum mechanics:

$$\widehat{R}(E)=\frac{1}{i\mathrm{}}_0^{\mathrm{}}𝑑t\widehat{U}\mathrm{exp}\left(\frac{i}{\mathrm{}}Et\frac{ϵ}{\mathrm{}}t\right),$$
(38)
where

$$\widehat{U}=\mathrm{exp}\left(\frac{i}{\mathrm{}}\widehat{}t\right)$$
(39)
is the evolution operator, leading to

$$\widehat{R}(E)=\frac{1}{E\widehat{}+iϵ}$$
(40)
for $`ϵ>0`$. The convergence of the integral (38) requires $`ϵ>0`$. In analogy the convergence of the sum (37) requires $`|z|>1`$ in (37). The discrete Laplace (one sided Fourier) transform of the correlation function (35) is

$$\stackrel{~}{C}^{(A,B)}(z)=\underset{n=0}{\overset{\mathrm{}}{}}C^{(A,B)}(n)z^n=𝑑𝐱A(𝐱)\left[z\widehat{R}(z)\right]B(𝐱),$$
(41)
as one finds from (37).

Let $`\psi _i`$ be an “eigenfunction” of $`\widehat{U}`$ with the eigenvalue $`z_i`$,

$$\widehat{U}\psi _i=z_i\psi _i.$$
(42)
By “eigenfunction” we mean here that it satisfies (42), but it may not be a function in the usual sense (for example it may be a distribution). Then $`z_i`$ is a pole of the matrix elements of the resolvent as one can see from (37). The correlation function, involving an “eigenfunction”, takes the form

$$\stackrel{~}{C}^{(A,\psi _i)}(z)=A|\psi _i\frac{z}{zz_i},$$
(43)
where the Dirac notation

$$A|B=𝑑𝐱A^{}(𝐱)B(𝐱)$$
(44)
is used ($`A^{}`$ is the complex conjugate of $`A`$). Here $`A`$ and $`B`$ were assumed to be real. Let us expand $`B(𝐱)`$ in terms of the “eigenfunctions” of $`\widehat{U}`$ as

$$B(𝐱)=\underset{i}{}b_i\psi _i(𝐱).$$
(45)
Then

$$\stackrel{~}{C}^{(A,B)}(z)=\frac{r_i}{zz_i}$$
(46)
where the poles $`z_i`$ depend only on $`\widehat{U}`$ while the residues

$$r_i=b_iA|\psi _iz,$$
(47)
depend also on $`A`$ and $`B`$.

Because of the unusual nature of the “eigenfunctions” $`\psi _i`$ the manipulations leading to (46) are only heuristic. The Ruelle-Pollicott theorem justifies the expression (46) for hyperbolic systems if $`A`$ and $`B`$ are typical smooth functions (for a precise statement of the theorem that is transparent for physicists see ). Moreover the theorem assures that $`|z_i|<1`$, except for the eigenvalue $`z_0=1`$ corresponding to the equilibrium density $`\psi _0`$, that is independent of $`𝐱`$. The $`z_i`$ are called the Ruelle-Pollicott resonances.

The existence of poles with $`|z_i|<1`$ implies the decay of correlations with the rate $`\mathrm{ln}|z_i|`$. This results of

$$C^{(A,B)}(n)=\frac{1}{2\pi i}_{|z|=1}\frac{dz}{z}\stackrel{~}{C}^{(A,B)}(z)z^n$$
(48)
and the application of the residue theorem to (46).

For any approximation of dimension $`N`$ of the Frobenius-Perron operator $`\widehat{U}`$, that will be denoted by $`\widehat{U}^{(N)}`$ in what follows, the eigenfunctions $`\psi _i^{(N)}`$, corresponding to the eigenvalues $`z_i^{(N)}`$, are well defined. For hyperbolic systems, where the Ruelle-Pollicott theorem applies, it is reasonable to assume that in the limit $`N\mathrm{}`$, the eigenvalues $`z_i^{(N)}`$ approach the Ruelle-Pollicott resonances, for a typical choice of the basis. Inspired by the Ruelle-Pollicott theorem the following heuristic scheme for the calculation of the decay rates of correlations is proposed:

1. Introduce an orthogonal basis where the basis states are ordered by resolution. These states, for example, may be orthogonal polynomials or trigonometric functions ($`\mathrm{cos}2\pi kx,\mathrm{sin}2\pi kx`$) or exponentials $`\mathrm{exp}i2\pi kx`$.
2. Calculate the matrix elements of $`\widehat{U}`$ in this basis.
3. Introduce a truncation of dimension $`N`$ of this matrix, that will be denoted as $`\widehat{U}^{(N)}`$.
4. Calculate the eigenvalues $`z_i^{(N)}`$ and the eigenfunctions $`\psi _i^{(N)}`$ of $`\widehat{U}^{(N)}`$. The eigenfunctions $`\psi _i^{(N)}`$ are finite linear combinations of the basis states, and therefore are smooth.
5. Take the limit $`N\mathrm{}`$. In this limit $`z_i^{(N)}z_i`$ and $`\psi _i^{(N)}\psi _i`$.

From the Ruelle-Pollicott theorem it is expected, for hyperbolic systems, that the Ruelle-Pollicott resonances are obtained in this way, if a typical basis is used. The limiting functions $`\psi _i`$ are singular . This scheme was applied to various systems, including also to mixed systems . Limiting eigenvalues $`z_i`$ satisfying $`|z_i|<1`$ were obtained in this way, implying exponential relaxation of correlations with the rate $`\mathrm{ln}|z_i|`$ that takes place for a very long time (for some nonhyperbolic systems it is known that eventually the correlations decay as a power-law, but the time required to obtain this power-law behavior may be too long for any physical relevance). Although the Hamiltonian dynamics of trajectories is reversible the phase space densities exhibit relaxation that is irreversible, and is similar to the behavior of probability densities of random walkers that follow irreversible dynamics.

In Sec. II the heuristic scheme for the calculation of the $`z_i`$ will be demonstrated for the baker map that is hyperbolic and consequently the Ruelle-Pollicott theorem applies. Exploration of a hyperbolic system that is a modification of the cat map is briefly mentioned in Sec. III. In Sec. IV the scheme will be applied to the kicked top and in Sec. V it will be applied to the kicked rotor, that are mixed systems, where the Ruelle-Pollicott theorem is not valid. Also the limitations on the validity of this scheme for such systems will be demonstrated. The main conclusions are summarized in Sec. VI. The review will follow closely references , that are marked by \** in the list of references. A detailed discussion on the relaxation of correlations in chaotic systems and on related topics can be found in . The scheme for the use of truncation in the calculation of the Ruelle-Pollicott resonances is discussed in a wider context in .

## II The Frobenius-Perron Operator for the Baker Map, a Demonstration

The heuristic prescription for the analysis of the Frobenius-Perron operator that was outlined in the Introduction will be demonstrated for the baker map (17), where all the results are exactly known. The analysis will make use of exact results obtained in .

The basis states that will be used are

$$xy|kl=\stackrel{~}{P}_k(x)\stackrel{~}{P}_l(y)$$
(49)
where

$$\stackrel{~}{P}_k(x)=\sqrt{2k+1}P_k(12x)=\frac{\sqrt{2k+1}}{k!}\frac{d^k}{dx^k}x^k(1x)^k$$
(50)
are the modified Legendre polynomials while $`P_k`$ are the Legendre polynomials. A basis of orthogonal polynomials is natural, since the order of the polynomial is not affected by the map (17). The orthonormality of the basis

$$kl|k^{}l^{}=\delta _{kk^{}}\delta _{ll^{}}$$
(51)
follows from the orthonormality of the modified Legendre polynomials

$$_0^1𝑑x\stackrel{~}{P}_k(x)\stackrel{~}{P}_k^{}(x)=\delta _{kk^{}}.$$
(52)
The basis is naturally ordered by increased resolution, since this is a property of orthogonal polynomials. The action of $`\widehat{U}`$ of (28) on any phase space density is

$$\widehat{U}\rho (𝐱)=\rho (𝐅^1(𝐱)),$$
(53)
since the points that are at $`𝐱`$ at time step $`n+1`$ were at $`𝐅^1(𝐱)`$ the time step $`n`$. For the baker map (17)

$$\rho (𝐅^1(x,y))=\{\begin{array}{cc}(\rho (x/2,2y)\hfill & \text{ for }0y<\frac{1}{2}\hfill \\ (\rho ((x+1)/2,2y1)\hfill & \text{ for }\frac{1}{2}y<1.\hfill \end{array}$$
(54)
From (53) and (54) one finds that the matrix elements of $`\widehat{U}`$ are

$$kl|\widehat{U}|k^{}l^{}=\frac{1}{2}\left[1+(1)^{k+k^{}+l+l^{}}\right]I_{kk^{}}I_{l^{}l}$$
(55)
where

$$I_{kk^{}}=_0^1𝑑x\stackrel{~}{P}_k(x)\stackrel{~}{P}_k^{}(x/2).$$
(56)
Using (50) and integrating (56) by parts $`k`$ times one finds that if $`k>k^{}`$,

$$I_{kk^{}}=0.$$
(57)
This results in the nonrecurrence property of $`\widehat{U}`$

$$kl|\widehat{U}|k^{}l^{}=0\text{for }k>k^{}\text{ or }l^{}>l.$$
(58)
During the evolution, probability is transformed from states with $`k^{}`$ to states with $`k`$ only if $`kk^{}`$, therefore after application of $`\widehat{U}`$ the density becomes more uniform because the weight of lower order Legendre polynomials is increased. This is expected since $`x`$ is the unstable direction, where stretching takes place, making the density more uniform. In the $`y`$ direction, on the other hand, the density is transformed from $`l^{}`$ to $`l`$ only if $`ll^{}`$, therefore during the evolution, the weight of the high order Legendre polynomials increases. This results of the fact that $`y`$ is the stable direction, where contraction takes place, resulting in complexity that increases with time.

For $`kk^{}`$ one finds

$$I_{kk^{}}=\frac{[(2k+1)(2k^{}+1)]^{1/2}}{2^k}\underset{l=0}{\overset{k^{}k}{}}\left(\frac{1}{2}\right)^l\frac{(k^{}+k+l)!}{(k^{}kl)!(2k+l+1)!l!}.$$
(59)
In particular $`I_{kk}=2^k`$ and the diagonal matrix elements are

$$kl|\widehat{U}|kl=\frac{1}{2^{k+l}}.$$
(60)
The matrix elements of $`\widehat{U}^n`$ can be calculated with the help of the nonrecurrence property (58). Introducing the resolution of the identity one finds

$$kl|\widehat{U}^n|k^{}l^{}=\underset{[k_i;l_i]}{}kl|\widehat{U}|k_1l_1k_1l_1|\widehat{U}|k_2l_2\mathrm{}..\mathrm{}k_{n1}l_{n1}|\widehat{U}|k^{}l^{},$$
(61)
where $`[k_i;l_i][k_1,k_2\mathrm{}..k_{n1};l_1,l_2,\mathrm{}.l_{n1}]`$ and the nonrecurrence property implies $`kk_1k_2\mathrm{}....k^{}`$ and $`ll_1l_2\mathrm{}.\mathrm{}l^{}`$. In particular the diagonal matrix elements satisfy

$$kl|\widehat{U}^n|kl=kl|\widehat{U}|kl^n=\left(\frac{1}{2^{k+l}}\right)^n.$$
(62)
In the limit $`n\mathrm{}`$ the off-diagonal matrix elements are dominated by powers of the diagonal matrix elements, since the nonrecurrence property limits the number of non-diagonal matrix elements in (61).

The diagonal matrix elements of the resolvent can be easily calculated with the help of (37). For $`|z|>1`$ the sum in (37) is convergent and using (62) one finds

$$kl|\widehat{R}(z)|kl=\underset{n=0}{\overset{\mathrm{}}{}}z^{(n+1)}kl|\widehat{U}|kl^n=\frac{1}{z2^{(k+l)}}.$$
(63)
These matrix elements are singular for $`z=2^{(k+l)}`$, and except for $`k=l=0`$ all these singular points are inside the unit circle in the complex plane. This is also the case for the off-diagonal matrix elements. Therefore the Ruelle-Pollicott resonances are

$$z_m=2^m$$
(64)
and their degeneracy is $`m+1`$. They are related to the decay of correlations by (46). They were obtained by analytic continuation from $`|z|>1`$ where the resolvent is defined. One can show  that the matrix elements exhibit a cut at $`|z|=1`$ and the physically relevant poles are on the Riemann sheet continued from $`|z|>1`$.

What happens if the Frobenius-Perron operator is restricted to $`k,k^{}k_{max}`$ and $`l,l^{}l_{max}`$ and $`\widehat{U}`$ is approximated by $`\widehat{U}^{(N)}`$, an $`N=k_{max}l_{max}`$ dimensional matrix, resulting of the truncation of $`\widehat{U}`$? The nonrecurrence property holds for $`\widehat{U}^{(N)}`$ and consequently (62) holds. Therefore for the $`N`$ dimensional matrix $`\widehat{U}^{(N)}`$ the diagonal matrix elements, are zeros of the characteristic polynomial, that is of order $`N`$. To see this, note that by the Hamilton-Cayley theorem, $`\widehat{U}^{(N)}`$ satisfies its characteristic polynomial, namely $`_{j=0}^Np_j\widehat{U}^{(N)j}=0`$, and calculate the diagonal matrix elements of this expression with the help of (62). Consequently the diagonal matrix elements are eigenvalues of $`\widehat{U}^{(N)}`$, taking the values $`2^m`$. The multiplicity of the eigenvalue $`2^m`$ is $`m+1`$ if $`mk_{max},l_{max}`$. By a similarity transformation the matrix can be transformed to the canonical Jordan form. Another way to see that the eigenvalues of $`\widehat{U}^{(N)}`$ coincide with the diagonal matrix elements is by using the indices $`l_{max}l`$ and $`l_{max}l^{}`$ instead of $`l`$ and $`l^{}`$. With these indices, $`\widehat{U}^{(N)}`$ is upper triangular, therefore the eigenvalues coincide with the diagonal matrix elements. The eigenvalues are independent of $`k_{max}`$ and $`l_{max}`$, therefore these are of the form $`2^m`$ also in the limit $`k_{max}\mathrm{}`$ and $`l_{max}\mathrm{}`$, in spite of the fact that $`\widehat{U}`$ is unitary. For $`\widehat{U}^{(N)}`$ the eigenfunctions are finite combinations of Legendre polynomials and are therefore smooth. In the limit $`k_{max}\mathrm{}`$ and $`l_{max}\mathrm{}`$ the right eigenfunctions are independent of $`x`$ and are polynomials of $`y`$, while the left eigenfunctions tend to distributions .

The main results that were obtained for the baker map are:

1. A basis that is ordered by increased resolution was introduced. A truncation of dimension $`N`$ in this basis was implemented.
2. The Ruelle-Pollicott resonances were calculated from $`\widehat{U}^{(N)}`$, the $`N`$ dimensional truncation of $`\widehat{U}`$. As $`N`$ increases more eigenvalues are revealed. The values of the eigenvalues do not depend on the truncation dimension $`N`$, and therefore also in the limit of infinite $`N`$ they remain at values that were obtained for finite $`N`$.
3. The poles of the matrix elements of the resolvent (37) and (63) were obtained by analytic continuation from $`|z|>1`$.
4. The eigenfunctions of $`\widehat{U}`$ are functions of $`y`$ only ($`y`$ is the stable direction).

In this section the procedure for the calculation of Ruelle-Pollicott resonances from $`\widehat{U}^{(N)}`$, a truncated matrix approximating $`\widehat{U}`$, was demonstrated for a system where the results are known exactly. In the following sections this method will be applied to systems where there is no exact theory.

The basis used here was of Legendre polynomials. If a basis of sines and cosines is used it is impossible to find the Ruelle-Pollicott resonances in this way, since each of the basis states collapses to $`0`$ after a finite number of applications of $`\widehat{U}`$ . On the basis of the Ruelle-Pollicott theorem, we believe that for a typical basis one does not encounter problems of this nature.

## III The Frobenius-Perron Operator of a Modified Cat Map

For the Arnold cat map, defined by (15), correlations in time decay faster than exponentially. A modification, where the function $`f(x_{n+1})`$, that is defined by

$$f(x)=\frac{K_0}{2\pi }\left[\mathrm{cos}(2\pi x)\mathrm{cos}(4\pi x)\right],$$
(65)
is added to the equation for $`y_{n+1}`$ in (15), is a hyperbolic system if $`K_0`$ is sufficiently small . The correlations for this system decay exponentially in time. It was studied with the help of the truncated Frobenius-Perron operator, according to the scheme that was outlined in the introduction. The leading Ruelle-Pollicott resonances were found with the help of a variational approach.

## IV The Frobenius-Perron Operator for the Kicked Top

We turn now to apply the method where $`\widehat{U}`$ is approximated by a finite dimensional matrix to mixed systems (where in some parts of phase space the motion is chaotic and in other parts it is regular). The Ruelle-Pollicott theorem does not apply to such systems. It will be demonstrated, however, that the Ruelle-Pollicott resonances are meaningful and describe the decay of correlations in the chaotic component of mixed systems for a time that may be very long. Asymptotically power-law decay takes place, due to sticking to regular islands in phase space.

In this section the kicked top , that was defined by (13) and (14), will be analyzed. It is a summary of results that were presented in  (see also  and ). The phase space is shown in Fig. 2.

A natural basis for this problem is of spherical harmonics

$$Y_{lm}(\theta ,\phi )=(1)^m\sqrt{\frac{2l+1}{4\pi }\frac{(lm)!}{(l+m)!}}P_l^m(\mathrm{cos}\theta )e^{im\phi }.$$
(66)
Alternatively one can choose a basis of real functions with $`e^{im\phi }`$ replaced by $`\mathrm{sin}m\phi `$ and $`\mathrm{cos}m\phi `$, resulting in a real matrix for $`\widehat{U}`$. It is naturally ordered by increased resolution. A truncation at $`l=l_{max}`$ is introduced. The dimension of the truncated matrix $`\widehat{U}^{(N)}`$, approximating $`\widehat{U}`$ is $`N=(l_{max}+1)^2`$. As $`l_{max}`$ is increased the eigenvalues of $`\widehat{U}^{(N)}`$ converge to values inside the unit circle in the complex plane, as can be seen from Table I.

The corresponding eigenfunctions are presented in Fig. 4.

Dark-shaded regions in phase space indicate large amplitudes of the eigenfunctions. In Figs. 5 (a) and (b)

it is demonstrated that as $`l_{max}`$ increases finer details of the eigenfunctions are revealed. This is a result of increased resolution. Comparison with Fig. 6 (b)

demonstrates that the eigenfunctions are large and tend to be uniform on the unstable manifold. In Fig. 7

the density-density correlation function (proportional to (33)) is plotted. It is demonstrated that also here the Ruelle-Pollicott resonances describe the decay of correlations for a very long time.

The behavior of the eigenfunctions localized in regular regions and the corresponding eigenvalues is very different. Eigenfunctions that are localized in islands of regular motion presented in Fig. 2 (a) are depicted in Fig. 8.

The eigenfunctions follow elliptic orbits inside the regular island and therefore do not involve the resolution of very fine structures. The corresponding eigenvalues approach the unit circle in the complex plane as $`l_{max}`$ increases.

It is instructive to examine the evolution operator $`\widehat{U}^1`$ corresponding to the inverse map $`𝐅^1`$. Unitarity implies

$$\widehat{U}^1=\widehat{U}^{}.$$
(67)
A basis can be chosen so that $`\widehat{U}`$ is real and $`\widehat{U}^{}=\widehat{U}^T`$, where $`T`$ denotes transpose. The eigenvalues of $`\widehat{U}^T`$ are equal to the ones of $`\widehat{U}`$. Introducing the truncation one finds

$$\widehat{U}^{T(N)}=\widehat{U}^{(N)T},$$
(68)
where $`\widehat{U}^{T(N)}`$ is the truncated $`\widehat{U}^T`$ and $`\widehat{U}^{(N)T}`$ is the transpose of $`\widehat{U}^{(N)}`$, the truncated $`\widehat{U}`$. The eigenfunctions of $`\widehat{U}^{T(N)}`$ are localized on the unstable manifold of $`𝐅^1`$, that is the stable manifold of $`𝐅`$, as is clearly demonstrated comparing Fig. 5 (c) with Fig. 6 (c). By eigenfunctions we meant so far the right eigenfunctions. The left eigenfunctions of $`\widehat{U}`$ (and $`\widehat{U}^{(N)}`$) are the right eigenfunctions of $`\widehat{U}^T`$ (and $`\widehat{U}^{T(N)}`$).

The main results that were obtained for the kicked top are:

1. As $`N`$ or $`l_{max}`$ increases the eigenvalues of $`\widehat{U}^{(N)}`$, corresponding to eigenfunctions localized in the chaotic component of phase space, converge to values inside the unit circle. Unlike the case of the baker map, the eigenvalues $`z_i^{(N)}`$ depend on $`N`$ and $`lim_N\mathrm{}z_i^{(N)}=z_i`$ with $`|z_i|<1`$, except the unit eigenvalue corresponding to the equilibrium density.
2. These eigenvalues determine the decay of correlations in the chaotic component of phase space.
3. In the chaotic component of phase space the right eigenfunctions of the truncated Frobenius-Perron operator $`\widehat{U}^{(N)}`$ are localized on the unstable manifold of $`𝐅`$, and tend to be uniform along this manifold. The left eigenfunctions are localized on the stable manifold. These manifolds are not as simple as for the baker map.
4. The eigenvalues of $`\widehat{U}^{(N)}`$ corresponding to eigenfunctions localized in regular regions approach the unit circle as $`N`$ or $`l_{max}`$ increases.

## V The Frobenius-Perron Operator for the Kicked Rotor

In this section the decay of correlations in time will be studied for the kicked rotor, defined by the Hamiltonian (7), and its dynamics is given by the standard map (10) . The phase space is presented in Fig. 1. It is a mixed system (where in some regions of phase space the motion is chaotic and in other regions it is regular). In this section we summarize results that were presented in  (see also ,  and ).

For sufficiently large values of the stochasticity parameter $`K`$ (that takes a typical value), the spread of angular momentum is diffusive, to a good approximation, namely

$$(J_nJ_0)^22Dn,$$
(69)
for large $`n`$. The average is over the initial points and is denoted by $`\mathrm{}..`$. A very crude way to get this result is to iterate (10) $`n`$ times to obtain

$$J_nJ_0=K\underset{i=1}{\overset{n}{}}\mathrm{sin}\theta _i,$$
(70)
then squaring and averaging with the assumption of absence of angular correlations

$$\mathrm{sin}\theta _i\mathrm{sin}\theta _j=\frac{1}{2}\delta _{ij},$$
(71)
one obtains (69) with $`D=\frac{1}{4}K^2`$. A more careful calculation, that takes some of the angular correlations into account results in

$$D(K)=\frac{K^2}{4}\left(12J_2(K)+\mathrm{}\text{}\right).$$
(72)
It is actually an expansion in powers of $`1/\sqrt{K}`$. Because of the mixed nature of phase space it is expected that diffusion may not take place asymptotically in time . Moreover, there are values of $`K`$ (near integer multiples of $`2\pi `$) where acceleration in momentum is found for some initial conditions, resulting in $`(J_nJ_0)^2n^2`$. In what follows the Frobenius-Perron operator will be used to study whether correlations in angular momentum decay as expected for a true diffusion process, described in the Introduction (see (20)-(27)). We will analyze also the decay of angular correlations.

The Frobenius-Perron operator, corresponding to the map (10) on the torus

$$(0J<2\pi s)$$
(73)
$`(0\theta <2\pi ),`$

where $`s`$ is integer, will be studied. For the usual kicked rotor $`s\mathrm{}`$. Noise is added to the standard map by the addition of a random variable $`\xi _n`$ with variance $`\sigma ^2`$ on the right hand side of the equation for $`\theta _{n+1}`$ in (10) following . The noise leads to truncation of the Frobenius-Perron operator. The natural basis for the analysis is

$$J\theta |km=\frac{1}{\sqrt{2\pi }}\frac{1}{\sqrt{2\pi s}}\mathrm{exp}(im\theta )\mathrm{exp}\left(i\frac{kJ}{s}\right).$$
(74)
where $`k`$ and $`m`$ are integers. Note that the functions $`J\theta |k,m=0`$ form the basis of eigenstates of the diffusion operator in the angular momentum $`J`$. In this basis, in presence of noise, the Frobenius-Perron operator is

$$km|\widehat{U}^{(\sigma )}|k^{}m^{}=J_{mm^{}}\left(\frac{k^{}K}{s}\right)\mathrm{exp}\left(\frac{\sigma ^2}{2}m^2\right)\delta _{kk^{},ms}.$$
(75)
For $`\sigma 0`$ the operator is not unitary. The effective truncation is $`m<1/\sigma ,|kk^{}|<s/\sigma ,|mm^{}|<k^{}K/s`$. The $`\sigma 0`$ limit , corresponding to the $`N\mathrm{}`$ limit of the previous sections, will be taken in the end of calculation. The analysis will be performed for

$`s`$ $``$ $`1`$ (76)
$`K`$ $``$ $`1`$ (77)
$`{\displaystyle \frac{kK}{s}}`$ $``$ $`1`$ (78)
and the limit $`\sigma 0`$ will be taken in the end of the calculation.

The Ruelle resonances are identified from the matrix elements of the resolvent (37). It is useful to define a variant of the resolvent

$$\widehat{R}^{}(z)=\underset{j=0}{\overset{\mathrm{}}{}}\widehat{U}^jz^j=\frac{1}{1z\widehat{U}}$$
(79)
that is related to $`\widehat{R}(z)`$ via

$$\frac{1}{z}\widehat{R}^{}\left(\frac{1}{z}\right)=\widehat{R}(z).$$
(80)
Also the matrix elements of $`\widehat{R}`$ and $`\widehat{R}^{}`$ satisfy (80). Let

$$R_{12}=k_1m_1|\widehat{R}(z)|k_2m_2$$
(81)
and

$$R_{12}^{^{}}=k_1m_1|\widehat{R}^{}(z)|k_2m_2$$
(82)
be the matrix elements of $`\widehat{R}`$ and $`\widehat{R}^{}`$ respectively. Therefore a singularity of $`R_{12}`$ at $`z_c`$ implies a singularity of $`R_{12}^{^{}}`$ at $`1/z_c`$ and vice versa. The series (37) for $`\widehat{R}(z)`$ converge for $`|z|>1`$ while the series (79) for $`\widehat{R}^{}(z)`$ converge for $`|z|<1`$. Therefore the Ruelle-Pollicott resonance (located inside) that is closest to the unit circle corresponds to the singularity of $`R_{12}^{^{}}`$ that is closest to the unit circle (located outside). This singularity is the radius of convergence of the series

$$R_{12}^{^{}}=\underset{j=0}{\overset{\mathrm{}}{}}a_jz^j,$$
(83)
where

$$a_j=k_1m_1|\widehat{U}^j|k_2m_2.$$
(84)
According to the Cauchy-Hadamard theorem (see ) the inverse of the radius of convergence is given by

$$r^1=\underset{j\mathrm{}}{lim}sup\sqrt[j]{|a_j|}.$$
(85)
and asymptotically

$$|a_j|\frac{\text{const.}}{r^j}.$$
(86)
Since the radius of convergence is the singularity of $`R_{12}^{}(z)`$ that is the closest to the unit circle it satisfies $`r=1/z_c`$ and $`\sqrt[j]{a_j}z_c`$ in the limit $`j\mathrm{}`$ (at least for some subsequence of $`\{j\}`$).

The limit of the series (84) was found in the leading order in perturbation theory in the small parameters implied by (76). The calculation can be performed separately for the various values of $`k`$, resulting in

$$z_k=\mathrm{exp}\left(\frac{k^2K^2}{4s^2}\left(12J_2(K)e^{\sigma ^2}\text{}\right)\right),$$
(87)
corresponding to diffusion modes in presence of noise. Taking the limit $`\sigma 0`$ one finds

$$z_k=e^{k^2D(K)/s^2}$$
(88)
where $`D(K)`$ is given by (72). This relaxation is similar to the one that is found for usual diffusion (see (20-27)). The relaxation rates are

$$\gamma _k=\frac{k^2}{s^2}D(K),$$
(89)
that correspond to (24). These describe the decay of correlations in angular momentum $`J`$. The equilibrium density corresponds to $`z_0`$. Exploring the subspace involving $`|k=0,m`$, one finds that the Ruelle-Pollicott resonances, corresponding to the slowest mode of the relaxation of correlations in angle, take four values

$$\pm \stackrel{~}{z},\pm i\stackrel{~}{z}$$
(90)
with

$$\stackrel{~}{z}=\sqrt{\left|J_{2m^{}}(m^{}K)\right|\mathrm{exp}\left(\sigma ^2m^2/2\right)}.$$
(91)
Taking the limit $`\sigma 0`$ one finds

$$\stackrel{~}{z}=\sqrt{\left|J_{2m^{}}(m^{}K)\right|},$$
(92)
where $`m^{}`$ is the integer $`m`$ that maximizes $`\left|J_{2m}(mK)\right|`$ for a given value of $`K`$. A similar result was found in  by a variational approach. The corresponding relaxation rate is

$$\stackrel{~}{\gamma }=\mathrm{ln}\stackrel{~}{z}.$$
(93)
In order to test the analytical results, the correlation function of the form (35)

$$C_{fg}(n)=f|\widehat{U}^n|g,$$
(94)
where $`f`$ and $`g`$ are taken to be basis states of (74), was computed numerically. The results are plotted in Fig. 9.

Clear exponential decay is found. From the slopes, the relaxation rates $`\gamma _k`$ are obtained for the various modes, $`k`$. The values of $`D`$ implied by $`\gamma _k`$ are calculated from (89) and compared to the ones obtained from (72) in Fig. 10.

Good agreement was found. In Figs. 9 and 10 large values of the stochasticity parameter, $`K10`$, were used. In this regime the regular regions are extremely small (they are invisible in Fig. 1), therefore the theory presented in this review is expected to work very well, as is indeed found in Fig. 10. In Fig. 11 the correlation function (94) is presented for $`K20`$.

From these results the relaxation rates $`\gamma _k`$ were calculated. The resulting values of the diffusion coefficient obtained from (89) are compared to the ones found from (72) in Fig. 12, for $`K20`$. Here, in contrast to the regime of large $`K`$ presented in Fig. 10, appreciable deviations are found.

These are related to sticking to regular structures, that is not taken into account in the theory presented in this review. To test the decay of correlations in the angle variable, the correlation function (94) was computed for $`|f`$ and $`|g`$ that are basis states of the form $`|k=0,m`$, and some of the results are presented in Fig. 13.

The relaxation is very fast, therefore it is very difficult to estimate numerically the relaxation rate. The numerical results are compared to the analytical prediction (93) in Fig. 14.

The analytical formula provides a reasonable estimate of the relaxation rate and of its dependence on the stochasticity parameter $`K`$. There are however deviations that are significant. Also here some of the most significant deviations were found to be related to regular structures.

The main results that were obtained for the kicked rotor are:

1. Addition of noise results in the effective truncation of the Frobenius-Perron operator. In the limit of vanishing noise the Ruelle-Pollicott resonances are found. The resonances depend on the noise, and therefore on the truncation, as is the case for the kicked top, and approach limiting values inside the unit circle in the limit of vanishing noise, except the unit eigenvalue corresponding to the equilibrium density.
2. The Ruelle-Pollicott resonances determine the decay of correlations.
3. Deviations from exponential decay of correlations may be found after some time as a result of the existence of regular structures.
4. The following physical picture emerges. First the correlations in the angle variable $`\theta `$ decay. Then the correlations in angular momentum $`J`$ decay as for a true diffusion process. For finite $`s`$ a uniform equilibrium density is reached. For $`s=\mathrm{}`$, after the decay of the angular correlations, diffusive spreading takes place.
5. At some time, that is increasing with the stochasticity parameter $`K`$, sticking to regular structures becomes important, and the picture based on Ruelle-Pollicott resonances and exponential decay of correlations, breaks down.

## VI Summary

It was demonstrated that Ruelle-Pollicott resonances are relevant even when the conditions for the Ruelle-Pollicott theorem do not hold. A method for the calculation of these resonances by a finite truncation of the Frobenius-Perron operator was presented. The relation between the resonances and the decay of correlations was established under these conditions. For mixed systems (as well as for other nonhyperbolic systems) the picture breaks down after some time, that may be very long, as a result of sticking to regular structures.

In order to evaluate the relevance of the Ruelle-Pollicott resonance picture for realistic systems it is instructive to introduce several time scales:

1. If $`\gamma _1`$ is the slowest decay rate (corresponding to the resonance that is closest to the unit circle), $`t_{chaos}=1/\gamma _1`$ is the decay time of correlations and the relaxation time to the invariant density.
2. The time scale when sticking to regular structures becomes important will be denoted by $`t^{}`$.
3. In presence of external noise correlations are destroyed on the time scale $`t_c1/\sigma ^2`$, where $`\sigma ^2`$ is the variance of the noise.

In order to observe the exponential decay of correlations it is required that

$$t_{chaos}t_c,t^{}.$$
(95)
Destruction of sticking to regular structures requires $`t_c<t^{}`$. In the regime

$$t_{chaos}t_c<t^{}$$
(96)
the decay of correlations is expected to be similar to the one of hyperbolic systems and one should be able to explore it in the framework of the method that was outlined in the present review. Application of the method to specific examples is of great interest.

###### Acknowledgements.

It is our great pleasure to thank Fritz Haake for many illuminating discussions during the preparation of the lectures and of the review, for the files of the figures from  and for the permission to use them. Part of the review is based on work in collaboration with Oded Agam and Maxim Khodas, that is acknowledged with great pleasure. We would like to thank Christopher Manderfeld and Joachim Weber for useful discussions and communications. SF thanks Andreas Buchleitner for the hospitality at the Max Planck Institute for the Physics of Complex Systems in Dresden, where the lectures were prepared. This research was supported in part by the US-Israel Binational Science Foundation (BSF), by the Minerva Center of Nonlinear Physics of Complex Systems, by the Max Planck Institute for the Physics of Complex Systems in Dresden, and by the fund for Promotion of Research at the Technion.

Most of the results presented in the review can be found in references marked by \**


format: text

text:
# SOLITON ORGANIZATION of THERMAL FIELD in a CHAIN at HIGH TEMPERATURE

## Introduction

A phonon (the extended mode) delocalized in the coordinate space, but localized in the wave space is the fundamental quasiparticle in crystal reality. Since any of a number of phonons occupies the entire space of the crystal, each individual particle takes part in the motion of all the phonons. In the result of accidental phase combination of some phonon groups the multiphase and single-phase waves may arise with spatial localization of the order of the lattice spacing. In consequence of disperse features of the lattice they will be scattered prior to their wave properties will be exhibited. With nonlinearity phonon groups may form a long-lived bound states and in this situation the wave properties will be exhibited in full measure. The soliton-like excitations of two types, such as kinks and bell-form solitons \- , on the one hand, or breasers (gap solitons, intrinsic localized modes and so on) , on the other hand, observed in computer experiments for different nonlinear lattices, directly confirm this fact. The former may exist below the linear phonon band, the latter ones above the upper harmonic phonon band edge, that is inside the gap  or inside the self-induced gap .

Nonlinear waves of different nature manifest a series of common properties including suppression and absorption of small waves by larger ones. In Ref.  such phenomenon was observed for small and large exact breasers. In Ref.  the chaotic breasers were observed in the numerical experiment for the FPU-chains, which manifest the effect of energy concentration by separate excitations. As the initial condition the $`\pi `$-mode was taken. It is shown that the $`\pi `$-mode decays into local excitations with chaotic behavior. In the initial stage the concentration of energy by separate breasers grows. Latter on these excitations vanish and the system approaches the equipartition state. The authors state that the ”pumping” of energy from the high-frequency $`\pi `$-mode to the low-frequency region of phonon spectrum takes place by this way.

After relaxation the chaotic breasers die out because of ”starvation”. This brings up the question: is such behavior of a system, when the chaotic breasers do not die out and rest after the relaxation, possible? To achieve such regime it is necessary, that the mean energy per one particle, that is, temperature, be high enough. The localization modes are regarded also for FPU-chains in numerical experiment of Ref. .

The aim of this paper is to investigate similar effects for kink solitons at high temperature and their influence on thermal field by means of computer simulation. To minimize the non-equilibrium pumping effect here the ”white noise” distribution is used as the initial state. The ”white noise” is energy distribution closest to thermodynamic equilibrium state. One may expect, that the relaxation passes fast enough.

## 1. Basic Statements

One may note that a FPU-chain is a very idealistic model for atomic systems and it deviates from realistic ones to a greater extend with the vibration amplitude increase. In our case such a situation presents the main interest for the investigation. It follows herefrom that one of the realistic potentials containing non-linearity of all orders must be used for such system at high temperature.

All simulation was performed for a chain built up from 200 particles interacting via the most popular 6-12 Lennard-Jones’ potential, that is:

(1.1)
$$U_{ij}=u_{min}\{(r_0/r_{ij})^{12}2(r_0/r_{ij})^6\}$$

where $`r_{ij}=\sqrt{(X_iX_j)^2+(Y_iY_j)^2}`$\- is the distance between particles with numbers $`i`$ and $`j`$, and Cartesian coordinates $`X_i,Y_i,X_j,Y_j,`$ $`r_0`$ \- - is the equilibrium distance for two-particle systems, $`u_{min}`$ \- is the binding energy. For the sake of simplicity the reduced units of length and time as $`\stackrel{~}{r_{ij}}=r_{ij}/r_0`$ and $`\stackrel{~}{t}=t/T`$, respectively, are defined. Here $`T=2\pi \sqrt{m/U_{ij}^{\prime \prime }(0)}`$ is the period of small vibrations, $`m`$ is the particle mass which equals unity thereafter. The binding energy in the reduced units turned out to be equal to $`2(\pi /6)^2`$. A smaller than $`T`$ conventional unit of time must be introduced for the description of wave phase. We take it to be equal to $`t_{con}=0.005`$ (time step in the computer experiment). Sixth-order Yoshida’s symplectic algorithm is used for calculation . The law of conservation of energy is fulfilled to relative accuracy $`10^6`$.

The equilibrium state is determined by coordinates $`X_i^{(0)},Y_i^{(0)}`$:

(1.2)
$$X_i^{(0)}=r_0i,Y_i^{(0)}=const$$

The initial distribution of particle velocities is chosen in the form of white noise with amplitude 4.078 for X- polarization only. The initial disturbances along Y-direction was not preset, however, for additional control over computer errors the movement in this direction was permitted. The average additional energy per one particle contributed by white noise turned out to be equal to 2.6895 of reduced units. This value is 4.9 times as much as the binding energy. Such correlation of energies corresponds to high enough temperature of the chain which is much higher than the melting point.

## 2. Results and Discussions

The evolution of white noise in time for a chain built up from 200 particles is presented by four graphic fragments. Every fragment is composed of 46 constitutive spatial profiles reflecting the square of particle velocity or kinetic energy versus its number in the chain. In its turn, every constitutive profile is formed by graphic superposition of ten elementary profiles taken during ten time steps. Such display permits one to improve the tracing of waves for a long time.

Really as would be expected, the early stage of evolution (the upper fragment in fig.1) is very close to homogeneous distribution. Nevertheless at this stage the local excitations with size of the order of lattice spacing are traced already. These local objects conserve their shape both at free spreading, and after interactions in between for a long time. Their velocities determined by slope of traces vary from 10 to 16 of reduced units and they turn out to be proportional to wave amplitude. A consequence of the above is the fact that we deal with the soliton-like waves.

From general considerations it is quite evident that they are the soliton-like objects of kink kind. Really, if spreading of excitation is presented as consecutive frontal collisions of particles, then in the result of a separate collision a moving particle stops and the resting particle starts to move with the same velocity. In the next collision this moving particle stops too in the new position, and so on. As a result, we have the switching wave of the kink type. Evidently, in our case we have a multitude of such waves, interacting in addition each with other.

By virtue of white noise proximity to the equilibrium state if may seem that the relaxation will pass sufficiently fast. However, the next, second fragment on the top of fig. 1 testifies that the system approaches equipartition state, but deviated towards even greater inhomogeneity of the spatial distribution of particle velocities. In the chain the high-amplitude waves arise, concentrating a significant part of system energy. As the law of energy conservation is fulfilled with high precision, this can be traced at the cost of energy re-distribution among the excitations, namely, high excitations have captured the energy of the weak ones. This phenomenon is analogous to the decay of the For calculation the $`\pi `$-mode into the chaotic breasers . The surprising thing is that almost identical effect accompanies the breakdown of white noise. Next two fragments (fig. 1) show, that the level of inhomogeneity of the velocity field varies in time, however in any event it is considerably above than at the starting stage (the uppermost fragment). This inhomogeneity don’t vanish at high time of system evolution and, it is obvious that it don’t vanish at all. Thus, while in the Ref.  chaotic breasers were an interface at system relaxation and they vanished after the achievement of thermal equilibrium, in our case the soliton-like excitations are the main objects of the equilibrium thermal field. One may conclude that solition organization of thermal field takes place in such situation.

To support this conclusion let us correlate calculated distribution of particle energy with theoretical Gibbs one. In fig. 2 the curves with circles correspond to ”real” data of the numerical experiments, solid ones correspond to theoretical distribution: Descending curves describe a distribution of particles by energy , determined as number of particles with energy $`(i1)\mathrm{}\epsilon \epsilon <i\mathrm{}\epsilon `$, where $`i`$ \- is the number of interval. The curves with one or many peaks are the energy spectra , determined as total energy of particles in the same interval. Theoretical equilibrium distribution of this quantities are of the form:

$`n_i=n_0\mathrm{exp}(\epsilon _i/\theta ),`$
(2.1) $`E_i=\epsilon _in_i`$
with the fulfillment of laws of conservation for particles and energy:

$`N={\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}n_i,`$
(2.2) $`E={\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}\epsilon _in_i`$
where $`n_0`$ \- is the number of particles with minimum energy, $`\theta `$ \- is temperature in the energy representation of reduced units.

From the last relations the expressions for $`n_0`$ and $`\theta `$ follow:

$`\theta =E/N,`$
(2.3) $`n_0=N^2/E`$
All graphs built with numerical experiments are obtained by means of averaging over 460 time steps and they correspond to fragments in fig 1. Energy interval $`\mathrm{}\epsilon `$ used at building of graphics was taken to be equal to 1.6 of reduced units.

One can see from the first graph that the white noise deviates from Gibbs distribution to a greater extend than any in fig. 2. On the energy spectra for real data peak has higher value than for theoretical ones at the cost of a fraction of low-energy excitations. The particles with high energy are absent in the system at all. With the evolution of white noise the energy is ”pumped” from the low-energy excitation to the high-energy ones. One can see from the second graph in fig. 2 that the first peak has been decreased below theoretical value while additional peaks have arisen in the high-energy region. In this region the energy curve lies markedly higher than the theoretical one. With latter evolution (the third graph in fig. 2) a peculiar kind of returning is observed, namely, real distribution achieves the theoretical one. Additional peaks have vanished. It should be noted that such equilibrium differs in the nature of its organization from a phonon variant . Really, in our case (see the third fragment in fig. 1) local excitations continue to exist during all the time. This is the principal new finding.

One may expect that the system having achieved, finally, the thermal equilibrium, will rest in this state for a long time. However, the further evolution of the system (see the fourth graph in fig. 2) shows that with time the system returns to the state with additional peaks again. At the same time the returning can’t be considered as simple fluctuation, since it develops directly and sequentially during a long time of about 300-500 reduced units of time. Latter on both of these states alternate in time.

The departure of equilibrium distribution from the Gibbs one is marked in Ref.  for chock waves in one-dimensional chains. If $`\theta _0`$ is the initial temperature and $`\theta `$ is the final temperature ahead and behind of the chock front, respectively, then the equilibrium distribution of particles velocities is determined by the following expression :

(2.4)
$$n_i=\frac{1}{2}(2k\theta )^{\frac{1}{2}}[\mathrm{exp}(\frac{(\nu _i\vartheta ^{\frac{1}{2}})^2}{2\theta _0})+\mathrm{exp}(\frac{(\nu _i+\vartheta ^{\frac{1}{2}})^2}{2\theta _0})]$$

where $`\vartheta =\theta \theta _0`$, $`\nu _i`$ \- is velocity of a particle, or with regard of potential energy contribution:

(2.5)
$$n_i=\overline{n}_0\mathrm{exp}(\frac{\epsilon _i}{\theta _0})\mathrm{coth}(2\epsilon _i\theta )^{\frac{1}{2}}$$

where $`\overline{n}_0=n_0\mathrm{exp}(\frac{\vartheta }{2\theta _0})`$ , $`n_0=(2\pi \theta _0)^{\frac{1}{2}}`$ \- is the same $`n_0`$ as in the expression (2).

One may see, that at $`\vartheta >0`$ $`\overline{n}_0<n_0`$ is true, that is, the number of particles with minimum energy is less then it for the Gibbs distribution. In the same time owing to factor $`\mathrm{coth}(2\epsilon _i\theta )^{\frac{1}{2}}`$ a share of the high-energy particles grows. The same is true for energy spectrum.

From third graph in fig. 2 one may see, that the ”experimental” distribution being qualitatively described by the expression (2.5) is formed in the system. However no peaks in the spectra in fig. 2b, d are described by this expression. The emergence of the ones is quantitatively new peculiarity of a non-linear system at high temperature. Its emergence obviously connects with the long-living local high-energy excitations.

Thus the described computer experiment has shown that the fundamental soliton organization of thermal wave field is possible at high temperature. Such phenomenon has led to spontaneous production of intensive nonlinear waves, which can be considered as a peculiar kind of self-organization of a thermal field. Investigation of this phenomenon is important to gain a better understanding of the break of solid.


text:
# The extreme ultraviolet excess emission in five clusters of galaxies revisited

## 1 Introduction

The Extreme Ultraviolet Explorer (EUVE) has detected emission from a few clusters of galaxies in the $``$70-200 eV energy range. By order of discovery, these were: Virgo (Lieu et al. 1996b; Berghöfer et al. 2000a), Coma (Lieu et al. 1996a), Abell 1795 (Mittaz et al. 1998), Abell 2199 (Bowyer et al. 1998), Abell 4059 (Berghöfer et al. 2000b) and Fornax (Bowyer et al. 2001). The problem of the physical origin of the EUV emission was immediately raised. Excess EUV emission relative to the extrapolation of the X-ray emission to the EUVE energy range was detected in several clusters, suggesting that thermal bremsstrahlung from the hot ($`10^8`$K) gas could not be entirely responsible for the EUV emission. Although there has been some controversy on the reality of this EUV emission, it now seems well established for the following clusters: Coma (Lieu et al. 1996a), Virgo (Lieu et al. 1996b; Berghöfer et al. 2000a; Bonamente et al. 2001a), Abell 1795 (Mittaz et al. 1998; Bonamente et al. 2001a), Abell 2199 (Bowyer et al. 1998; Kaastra et al. 1999; Lieu et al. 1999b) and Abell 4038 (Bowyer et al. 1998).

This excess EUV emission can be interpreted as due to two different mechanisms: thermal radiation from a warm ($`10^510^6`$K) gas, as first suggested by Lieu et al. (1996a), or inverse Compton emission of relativistic electrons either on the cosmic microwave background or on stellar light originating in galaxies (Hwang 1997; Bowyer & Berghöfer 1998; Ensslin & Biermann 1998; Sarazin & Lieu 1998; Ensslin et al. 1999), or both. Note that although the existence of a multiphase intracluster medium is plausible (Lieu et al. 2000; Bonamente et al. 2001b), there are several difficulties with the thermal model. The most serious one is that since gas in the temperature range $`10^510^6`$K cools very rapidly, a source of heating for this gas is necessary, and only one model has been proposed to account for this gas heating (Fabian 1997). Besides, the non-detection by the Far Ultraviolet Spectroscopic Explorer of the OVI $`\lambda \lambda `$1032,1038 lines in the spectra of the Coma and Virgo clusters could exclude the presence of warm gas at a 2$`\sigma `$ level (Dixon et al. 2001). On the other hand, non thermal models (Sarazin & Lieu 1998; Bowyer & Berghöfer 1998; Lieu et al. 1999a; Ensslin et al. 1999; Atoyan & Völk 2000; Brunetti et al. 2001) seem able to account for the excess EUV emission, but require very high cosmic ray pressure.

Following a controversy that the extended EUV emission from some clusters could be simply an artefact of bad background subtraction (Bowyer et al. 1999; Lieu et al. 1999b) we have performed a wavelet analysis and reconstruction of the EUVE images of several clusters. Our aims were to confirm the detection of EUV emission as far as possible from the cluster center, to derive accurate EUV and X-ray profiles and to confirm the existence and radial distributions of the soft excesses over thermal bremsstrahlung in these clusters.

We present here the wavelet reconstruction of the images of these five clusters, both in the EUV (EUVE images) and in X-rays (ROSAT PSPC images). Abell 4038 was also analyzed but was discarded due to its low signal to noise ratio in the EUV. We then compare the EUV and X-ray profiles of the images thus obtained, arguing that if both emissions have the same physical origin it should not be possible to distinguish the shapes of these profiles. In fact, we will see that they are found to be different in all clusters, and that soft excesses over thermal bremsstrahlung are detected in all five clusters with various radial distributions.

The data and method are presented in Sect. 2. Results are presented in Sect. 3 and discussed in Sect. 4.

## 2 The data and image analysis

### 2.1 The data

Five clusters were observed in the EUV range with the EUVE satellite and in X-rays with the ROSAT PSPC satellite (archive data). The exposure times and main cluster characteristics are given in Table 1.

In order to have comparable spatial resolutions in the EUV and X-rays, the EUV images were rebinned 4$`\times `$4, leading to a “superpixel” size of 0.3077 arcminutes (18.5 arcseconds). X-ray images were also rebinned to the same 18.5 arcsec pixel size in order to allow a direct comparison of the EUV and X-ray profiles.

The linear scale per superpixel at the cluster distance, estimated with H$`{}_{0}{}^{}=50`$ km s-1 Mpc-1 and q0=0 is given in Col. 7 of Table 1.

### 2.2 The EUVE data

The EUVE/Deep Survey photometer (DS) images were extracted from the raw event lists with the customary linear bin size of 13 pix/arcmin. The rectangular shape of the Lex/B filter ($``$70-200 eV) results in images where the length of one side of the image greatly exceeds the other side (see Fig. 1). When more than one observation was available for each target, images were coadded to improve the S/N; the rectangular shape often results in only partial overlap between the exposures leading to diamond-shaped geometries (see Fig. 2). Without taking into account any edge corrections, this peculiar shape of the EUVE images may be a problem for any analysis involving an isotropic distribution of the data at scales larger than about one-third the size of the shorter side.

Providing that the image does not truncate the cluster emission in any direction, this shape problem can be treated as a classical edge effect. A simple solution is then to use the asymptotic background values to obtain a simulated larger square image. From a practical point of view only a small truncation of the data occurs in our images along the shorter axis. Its impact on an azimuthal profile can therefore be considered as negligible and we thus decided to apply the above procedure to modify the geometry of the EUVE images. We now give details on the required background estimation.

Fig. 1 shows the rebinned image of a long observation of the Virgo cluster. The cluster is clearly visible in the central area of the field of view. Three rectangles are superimposed on the raw image. The largest one, denoted by (A,B,C,D) follows the outer detector shape closely. Its edges are detached from the detector boundaries to avoid the steepening of the background near the edge, which would bias the noise estimation (see the figure 2 of Lieu et al. 1999c). In this box, we define two rectangles (denoted $`r_1`$ and $`r_2`$), one towards each edge, as far away from the central source as possible. Both $`r_1`$ and $`r_2`$ are assumed to contain only noise and to be representative of the noise level far from any celestial source. This is ensured because the centers of both rectangles lie at $`40`$ arcmin of the central source, a distance notably greater than the detection radius of the extended EUVE cluster emission, as seen with classical methods. The second order moment of the count distribution is used to estimate the level of the noise in each of these rectangles. The noise level in the large box is then computed as the mean of the noise levels in both rectangles. Once this has been done, each pixel outside the larger box is replaced by a random number drawn from a Poisson distribution, the variance of which is obtained by adding in quadrature the variances in each of the rectangles. The accuracy of the noise estimation is readily observable in the product image: if the noise is underestimated, the central box will obviously rise against the rest of the image, with its boundaries showing a clear discontinuity. On the other hand, if the noise has been overestimated, the outside parts of the new image will rise against the central box, again showing a discontinuity on the boundaries.

When more than one EUVE observation of a cluster was available, we used the co-added images in order to increase the signal to noise ratio. Fig. 2 illustrates the case of Abell 1795. This image is typical of the co-added images we dealt with in the case of other clusters. The difference in roll-angle of the satellite in both observations is apparent. The central intersection of both observations obviously has the highest signal to noise ratio. It corresponds to a region indicated by the circle of radius 16 arcmin in Fig. 3, which we will take as the radial limit for our analysis. For each of the separate observations, we drew again a box closely following the outer detector shape. Both boxes are labelled (A,B,C,D) and (A,B,C,D) in Fig. 2. In each of the boxes, we define two small rectangles, where the noise is evaluated. The noise level in a box is computed as the mean of the noise levels in both rectangles. Once this has been done for each box, each pixel outside the intersection of the boxes (the central diamond) is again replaced by a random number drawn from a Poisson distribution, as discussed before. Once again, the under or overestimation of the noise will be clearly visible. An example of such an image is shown in Fig. 3.

### 2.3 The ROSAT data

The ROSAT X-ray data for each cluster were obtained from the HEASARC web archive111http://heasarc.gsfc.nasa.gov/. Owing to the higher sensitivity, we preferred PSPC to HRI data and pointed observations to all-sky survey scanning mode observations. When multiple PSPC observations were available, we always used the one with the smallest offset radius compared to the center of the cluster, and then the longest one.

S. Snowden’s Extended Source Analysis Software (ESAS222ftp://legacy.gsfc.nasa.gov/rosat/software/fortran/sxrb) was used to perform the data reduction, i.e., rejection of high-background times, modelling of the different background types, energy-dependent background subtraction, exposure and vignetting corrections (see Snowden et al. 1994). We adopted a conservative value of $`170\mathrm{cts}\mathrm{s}^1`$ for the maximum Average Master Veto rate allowed. We then carefully examined the light curves of the total counts in the entire image per energy band, and checked that no short time scale glitches were present. We limited our analysis to the \[0.5,2.0\] keV energy band (bands R4-R7 as defined by Snowden). This preliminary reduction produced a surface brightness image with $`512\times 512`$ pixels, 15 arcsec per pixel (roughly the FWHM of the PSPC PSF at 1 keV), as well as an exposure image for each of the bands and a background image. The image was then divided by the mean exposure map. We thus obtained a $`512\times 512`$ pixel image, in units of counts s-1 pixel-2, which was then rebinned to 18.5$`\times `$18.5 arcsec2 pixels to match the EUV superpixel size.

The ROSAT images do not present the same problems as the EUVE ones, because in this case, the geometry is circular around the central cluster (all our observations were pointed). The only geometrical constraint arises from the presence at an offset angle of approximately 20 arcmin of the structure support of the telescope. However, since this value is smaller than the EUVE limiting radius, it does not place any further constraint here. We therefore just extracted a square image inscribed in the 20 arcmin radius circle.

### 2.4 Wavelet analysis and reconstruction

The comparison between the EUVE observed surface brightness profiles of clusters and their ROSAT X-ray predicted counterparts is the key step in our analysis. However a direct comparison of the X-ray and EUVE data would provide misleading results due to the noise in the raw data. To remove most of this noise while preserving information on the astronomical emission at all spatial scales of interest we decided to apply part of the wavelet vision model described in detail in Rué & Bijaoui (1997).

Let us hereafter briefly summarize the main steps of the image processing involved. First, a discrete wavelet transform of the image $`I_0`$ is performed using the so-called à trous implementation, which is nothing more than the classical multiresolution algorithm of Mallat (1989) but without any decimations since “holes” are introduced in the required convolutions. The initial sampling is kept for each scale; this insures shift-invariance and allows straightforward interscale computations. The scaling function leading through a dilation factor $`a`$ to the involved sets of low-pass and band-pass (wavelets) filters is a cubic B-spline (defined in the Fourier space as sinc$`{}_{}{}^{4}(\omega )`$), very well suited for enlightening gaussian-like patterns of various sizes in an image. Owing to the properties of the underlying convolutions a stationary white noise signal leads to a statistical distribution of the wavelet coefficients with a zero mean value at each computed scale $`a_j`$ and a standard deviation decreasing as the spatial scale increases. Structures in the transform which are related to objects in direct space can then be identified and located by selecting high enough wavelet coefficients at each scale.

Removing noise from the data without any loss in resolution can thus be achieved quite easily in wavelet space by keeping only those coefficient values with a probability to be due to a chance fluctuation of the noise process lower than a chosen decision level. This second step of the noise removal process requires the knowledge of the probability density function of the coefficients at scale $`a_1=1`$ pixel (here 18.5 arcsec) for the noise process and of the decreasing law of the rms scatter $`\sigma `$, the asymptotic value of which appears to be 0.5 for a bidimensional image $`I_0`$ and a dyadic scale analysis ($`a_j=2^{j1},j𝒵`$) : $`\sigma _j=\sigma _1/2^{j1}`$. In case of Gaussian noise, the probability density function is easily obtained from a histogram of the wavelet coefficients providing that this first scale of one pixel in size is noise dominated. Due to the stationarity of the Gaussian process, the selection of the so-called statistically significant coefficients $`\nu _j`$ at scale $`a_j`$ can then be done using the same k-$`\sigma _j`$ threshold. We end up in this way with a subset $`\nu ^{(0)}`$ of coefficients at different scales, the pixels of which define what is called hereafter the multiscale significant mask $`M_0`$. In case of Poisson noise stationarity can be obtained by applying an Anscombe (1948) transform on the data : $`f(x)=2\sqrt{x+0.375}`$. When the mean parameter $`\mu `$ of the Poisson process is greater than about 5, this transform indeed leads to a new process with a very small bias in the expectation value ($`\frac{<f(x)>f(\mu )}{f(\mu )}<2\%`$) and a nearly perfect variance stabilization. The same decision scheme as in the Gaussian case can therefore be used for defining the multiscale significant mask from the wavelet transform of this image with a constant variance. Only second order errors in the wavelet coefficients would arise at the very first scales in areas with very low counts. This is not a matter of concern for the present application since the computation of a correct azimuthal brightness profile necessitates only a proper account of the medium and large-scale features in these faint regions.

The third and final step is to restore a positive image $`I`$ where the noise is strongly reduced from the subset $`\nu ^{(0)}`$ of the significant wavelet coefficients of $`I^{(0)}I_0`$. This is an inverse problem which has in general an infinity of solutions and the choice among these solutions must be done with respect to a regularization constraint. To do so, an iterative conjugate gradient algorithm is used. This algorithm minimizes the distance $`\nu ^{(0)}\nu ^{(n)}`$ where $`\nu ^{(n)}`$ is the restriction to $`M_0`$ of the wavelet transform of the solution $`I^{(n)}`$ at step $`n`$. Stating A the operator associated to the wavelet transform followed by the projection onto $`M_0`$ and Ã the joint operator, $`I^{(n+1)}=`$ Ã$`(\nu ^{(n)})`$ is obtained from the previous solution by : $`I^{(n+1)}=I^{(n)}+\alpha ^{(n)}I_r^{(n)}`$ where the residual image $`I_r^{(n)}`$ is defined as Ã$`(\nu _r^{(n)})+\beta ^{(n)}I_r^{(n1)}`$ with $`\nu _r^{(n)}=\nu ^{(0)}`$ A$`(I^{(n)})`$ ; $`\alpha `$ and $`\beta `$ are convergence parameters. Convergence ($`\frac{\nu ^{(0)}\nu ^{(n)}}{\nu ^{(0)}}<1\%`$) is reached after very few steps, leading to an image $`I`$ without any discontinuities despite the hard thresholding performed in the wavelet space. Note also that any flat instrumental or astronomical component or gradient in the initial image is automatically removed in image $`I`$ taking benefit of the mathematical properties of the wavelet function related to our scaling function.

We applied such an image processing to both our EUVE and ROSAT square images considering seven or eight spatial scales for the wavelet transform according to the size of the image. Pure Poisson noise was considered in both cases. This is not strictly the case for the ROSAT data because of the preprocessing of the raw counts (exposure correction, modelled background subtraction), but this approximation appears to be good enough for detecting correctly all the large-scale components we are interested in. Significant wavelet coefficients were thus selected according to a probability of false detection of 4$`\sigma `$, 3.5$`\sigma `$ and 3$`\sigma `$ for scales 1 (18.5 arcsec), 2 (37 arcsec) and 3 (74 arcsec) and above respectively. The rather high levels for the first two scales automatically eliminate most of the small-scale structures if any (these might distort the overall shape of the brightness profile or might have arisen from an error in the probability density function of the noise or from a loss in control over the bias in the Anscombe transform in regions with very few counts). We have checked in several cases that the profiles do not change significantly when the thresholding for scales 1, 2 and 3 is modified.

### 2.5 Radial profiles

After wavelet reconstructing both EUVE and ROSAT images for each cluster, we proceed to compare the decrease of the surface brightness in the two observational bands. Profiles were thus derived from the 2D reconstructed images within concentric elliptical rings using the STSDAS.ANALYSIS. ISOPHOTE.ELLIPSE task in IRAF. The cluster ellipticities and major axis position angles were estimated from the X-ray images, which have better signal to noise ratio, except for Coma because it fills the 20 arcmin PSPC ring and its ellipticity had to be derived from the EUVE reconstructed image. The ellipticities and major axis position angles were then fixed to be the same in the EUV and X-rays, in order to measure intensities in the same spatial regions.

In order to see how the instrument point spread functions (PSF) could influence our results, we also analyzed a point source in a similar way: the star $`\gamma `$Tau, for which both EUV and X-ray data were available. The EUV and X-ray exposure times for this star were 89924s and 8547s respectively, and the ROSAT PSPC image of $`\gamma `$Tau was divided by its corresponding exposure map, as for clusters. The EUV and X-ray images were rebinned to a pizel size of 18.5$`\times `$18.5 arcsec2 as previously, and the EUV and X-ray profiles were then derived in a similar way. Note that this point source is about 15 arcmin offset from the detector center in the ROSAT image, at a radius where the PSF of PSPC is already degraded. Therefore, the PSF effect in X-rays will in fact always be smaller or equal than what we measure for $`\gamma `$ Tau. This is not the case for EUVE, where the point source is close to the detector center.

The problem of a possible non-flatness of the background was already addressed in the case of EUVE by Lieu et al. (1999c). A comparable wavelet analysis was performed on a blank field totalling 85 ksec exposure time, and the backgroud was found to vary by at most 10% within the central circle of 16 arcmin radius defined above taken into account here (see Fig. 1 by Lieu et al. 1999c). The mean value of another blank field observed for 21 ksec is found to be 1.21 counts/s/pixel2, excluding the outermost regions (that is in a rectangle similar to ABCD from Fig. 1), which is comparable to the mean background value in Abell 1795 (1.19 counts/s/pixel2).

In the case of the X-ray images, they were divided by their corresponding exposure maps, and therefore the non-flatness of the ROSAT PSPC detector should be corrected for. We are therefore confident that the non-flatness of the detector backgrounds cannot change influence our results by more than a few %.

Error bars on all profiles are those directly given by the ELLIPSE task. However, since we were concerned that they might be underestimated, we also computed them from the Poisson fluctuation in the number of counts in each ring ; these error bars were found to be commensurate with those given by ELLIPSE (the rms fluctuation of the pixel intensity within each concentric ellipse). For reconstructed images, they are typically smaller than $`\pm 0.01`$ in logarithmic scale, i.e. too small to be clearly visible on Figs. 4 to 8).

## 3 Results

### 3.1 EUV and X-ray profiles

For co-added images, we do not expect to detect any structure outside a certain radius since the outer part of the image is purely noise. The central circle in Fig. 2 shows a conservative estimate of this radius, i.e. the half width of the EUVE detector $``$16 arcmin, which we will take as the limit of our analysis.

We first derived the profiles of the EUVE and X-ray images of all clusters before and after wavelet analysis and reconstruction, as a test to check that our method conserves the overall profile shape and intensity level. The corresponding profiles are shown in Fig. 4 for Abell 1795. Note that all the profiles drawn in Figs. 4 to 8 are obtained after subtracting the background. For the raw data, this background is estimated as a mean value far from the edges of the detector; for wavelet reconstructed images, the background is modelled during the wavelet analysis itself and subtracted. There is good consistency between the two EUVE curves on the one side and between the two X-ray curves on the other side, implying that the wavelet analysis and reconstruction applied here does not affect fluxes, but simply provides a means of detecting cluster emission out to a somewhat larger radius (as in the case of the EUVE data). The small feature apparent in the X-ray profile at $``$6 arcmin is due to a Seyfert galaxy superimposed on the cluster; it was removed from the wavelet reconstructed image and therefore does not appear on the corresponding profile. Note however that Abell 1795 is the only cluster for which the raw and wavelet reconstructed X-ray profiles do not coincide exactly (in the other clusters they can hardly be distinguished). Although we have tried various thresholding values for the small scales and various background subtractions of the raw data to see if they could coincide better, we have not been able to superimpose them exactly.

The EUV and ROSAT PSPC PSFs drawn from the images of a point source (see Sect. 2.5.) for comparison are also plotted in Fig. 4. They indicate that the PSFs have no influence on the profile shapes for radii larger than about 1arcmin and 2.5 arcmin in the EUV and X-rays respectively. Therefore any features found in the profiles within regions smaller than these radii are probably not real.

The EUV and X-ray radii of Abell 1795 are at least 16 arcmin, corresponding to a physical extent of about 1700 kpc. This extent is notably larger than the previously found value of 10 arcmin (Mittaz et al. 1998, Bonamente et al. 2001b).

The EUV and X-ray radii of Abell 2199 extend to at least 16 arcmin, or 820 kpc (Fig. 5), comparable to the value of 15 arcmin measured by Kaastra et al. (1999) both with EUVE and Beppo-Sax.

In Abell 4059 (Fig. 6), the EUV and X-ray emissions also reach 16 arcmin or 1300 kpc, much further than the 4-5 arcmin reported by Berghöfer et al. (2000b). The slopes of the profiles in the EUV and X-rays differ notably, at least in the inner 8 or 10 arcmin. The difference in the central regions between the EUV profiles derived from the raw and wavelet reconstructed images may be due to the fact that structures only present at small scales were eliminated, and therefore signal may be missing in the very center. This cannot be due to the shape of the EUVE PSF, since it is flat for radii larger than about 1 arcmin while this discrepancy is observed up to a radius of about 3 arcmin. Such a behaviour is not observed in Abell 1795, Abell 2199 or Coma, probably because they are brighter.

The EUV and X-ray radii in Coma reach at least 16 arcmin (Fig. 7), corresponding to 630 kpc and slightly larger than reported by Bowyer et al. (1999). The two profiles appear quite similar at least up to a radius of about 10 arcmin.

Virgo is the most nearby cluster in our sample. EUV and X-ray emissions reach at least 16 arcmin (Fig. 8), that is slightly larger than detected by Lieu et al. (1996a) and Berghöfer et al. (2000a), but smaller than the 20 arcmin claimed by Bonamente et al. (2001b); however, this extent corresponds to the very small physical value of $`80`$ kpc, due to the proximity of Virgo. Both profiles appear quite parallel except in the inner 4 arcmin; note however that since the X-ray PSF is not flat for radii smaller than 2.5 arcmin, the difference between the EUV and X-ray profiles in the inner 4 arcmin may not be real.

In order to quantify the probability for the EUV and X-ray profiles to be similar, we performed a Kolmogorov-Smirnov (hereafter K-S) test on the EUV and X-ray profiles drawn from the wavelet reconstructed images, both in concentric elliptical annuli, and in ellipses with increasing radius, after normalizing the EUV and X-ray profiles to the same innermost pixel value.

In all clusters, the probability that the two profiles are issued from the same parent population varies with radius: it is a few % in the very center (1.5 arcmin radius), and less than 2% in concentric annuli of width 1.5 arcmin (i.e. for radii between 1.5 and 3.0, 3.0 and 4.5 arcmin etc.). If we consider ellipses of radii 3.0, 4.5, 6.0 arcmin etc. instead, the probability that the two profiles are issued from the same parent population is smaller than 0.1% in all cases. Therefore even if the EUV and X-ray intensity profiles may appear similar to the eye, as for example in Coma, it is clear that statistically they are not.

### 3.2 EUV to X-ray intensity ratios

In order to compare the EUV to X-ray profiles, we computed the EUV to X-ray intensity ratios in the same concentric ellipses; these are shown in Figs. 9 to 13. The ratio of the EUV to X-ray PSFs was also derived (see Fig. 9). As expected from the individual EUV and X-ray PSFs, this ratio is only constant for radii larger than 2.5 arcmin, while below it shows a strong dip. Therefore the EUV to X-ray ratios below this radius are not reliable. Such ratios have to be compared to those expected in the hypothesis that the EUV emission is the low energy tail of the hot ICM emission. We therefore employed spatially resolved hot ICM models for each cluster, in order to account for the effects of temperature and abundance gradients on the predicted ratios. Spectra of all cluster regions were modelled with a thin-plasma emission code (MEKAL in XSPEC, Mewe et al. 1985, 1986, Kaastra 1992) modified by photoelectric absorption (WABS in XSPEC, Morrison & McCammon 1983). In detail, the hot ICM in Abell 1795 was modelled with kT increasing from 3.5 keV at the center to 6 keV at radii $``$ 3 arcmin, and abundances in the 0.3-0.5 solar range (Briel & Henry 1996); Abell 2199 was modelled with kT in the 2.5-4 keV range, and abundances decreasing from $``$ 0.7 solar at the center to 0.3 solar for radii $``$ 3 arcmin (Siddiqui et al. 1998); Abell 4059 was modelled with kT increasing from 2 keV at the center to 5 keV at radii $``$ 3 arcmin, and abundances in the 0.3-0.6 solar range (Hwang & Sarazin 1998); Coma was considered as isothermal with kT=8.2 keV and abundances 0.21 solar (Hughes et al. 1993); finally, the hot ICM parameters of Virgo were directly fit to our PSPC data, with kT in the 1.2-3 keV range, and metal abundances of 0.4-0.6 solar. Using these models, we were able to derive the expected count rates in the EUV DS band and in the PSPC 1/4 keV band, taking into account the effective area of the instruments and the performances at all detector positions. The points corresponding to these calculated count ratios are shown in the following figures.

In Abell 1795, the EUV to X-ray intensity ratio (Fig. 9) is smaller than 1 in the center (for radii smaller than 1.5 arcmin), but this is most probably due to the shape of the ratio of the EUV to X-ray PSFs. It then increases with radius to reach a roughly constant value of 4. The shape of this curve is consistent with that given by Bonamente et al. (2001b), while the value of 4 that we find in the outer regions is intermediate between that of 6 given by Mittaz et al. (1998) and that of 2 given by Bonamente et al. (2001b). Error bars are typically smaller than $`\pm 0.01`$ in logarithmic scale, and are therefore difficult to see on the graphs. We therefore confirm the strong soft excess found in this cluster by these two groups of authors, and disagree with Bowyer et al. (1999) who claimed that no excess EUV emission was seen in this cluster. Note that the small disagreement between the raw and wavelet reconstructed X-ray image profiles (Fig. 4) cannot be at the origin of the EUV excess that we observe, since on the contrary the raw X-ray flux is higher that the reconstructed one, and therefore using the raw X-ray data would lead to an even larger EUV excess.

The EUV to X-ray intensity ratio in Abell 2199 is smaller than 1 for radii below 5 arcmin, but as for Abell 1795 this is most probably due to the shape of the ratio of the EUV to X-ray PSFs. It then increases and reaches a value of about 3.7 between radii of 10 and 16 arcmin. We therefore agree with the soft excess observed in this cluster from various sets of data by Kaastra et al. (1999), and here also disagree with Bowyer et al. (1999) who found no excess EUV emission in this cluster.

The EUV to X-ray intensity ratio in Abell 4059 is smaller than 1 for radii below 2 arcmin (Fig. 11), in agreement with Berghöfer et al. (2000b), but again this is most probably due to the shape of the ratio of the EUV to X-ray PSFs; for increasing radii, it strongly increases and reaches a value of almost 20 between radii of 12 and 16 arcmin, in total disagreement with Berghöfer et al. (2000b), who find no EUV emission beyond 4 arcmin and no EUV excess beyond 2 arcmin.

Contrary to what happens in the three previously studied clusters, the EUV to X-ray intensity ratio in Coma appears rather flat up to a radius of at least 10 arcmin (Fig. 12), but as shown above the two distributions are not similar. The fact that no dip is seen within the central few arcminutes, as would be expected from the shape of the ratio of the EUV to X-ray PSFs could indicate that the EUV excess in the center is even stronger than seen in Fig. 12. There is an EUV excess by a factor of $`1.4`$ over a thermal bresmsstrahlung tail.

The EUV to X-ray emission ratio in Virgo is in agreement with the EUV excesses previously reported by Berghöfer et al. (2000a) and Bonamente et al. (2001b), therefore implying that here also the EUV emission is not a thermal bremsstrahlung tail. The rise and fall of the EUV to X-ray ratio for radii increasing between 0.3 and 2.5 arcmin (Fig. 13) is rather surprising since it varies in an opposite way to the ratio of the EUV to X-ray PSFs; however images with higher spatial resolution are obviously necessary to push the analysis further at these radii.

The superimposed EUV to X-ray intensity ratios for all the clusters of our sample are plotted in Figs. 14 and 15 for radii expressed in arcminutes and in kpc respectively. Although the shapes of these curves in the central 2.5 arcmin most probably reflect that of the ratio of the EUV to X-ray PSFs (Fig. 14), we can see that at least beyond this radius all the clusters do not have the same behaviour. This is particularly true when radii are expressed in kpc, and certainly must be linked to a difference in the emission mechanisms of the EUV emission in these clusters. Note that the radial extent of the EUV and X-ray emission in Virgo is comparable to the other clusters when the radius is expressed in arcminutes, but becomes extremely small in physical distance units, due to the very small redshift of the cluster: slightly more than 80 kpc.

## 4 Discussion and conclusions

In this paper, we have shown unambigously the existence of a EUV excess in all five clusters of our sample.

In the first three clusters (Abell 1795, Abell 2199 and Abell 4059), the EUV to X-ray intensity ratios have comparable behaviours: they show a possible deficiency of EUV emission over a bremsstrahlung tail in the very central regions (but this may be an effect of the ratio of the EUV to X-ray PSFs), and a EUV excess beyond a few arcmin. As suggested by several authors (Bowyer et al. 1999; Lieu et al. 2000), the first of these features, if real, can be interpreted as due to excess absorption within the cluster core due to the fact that in the cooler central regions some metals are not fully ionized and these ions absorb part of the soft X-ray flux. The presence of such excess absorption would agree with the previously claimed existence of cooling flows in these three clusters. However, XMM-Newton has detected much weaker emission lines than expected from bona fide isobaric cooling flow models, implying that there is significantly less “cool” gas than predicted (e.g. Kaastra et al. 2001). Therefore it is no longer straightforward to interpret excess absorption at the center of these clusters as due to the presence of cooler gas.

One may however note that the only cluster (Coma) for which no EUV dip is seen in the very center is both the hottest one by far and the only one with no “cooling flow” whatsoever. Since the ratio of the EUV to X-ray PSFs shows a dip for radii smaller than about 2.5 arcmin, the absence of such a dip in Coma suggests that in fact there is a significant EUV excess in the central regions of Coma.

Virgo shows a different behaviour, since the EUV excess is the strongest in the very center, within radii smaller than about 4 arcmin (corresponding to a physical size of $``$20-30 kpc, that is roughly the size of M87), then remains roughly constant with radius between 4 and 7 arcmin, and finally increases with a shallow slope for radii larger than 7-8 arcmin. As for Coma, the presence of a bump in the central regions suggests that the EUV excess is in fact very stronge in this zone. For these last two clusters, images with higher spatial resolution are obviously necessary to analyze the EUV excess at these small radii.

The mere presence of an EUV excess in all five clusters of our sample indicates that a mechanism other than bremsstrahlung is needed to account for this emission. Our data does not allow us to discriminate between the various mechanisms proposed in the literature for the soft excess. However, in view of the most recent results obtained with XMM-Newton suggesting that there is much less warm gas in the central regions of clusters than previously believed, it seems likely that this EUV excess is probably of non thermal origin.

###### Acknowledgements.

We are very grateful to M. Traina, from the Observatoire de la Côte d’Azur, for her invaluable help in obtaining the up-to-date version of the multiscale vision model package used throughout the technical part of this paper. We acknowledge discussions with D. Gerbal, G. Lima Neto, G. Mamon and J. Mittaz. We also warmly thank the referee, Jelle Kaastra, for several pertinent suggestions which helped us to improve the paper. Finally, R. Lieu thanks the Institut d’Astrophysique de Paris, CNRS, and the Université Pierre et Marie Curie for their hospitality.