text:
# Self-organized criticality in the Kardar-Parisi-Zhang-equation

## Abstract

Kardar-Parisi-Zhang interface depinning with quenched noise is studied in an ensemble that leads to self-organized criticality in the quenched Edwards-Wilkinson (QEW) universality class and related sandpile models. An interface is pinned at the boundaries, and a slowly increasing external drive is added to compensate for the pinning. The ensuing interface behavior describes the integrated toppling activity history of a QKPZ cellular automaton. The avalanche picture consists of several phases depending on the relative importance of the terms in the interface equation. The SOC state is more complicated than in the QEW case and it is not related to the properties of the bulk depinning transition.

PACS: 05.70.Ln, 64.60-i, 45.70.Ht

One of the more interesting developments in non-equilibrium statistical physics has been the realization that inhomogeneous systems give rise to new, rich behavior compared to those in which translational invariance is present. A prime example thereof is the criticality found in so-called sandpile models, coined ’self-organized’ due to the absence of an obvious control parameter .

Usual ’self-organized critical’ SOC sandpiles have two central features. First, the system is driven very slowly so that individual avalanches occur on a separate, fast time-scale. Second, the drive (like the addition of grains) and the dissipation (loss of grains) occur non-uniformly. Typically the system is forced, on the average, in an uniform fashion but open boundaries create dissipation, in sandpiles via the loss of grains. These two characteristics are the keys to describe, and the earmarks of SOC (for early work see ) .

Essential to a wide class of SOC models is that the movement of grains has a Laplacian character. Consider now the history of activity, topplings in such models. The activation of a site means that it has received enough grains, $`z(x,t)`$, from outside or through redistribution, to overcome a local activation threshold: $`z(x,t)>z_c(x,t)`$. For a large class of models, this implies that if $`H(x,t)`$ defines the time-integrated activity, $`H=^t\rho (x,\tau )𝑑\tau =^t\mathrm{\Theta }(z(x,\tau )z_c(x,\tau ))𝑑\tau `$, then after discretization $`H`$ obeys the discrete interface equation

$`H(x,t+1)`$ $`=`$ $`\{\begin{array}{cc}H(x,t)+1,\hfill & f(x,t)>0,\hfill \\ H(x,t),\hfill & f(x,t)0,\hfill \end{array}`$ (3)
where $`f`$ denotes the local ’force’ acting on the interface . This can be recast in the form

$$\frac{\mathrm{\Delta }H}{\mathrm{\Delta }t}=\mathrm{\Theta }\left(\nu ^2H+\eta (x,H)+F(x,t)F_c\right).$$
(4)
The Laplacian comes from the average redistribution of grains, and a constant $`\nu `$ takes into account the specific updating rules of the system. $`\eta `$ is a quenched noise term that originates from the random rules of sandpiles. $`F(x,t)`$ measures the external drive and changes on the slow-time scale, and $`F_c`$ is the toppling threshold written as a force acting on the interface. By casting the dynamics of a sandpile in the form of a driven interface in a random enviroment ($`\eta `$) one sees that the SOC critical point is reached under the combined action of a slowly-increasing external force $`F(x,t)`$ and the boundary conditions, $`H=0`$, which amount to the loss of grains. This is a slightly unusual case of a depinning transition of the quenched Edwards-Wilkinson universality class, since the macroscopic shape of the interface is parabolic (in 1D) . The avalanche picture in an ordinary depinning transition  is however different since the translational invariance is here broken as manifested already in the average interface profile - which tells us that the activity density peaks in the bulk, in the center of the system.

The connection of SOC to non-equilibrium phase transitions can be made both in the interface context and in the realm of absorbing state phase transitions, exemplified by the contact process and various reaction-diffusion models . One way to classify SOC models is to consider the noise correlator in the presentation at hand. Such analysis however assumes certain symmetries included in the representation (the QEW Langevin equation; the Langevin equation describing the active site density). A pertinent question thus arises: what happens if the ensemble, in which SOC is observed, is varied with respect to its fundamental symmetries? In this work we illustrate this by analyzing the outcome for the Kardar-Parisi-Zhang universality class. The difference between the QKPZ and the EW equations (or a non-linear Langevin equation compared to a linear diffusion equation) is that the diffusion equation is insensitive to most initial profiles due to the fundamental symmetries of the equation, while in the QKPZ case the growth component perpendicular to the interface slope is always non-zero. This means that there will be a new history dependence in the QKPZ interface which is absent from the EW dynamics and from usual absorbing state phase transitions.

Consider, in analogy to the sandpile interface representation, the quenched KPZ equation  which reads

$`{\displaystyle \frac{\mathrm{\Delta }H}{\mathrm{\Delta }t}}=\mathrm{\Theta }(\nu ^2H+{\displaystyle \frac{\lambda }{2}}(H)^2+\eta (𝐱,H)+F(x,t)),`$ (5)
where the nonlinearity proportional to $`\lambda `$ has been added to the QEW equation. Note the step function that enforces the constraint that the interface cannot propagate backwards. In the case of the sandpile equation, it also discretizes the interface velocity to zero or unity which in the absence of mixing by randomness creates an important Abelian symmetry and thus becomes crucial in that case . The external drive is in the following taken to be uniform in space, $`F(x,t)F(t)`$ and the noise to be quenched and delta-correlated, $`\eta (𝐱,H)\eta (𝐱^{},H^{})=2D\delta ^d(𝐱𝐱^{})\delta (HH^{})`$. Equation (5) has to be augmented with suitable boundary conditions in order to reach a steady-state. As in the QEW or sandpile ensemble, we use the selection $`H(x,t)=0`$ at $`x=0`$, $`x=L+1`$ (in 1+1 dimensions) instead of e.g. periodic boundary conditions. In the case of such ordinary boundary conditions Eq. (5) has a depinning transition at a critical $`F_c`$. The exponents depend on whether $`\lambda `$ approaches zero with $`\mathrm{\Delta }FFF_c0`$ or not, and in the latter case we have a true QPKZ depinning transition .

An interface evolving from a flat boundary condition will develop a non-zero roughness and finally stops due to the influence of boundary conditions. The roughening, the development of average local slopes, depends on the sign of $`\lambda `$. To reach SOC-like conditions, the separation of timescales, or a vanishingly small but positive interface velocity (average activity level) $`v=0^+`$ one has to choose the right way of driving the system, $`F(t)`$. For $`F(t)>0`$ the interface becomes unstable if $`\lambda >0`$ and $`\frac{dv}{dt}>0`$. This breaking of up/down symmetry, absent in the continuum QEW model, shows the first sign that the possible SOC state in the QKPZ equation is different. If $`\lambda <0`$ it becomes posssible to take the SOC limit of vanishingly slow external drive,

$$F(t)=\zeta t,\zeta 0$$
(6)
from which follows that a separation of time-scales happens, similarly to the QEW or normal sandpile case. In the following we investigate the ensuing interface dynamics.

To this end, a cellular automaton  version of Eq. (5) is defined as follows. Let $`f_i`$ be a microscopic force at the discretized sites $`xi`$, such that for $`f_i>0`$ $`i`$ topples. $`f_i`$ reads in 1D

$`f_i(t)=\nu \left[H_{i+1}(t)+H_{i1}(t)2H_i(t)\right]`$ (7)
$`+{\displaystyle \frac{\lambda }{4}}\left([H_i(t)H_{i1}(t)]^2+[H_{i+1}(t)H_i(t)]^2\right)`$ (8)
$`+F(t)+\eta _{i,H_i}.`$ (9)
This can be interpreted, in analogy to SOC ones, as a QKPZ sandpile where a site gives $`2\nu `$ particles or units of energy to its neighbors once it topples, and the toppling criterion also compares the integrated activity at site $`i`$ to its neighbors. In the interface picture $`H_i(t+1)=H_i(t)+\mathrm{\Theta }(f_i(t))`$. With this definition, $`H(x,t)=^t\rho (x,\tau )𝑑\tau `$ maps the discrete QKPZ equation to the ’sandpile’ and vice versa. Notice that the nonlinearity, proportional to $`\lambda `$, above is discretized in a particular way, to avoid the numerical instabilities of the QKPZ equation in a reasonable fashion .

Fig. 1 shows a series of snapshots from a simulation of the QKPZ equation with $`\eta _{i,H_i}=\pm g`$ with equal probability. This describes random toppling thresholds that change after each toppling in the sandpile language. The system is driven in such a way that $`F(t)`$ is increased after an avalanche has relaxed (defined in the usual way by the stoppage of all activity), so that one site becomes active, and the interface advances. It is evident that the dynamics of such a system consists of at least of the development of an initial profile, and the subsequent phase(s). The macroscopic shape of the QKPZ interface is to first order triangular in 1+1-D as expected with $`\lambda <0`$. The slope follows from the different terms in Eq. (5). For a flat slope the surface tension is negligible, and the slope arises from the balance of the nonlinearity and the driving force. This gives $`|h|=\sqrt{2\zeta t/\lambda }`$. The development of such a profile, determined by the symmetries of the interface equation, is analogous to the development of the steady-state grain configuration in ordinary SOC sandpiles from an empty system. Next we discuss the phase that appears after the initial transient.

As in ordinary SOC, avalanches are characterized by the duration $`T`$, area $`l`$, and number of topplings $`s`$ (the number of times sites advance). The initial growth, giving rise to the triangular shape, is neglected. We obtain for the avalanche sizes in this regime the result shown in Fig. 2 for four system sizes. The size distribution has *asymptotic power law decay* with a measured exponent of $`\tau _s=2.5\pm 0.1`$ with $`P(s)s^{\tau _s}`$. Due to the large $`\tau _s`$ exponent the average avalanche size is not dependent on $`L`$, unlike in usual SOC sandpile models.

However, the avalanches are time-translation invariant only in the average sense, due to the role of the average interface slope $`m`$. The avalanche sizes are plotted versus time on Fig. 3. On the same diagram, the standard deviation $`D`$ of the local slope is also shown.

The largest avalanches take place when the roughness of the interface $`D`$ is the largest. Smaller, but still relatively large avalanches develop where the sides of the triangular interface are smooth. The local slopes, i. e. $`|H_iH_{i1}|`$ alternate between $`n`$ and $`n+1`$, where $`n`$ is a positive integer. It is not surprising that a smooth interface evokes large avalanches, because there are no fluctuations along the surface that would stop an upwards propagation of the avalanche. When $`D`$ is the smallest, the first few avalanches start from either bottom of the triangle and propagate up to the top, increasing the interface height by one at every site. The really large avalanches occur at the maxima of $`D`$. They typically do not increase the average gradient of the side, but only the height of interface evenly along the side. In spite of the extra slope-dependence, there is no characteristic size as is typical of SOC and the correlation length equals the system size. The period of the avalanches is simply related to the squared system size. This is since the area covered by all the avalanches between two integer values of $`m`$ is proportional to $`L^2`$, and, the average avalanche size and duration are of the order of unity: $`s𝒪(1)`$ and $`T𝒪(1)`$.

The support or the area of the avalanche is as usual defined by the number of sites that have toppled at least once during an avalanche. Notice that an avalanche always tends to spread upwards, and thus in 1D the starting and final, extremal sites usually give the area directly. Sometimes, the Laplacian makes it so that an avalanche can propagate over the top of the ’triangle’. We find that the probability of starting an avalanche linearly increases with the location of a site and the probability of an avalanche to stop grows linearly also as the top is approached.

The area of the avalanches, defined as $`l=i_ri_l+1`$, is expected to follow a similar power law as the avalanche sizes do and one obtains that $`\tau _l3.0\pm 0.1`$. The avalanche lifetime $`T`$ is the time that elapses between two successive stationary states in the system and also behaves as a power-law, $`P(T)T^{2.9\pm 0.2}`$. Due to the fact that the average avalanche properties are defined by the small avalanches no analogy of the usual scaling relations between the exponents seems to exist.

The above results have, however, been obtained in just a window of time. Consider an interface (side of the triangle) with an uniform slope $`m`$ and an avalanche which has just started to propagate and reached site $`i`$. The local force at site $`i`$, at time $`t`$, becomes $`f_i(t)=\nu (m1m)+\frac{\lambda }{4}\left[(m1)^2+m^2\right]+F+\eta _i2=\nu +\frac{|\lambda |}{4}(2m1)\frac{|\lambda |}{4}2m^2+F+\eta _i`$ and at $`i+1`$ $`f_{i+1}(t)=\nu [m(m1)]+\frac{\lambda }{4}\left[m^2+(m1)^2\right]+F+\eta _{i+1}=\nu +\frac{|\lambda |}{4}(2m1)\frac{|\lambda |}{4}2m^2+F+\eta _{i+1}.`$ Assuming that on average on an even interface $`F`$ balances the local forces for slope $`m`$ ($`F\frac{|\lambda |}{4}2m^2`$) and neglecting the noise fluctuations, we find a critical slope

$$m_{cr}\frac{2\nu }{|\lambda |}+\frac{1}{2}.$$
(10)
This means that the cell at $`i`$ may stay active for the next updating step at time $`t+1`$ when $`m>m_{cr}`$. The natural progress of avalanches is now such that the site at $`i+1`$ is normally always activated if the interface is smooth enough. In this phase the site at $`i+1`$ gets active and topples at $`t+1`$, but at the same time $`H_iH_i+1`$ as well. Thus, a short region of the interface with a smaller gradient will be formed and starts to propagate upwards, giving rise to very large avalanches. Since these are more difficult to induce, the majority of the avalanches are of size unity, and tend to smoothen out the interface fluctuations. The avalanche statistics data presented previously is restricted to systems with a maximum slope of 10: the critical slope for those systems is $`m_{cr}220/4+1/2=10.5`$. In this regime above the critical slope $`m_{cr}`$, the scaling structure of avalanche behavior is lost, and the growth process produces either extremely large or very small avalanches. This feature makes it prohibitively difficult to take this regime under closer scrutiny with numerical simulations. The time to reach $`m_{cr}`$ can be estimated in two ways depending on whether one uses macroscopic time or ’SOC’ time (counting the number of avalanches). With a uniformly, but very slowly increasing force $`F`$ the critical slope is reached at the same time regardless of the system size $`L`$. In contrast, since the average size of an avalanche is independent of $`L`$ (recall that $`\tau _s>2`$), it follows that the number of avalanches $`n_{cr}`$ to reach $`m_{cr}`$ will depend on $`L`$ as $`n_{cr}L^2`$.

The previous results illustrate the nature of SOC in the QKPZ universality class. The combination of slow drive and boundary pinning chooses a certain symmetry for the average interface profile. This can be further elaborated by considering the continuum equation, and writing the interface field on one of the sides of the triangle as $`H(x,t)=m(t)x+\delta H(x,t)`$. Inserting this into the QKPZ equation we obtain for $`\delta H`$

$`{\displaystyle \frac{\delta H(𝐱,t)}{t}}=\nu ^2\delta H+\lambda m(t)\delta H+{\displaystyle \frac{\lambda }{2}}(\delta H)^2`$ (11)
$`+\eta (𝐱,H),`$ (12)
where the mean-field solution for $`H`$ has been subtracted from the equation and the result is valid for intermediate times $`\delta tt`$. We discover that the externally imposed slope (by the combination of $`F`$ and $`H=0`$ at the boundaries) changes the effective interface equation by creating a term linear in $`\delta H`$, that will dominate over the KPZ-nonlinearity. The situation is thus exactly analogous to depinning transitions in the presence of an anisotropy, discussed by Kardar, Dhar, and Tang . This explains qualitatively why the avalanche behavior at criticality, at the SOC depinning transition, is different from that of the bulk depinning transition for the original equation. Regardless of whether the SOC ensemble differs from the translationally invariant one for ordinary SOC sandpiles, it is thus clear that in the case of SOC obtained with a similar recipe in QKPZ-systems the two ensembles do not share the same exponents. Should it be possible to write an activity-based description for the QKPZ universality class, this implies that the field theory should have a different structure in the two cases, in contrast to the connection between the QEW and absorbing state phase transitions with a conserved field .

Above, we have concentrated on a cellular automaton for the QKPZ equation since as discussed it is related to a microscopic ’activity’ picture, or a “QKPZ sandpile”. The slow-drive SOC-depinning limit can be taken also in the continuum version. Our investigations indicate that the same phenomenology of symmetry breaking persists in that case as well, though the actual exponents in the second phase before the critical slope is reached differ as $`\mathrm{\Delta }h`$ can now have arbitrary values and the avalanche properties change. The analysis of the continuum model is hampered by well-known numerical instabilities that develop once the smoothing effect of the Laplacian is reduced - which is inevitable here due to the increasing average slopes. We have also studied the two-dimensional case. While there, with periodic boundary conditions in one direction, and SOC ones in the other, are interesting details concerning the geometry of avalanche spreading, the basic physics is the same as in 1D.

In conclusion, in this paper we have analyzed the behavior of interfaces in the QKPZ universality class in the presence of conditions that produce self-organized criticality. Spatial and temporal invariance are broken by the slow drive and the dissipative pinning boundary conditions, giving rise to a macroscopic interface shape. The shape, in combination with the nonlinearity, results in the creation of a term linear in $`\delta h`$ in the effective equation that governs the physics of avalanches. In contrast to usual SOC sandpiles we find that there are two avalanche regimes depending on the microscopic dynamics, where the second one fails to show critical scaling seen in the first. The ensemble in which the criticality appears is fundamentally different from the usual depinning one and is related to the anisotropic QKPZ class .

Partial support of OTKA T029985 is acknowledged.


format: text

text:
# Synchronization in Small–world Systems

## Abstract

We quantify the dynamical implications of the small-world phenomenon. We consider the generic synchronization of oscillator networks of arbitrary topology, and link the linear stability of the synchronous state to an algebraic condition of the Laplacian of the graph. We show numerically that the addition of random shortcuts produces improved network synchronizability. Further, we use a perturbation analysis to place the synchronization threshold in relation to the boundaries of the small-world region. Our results also show that small-worlds synchronize as efficiently as random graphs and hypercubes, and more so than standard constructive graphs.

Recently, Watts and Strogatz Watts and Strogatz (1998) showed that the addition of a few long-range shortcuts to an otherwise locally connected lattice (the ”pristine world”) produces a sharp reduction of the average distance between arbitrary nodes. The ensuing semi-random lattice was denoted a small-world (SW) because the sudden appearance of short paths occurs early on, while the system is still relatively localized. This concept has wide appeal: the SW property seems to be a quantifiable characteristic of many real-world structures Watts and Strogatz (1998); Strogatz (2001); Watts (1999); Jeong et al. (2000), both human generated (social networks, WWW, power grid), or of biological origin (neural and biochemical networks).

A spur of ongoing research Strogatz (2001) has concentrated on static and combinatoric properties Barthelemy and Amaral (1999); Barrat and Weigt (2000); Newman et al. (2000); Newman and Watts (1999); Moore and Newman (2000); Kulkami et al. (2000) of a tractable SW model Watts and Strogatz (1998); Monasson (1999). Monasson Monasson (1999) considered the SW effect on the distribution of eigenvalues of the connectivity matrix (the graph Laplacian) which specifies the coupling between nodes—a relevant topic for polymer networks Jespersen et al. (2000). However, despite their central role in real-world networks, there are fewer studies of dynamical processes taking place on SW lattices. Among those, automata epidemics simulations Kuperman and Abramson (2001) and Web-browsing studies Kleinberg (2000) have revealed the importance of shortcuts. Numerical work on synchronization of Kuramoto oscillators Watts (1999), discrete maps Gade and Hu (2000) and Hodgkin-Huxley neurons Lago-Fernández et al. (2000) has shown improved SW synchronizability, as intuitively expected. However, these numerical examples are not generic, and fail to provide insight into how the SW property influences the dynamics.

In this paper, we explicitly link the SW addition of random shortcuts to the synchronization of networks of coupled dynamical systems. This is an example of dynamics on networks—leaving aside the distinct problem of evolution of networks here. By using a generic synchronization formulation Pecora and Carroll (1998); Fink et al. (2000) to factor out the connectivity of the network, we identify the synchronization threshold with an algebraic condition of the graph Laplacian. Through numerics and analysis, we quantify how the SW scheme improves the synchronizability of the pristine world, mainly as a result of the steep increase of the first-non-zero eigenvalue (FNZE). The synchronization threshold is found to lie in the SW region Kuperman and Abramson (2001); Watts (1999), but does not coincide with its onset—it can in fact be linked to the effective randomization that ends SW. Within this framework, we show that the synchronization efficiency of semi-random SW networks is higher than standard deterministic graphs, and comparable to both fully random and ideal constructive graphs.

Consider $`n`$ identical dynamical systems (placed at the nodes of a graph) that are linearly and symmetrically coupled (as represented by the edges of the undirected graph) with global coupling strength $`\sigma `$. The topology of the graph can be encoded in the Laplacian matrix $`G`$, a symmetric matrix with zero row-sum and real spectrum $`\{\theta _k\},k=0,1,\mathrm{},n1`$. A general linear stability criterion for the synchronized state of the system Pecora and Carroll (1998); Fink et al. (2000) is given by the negativity of the master stability function $`\lambda _{\mathrm{max}}(\sigma \theta _k)<0,k`$. This $`\lambda _{\mathrm{max}}`$ is a characteristic of the particular dynamics at the nodes but, crucially, a large class of oscillatory systems (chaotic, periodic and quasiperiodic) have master stability functions with generic features Fink et al. (2000). In particular, for several chaotic systems $`\lambda _{\mathrm{max}}`$ has a single deep well, as depicted in Fig. 1. (We remark that this analysis is quite general: it can be extended Fink et al. (2000) to eliminate the zero row-sum constraint, and to comprise nonlinear coupling and more general synchronization criteria.) Stability is thus ensured by tuning the coupling $`\sigma `$ to try and place the entire spectrum of transverse eigenvalues (times $`\sigma `$) in the deep, stable region: $`\sigma \theta _k(\alpha _1,\alpha _2)`$. This leads to an algebraic condition for the existence of a linearly stable synchronous state: a network is synchronizable if

$$\theta _{\mathrm{max}}/\theta _1<\alpha _2/\alpha _1\beta ,$$
(1)
where $`\theta _1`$ is the FNZE and $`\theta _{\mathrm{max}}`$ is the maximum eigenvalue of the Laplacian $`G`$. The figure of merit ($`\beta `$) ranges from 5 to 100 for a variety of oscillators (e.g., Lorenz, Rössler, double scroll).

Small-worlds are generated from a pristine world: a $`k`$-cycle of $`n`$ nodes and range $`k`$, each node coupled to its $`2k`$ nearest neighbors for a total of $`nk`$ edges Monasson (1999). The Laplacian of this graph $`G^0`$ is a banded circulant matrix with non-zero elements on the main diagonal and the $`2k`$ adjacent diagonals: $`G_{ii}^0=2k`$ and $`G_{ij}^0=G_{ji}^0=1`$ with $`(i+1)modnj(i+k)modn`$, and $`1in`$. The SW scheme dopes the pristine world by adding $`ns`$ edges picked at random from the $`n(n2k1)/2`$ remaining pairs. Each new edge between nodes $`l`$ and $`m`$ adds off-diagonal $`\mathrm{\Delta }G_{lm}=\mathrm{\Delta }G_{ml}=1`$ and on-diagonal $`\mathrm{\Delta }G_{ll}=\mathrm{\Delta }G_{mm}=1`$ contributions to the Laplacian, thus preserving the null row-sum and the bidirectional coupling. The average number of shortcuts per node ($`s`$) is related to other measures of randomness ($`p`$ and $`q`$) used previously Watts and Strogatz (1998); Monasson (1999): $`skpq(n2k1)/(2n)`$.

The numerical results in Fig. 2 illustrate the SW effect on the synchronization of networks of different size and range. For concreteness, all our numerics have been performed for a network of identical $`x`$-coupled Rössler chaotic oscillators with $`\beta 37.85`$. Similarly to other locally connected networks, pristine worlds have a large eigenratio $`\theta _{\mathrm{max}}/\theta _1`$ (i.e., they are difficult to synchronize). However, as $`s`$ is increased the eigenratio falls sharply until, at a value $`s_{\mathrm{sync}}`$, the condition (1) is reached (i.e., the addition of shortcuts makes it synchronizable). The dependence of $`s_{\mathrm{sync}}`$ on the network parameters $`\{n,k\}`$ is notably complicated. First, there appears to be an optimal range $`k4`$ for which the SW is most efficient. Moreover, the synchronization threshold $`s_{\mathrm{sync}}`$ lies in the small-world region but does not seem to coincide with its onset. The SW onset ($`s_L`$) is defined Watts and Strogatz (1998); Watts (1999) by the decay of the average graph distance Newman et al. (2000) $`L(s)(n/k)f(ns)`$, where $`f(x)=\mathrm{arg}\mathrm{tanh}\left(x/\sqrt{x^2+2x}\right)/\sqrt{4(x^2+2x)}`$. Fixing $`L(s_L)/L(0)=3/4`$, we obtain $`s_L1.061/n,nk1`$. The end of the SW region ($`s_C`$) corresponds to the effective graph randomization through the loss of transitivity Watts and Strogatz (1998); Newman et al. (2001); Watts (1999), as given by the decay of the clustering coefficient $`C`$ Barrat and Weigt (2000); Newman : $`C(s)/C(0)(2k1)/(2k(1+s/k)^21),n1`$. Again, fixing $`C(s_C)/C(0)=3/4`$, we obtain $`s_Ck\left(1+[(8k1)/(6k3)]^{1/2}\right)0.155k`$. The synchronization threshold generally lies between these two boundaries which scale differently with $`n`$ and $`k`$: $`s_L1.061/n<s_{\mathrm{sync}}<s_C0.155k`$.

We can gain insight into the synchronization threshold through an analytical perturbation of the eigenratio of the SW Laplacian $`G=G^0+G^r`$. Here, $`G^0`$ is the deterministic Laplacian of the pristine world, and $`G^r`$ is the stochastic Laplacian for the random shortcuts: $`G_{ij}^r=G_{ji}^r=\xi _{ij}`$ (for $`i+k+1j\mathrm{min}\{n,nk+i1\}`$ with $`1ink+1`$); $`G_{ij}^r=0`$ (otherwise); and $`G_{ii}^r=_{j=1}^nG_{ij}^r`$ (for $`1in`$). The $`\xi _{ij}`$ are $`n(n2k1)/2`$ i.i.d. Bernoulli random variables which take the value 1 with probability $`q/n2s/(n2k1)`$ (and the value 0 with probability $`1q/n`$). The circulant $`G^0`$ is Fourier-diagonalizable Monasson (1999) with spectrum $`\theta _j^0=(2k+1)\mathrm{sin}[(2k+1)\pi j/n]/\mathrm{sin}[\pi j/n],\mathrm{\hspace{0.17em}\hspace{0.17em}1}jn1`$, (plus $`\theta _0^0=0`$ of any Laplacian). The FNZE and the maximum eigenvalue of the unperturbed lattice are:

$`\theta _1^0`$ $``$ $`2\pi ^2k(k+1)(2k+1)/(3n^2),kn`$ (2)
$`\theta _{\mathrm{max}}^0`$ $``$ $`(2k+1)+\mathrm{csc}\left[{\displaystyle \frac{3\pi /2}{2k+1}}\right]`$ (3)
$``$ $`(2k+1)[1+2/3\pi ],k1,`$ (4)
where (3) follows from a continuum approximation.

Following an “honest” treatment Boyce (1968) with $`G^r`$ as the perturbation, we treat the analytical expressions of the doubly degenerate eigenvalues as random variables to obtain their expectations. We postpone the detailed calculations Barahona and Pecora  and sketch here the main results. After some stochastic calculus, the expectations of the eigenvalues of the SW Laplacian to second order are shown to be:

$`\theta _i^{(1)}`$ $``$ $`q\pm \sqrt{3\pi q/4n}`$ (5)
$`\theta _i^{(2)}`$ $``$ $`{\displaystyle \frac{2q}{n}}{\displaystyle \underset{m=1}{\overset{n}{}}^{}}(\theta _i^0\theta _m^0)^1,`$ (6)
for $`q/n`$ and $`k/n`$ small. To improve the accuracy of $`s_{\mathrm{sync}}`$, we have also obtained an approximation to (6) for FNZE:

$`\theta _1^{(2)}{\displaystyle \frac{2q}{K^3}}\left[{\displaystyle \frac{9n}{\pi ^2}}+K^2\left({\displaystyle \frac{7}{5}}+{\displaystyle \frac{6}{\pi ^2}}\right)K{\displaystyle \frac{2}{\pi }}\right],`$ (7)
where $`K=2k+1`$. Eqns. (1),(2),(4),(5), and (7) are then used to obtain an estimate of $`s_{\mathrm{sync}}`$ as the solution of an algebraic equation involving only $`n`$ and $`k`$:

$$\theta _{\mathrm{max}}^0+\theta _{\mathrm{max}}^{(1)}=\beta \left(\theta _1^0+\theta _1^{(1)}+\theta _1^{(2)}\right).$$
(8)
As shown in Fig. 2, this approximates well our numerics for $`n^{1/3}<kn`$, where the Rayleigh-Schrödinger perturbation expansion is valid.

Using (8), we can obtain Barahona and Pecora  a first order estimate for the maximum of the synchronization threshold $`s_{\mathrm{sync}}^{}`$ in the valid range. The maximum occurs at $`k^{}n\sqrt{(1+2/3\pi )/2\pi ^2\beta }`$ with the asymptotic value $`s_{\mathrm{sync}}^{}(2+4/3\pi )(12k^{}/n)/3(\beta 1)`$. Therefore, $`s_{\mathrm{sync}}<s_C[2(1+2/3\pi )/(2\sqrt{3}3)(\beta 1)]`$ is bounded by the end of the SW region but linked to it. For $`k<n^{1/3}`$, the eigenvalue bunching (quasi-degeneracy) in the pristine lattice renders the doubly degenerate perturbation invalid. We are developing another approximation to quantify the behavior in this limit, but our numerics Barahona and Pecora  indicate that the dependence of $`s_{\mathrm{sync}}`$ with $`n`$ is sub-logarithmic. This confirms that the synchronizability is most effectively improved for small-range networks (Fig. 2).

How efficient is the addition of random shortcuts for synchronization? We have compared the semi-random SW approach with purely random and purely deterministic schemes. An example of the latter is the synchronization of pristine worlds through the increase of the range $`k`$. From (2) and (4), the eigenratio of a pristine lattice is

$$\frac{\theta _{\mathrm{max}}^0}{\theta _1^0}\frac{3\pi +2}{2\pi ^3}\frac{n^2}{k(k+1)}.$$
(9)
Therefore, $`n`$ nodes can be synchronized in a $`k`$-cycle if $`k>k_{\mathrm{min}}n\sqrt{(3\pi +2)/2\pi ^3\beta }`$. (Note in Fig. 2 the consistency of our analytical approximation (8): $`s_{\mathrm{sync}}=0`$ precisely at $`k_{\mathrm{min}}`$.) For purely random graphs $`𝒢_{n,f}`$ Mohar (1997),

$$\frac{\theta _{\mathrm{max}}^{RG}}{\theta _1^{RG}}\frac{nf\sqrt{2f(1f)n\mathrm{ln}n}}{nf+\sqrt{2f(1f)n\mathrm{ln}n}},$$
(10)
where $`f`$ is the number of edges measured as a fraction of the complete graph. These are compared with SW graphs in Fig. 3.

We remark on several observations regarding Fig. 3. First, the SW addition of shortcuts is more efficient than the deterministic addition of short-range layers. Second, the effective randomization of SW lattices with edge addition translates into converging synchronization behaviors of random and SW networks at large $`f`$. The $`f1`$ region is thus robustly stable: cutting connections from the fully connected graph has very little effect on synchronization stability—not until over 90% are cut (for $`k`$ small) does the eigenratio begin to change drastically. Moreover, if we interpret the number of edges needed to synchronize $`n`$ nodes as a simple measure of “cost”, adding many connections buys little extra stability beyond the small-world regime. Finally, it can be shown Barahona and Pecora  that the general trends of the eigenratio in Fig. 3 (namely, “hyperbolic” dependence for $`f`$ small, and near independence for $`f`$ large) can be predicted with the naive perturbation result that both $`\theta _1`$ and $`\theta _{\mathrm{max}}`$ change linearly with the number of added connections.

We have also compared the synchronization cost (in edges) for SW systems and regular (constructive) lattices (Fig. 4). For $`x`$-Rössler systems in a $`k`$-cycle, this cost tends to $`f=2k_{\mathrm{min}}/(n1)0.140`$. Other constructive lattices Schwenk (1974); Barahona and Pecora  also tend to constant fractions: $`f=0.252`$ (for $`k`$-wheels), $`f=0.070`$ (for $`k`$-Möbius ladders), and $`f=0.053`$ (for the most economical bipartite graph). In all those cases, the cost of synchronization is high: the necessary number of edges scales like $`n^2`$, just like the complete graph. At the other end of deterministic graphs lies the quasi-optimal (though virtually unrealizable) hypercube, which is always synchronizable with a number of edges $`f\mathrm{log}_2n/n`$. This behavior is similar to that of random graphs: from Eq. (10) almost sure synchronization of $`𝒢_{n,f}`$ is asymptotically achieved when $`f\mathrm{ln}n/n`$. Remarkably, Fig. 4 shows that the cost of synchronizing small-$`k`$ SW networks is low, i.e., comparable to that of random graphs and hypercubes Bollobás and Chung (1988).

These results hint at research that could deepen the connections between topology and dynamics on networks. With a view to improved design, the dynamic eigenratio criterion can be related to other graph-theoretical properties (e.g., connectivity, diameter, and convergence of Markov chains) Chung (1997). Moreover, other measures of cost (e.g., robustness under edge deletion) should be considered as possible design constraints. Finally, recent results could lead to extensions of this work to incorporate more general concepts of stability Parrilo (2000), and broader definitions of small-world lattices Newman et al. (2001).

We thank Steve Strogatz for his deep and insightful involvement in this work, and Mark Newman for sharing computer code and unpublished results.


text:
# Quantum oscillations without quantum coherence

## Abstract

We study numerically the damping of quantum oscillations and the increase of entropy with time in model spin systems decohered by a spin bath. In some experimentally relevant cases, the oscillations of considerable amplitude can persist long after the entropy has saturated near its maximum, i.e. when the system has been decohered almost completely. Therefore, the pointer states of the system demonstrate non-trivial dynamics. The oscillations exhibit slow power-law decay, rather than exponential or Gaussian, and may be observable in experiments.

For a quantum system prepared in a linear superposition of its eigenstates, some observables can oscillate with time. Interaction of the system with its environment leads to a decay of the system’s initial pure state into a mixture of several ”pointer states”; it causes an increase of the system’s entropy (decoherence) and damping of quantum oscillations (dephasing) with time decgen ; decgen1 . Both effects, decoherence and dephasing, are often considered as equivalent results of the mixed state of the system. But careful analysis shows important differences decdeph ; zeh1  originating from the fact that the same density matrix can describe both an ensemble of similar systems, and a single decohered system. Thus, e.g., dephasing can appear in an ensemble of pure, non-decohered systems with slightly differing dynamics (this is an idealized picture of $`T_2`$ processes in NMR). Decoherence and dephasing are hard to distinguish in experiments which employ ensembles of quantum systems. However, recently it has become possible to study single quantum systems, such as trapped ions ions , atoms in cavities atoms , or even mesoscopically big Cooper-pair boxes nakamura . As a result, theoretical consideration of the relation between dephasing and decoherence has become experimentally relevant and important.

In this work, we compare dephasing and decoherence in single systems of interacting $`s=1/2`$ spins coupled to a bath of $`s=1/2`$ spins. We show that in some cases, the oscillations can survive long after the entropy has almost saturated, i.e. that the quantum oscillations can take place for a long time even in an almost completely decohered system. These oscillations do not decay according to usual exponential ($`\mathrm{exp}(t/\tau )`$) or Gaussian ($`\mathrm{exp}(t^2/\tau ^2)`$) law, but exhibit long power-law ($`1/\sqrt{t}`$ or $`1/t`$) tails. This result has interesting consequences. The standard picture of a decoherence process assumes that as soon as the system has decayed into a mixture of pointer states, the fast quantum mechanical motion is over, i.e. the pointer states are essentially static. This has been confirmed by numerous studies of different types of pointer states decgen . However, we observe that after the system has decayed into a mixture of the pointer states, and its entropy has reached maximum, the oscillations still persist. It means that the pointer states are not static: they exhibit non-trivial dynamical behavior. We show this explicitly by analyzing the structure of the density matrix. The models considered here may be relevant for a number of experimental systems stamprok .

Studies of the dynamics of quantum oscillations have a long history. In most situations considered so far, quantum oscillations exhibit very fast exponential or Gaussian decay (it is the main reason for the classical behavior of the world around us) decgen . Power-law-damped oscillations have been mentioned stamprok  for a single spin decohered by a spin bath, but the dynamics of decoherence and the structure of pointer states has not been analyzed in detail. In contrast, in this study the properties of the pointer states are of primary importance, especially for a central system containing several spins, where the non-trivial dynamics of the pointer states is even more pronounced. Generally, to our knowledge, the possibility of quantum oscillations without quantum coherence, i.e. pointer states with non-trivial dynamics, have not been discussed before.

For long-lasting oscillations, the system should be coupled weakly to the environment, and the characteristic energies of the system should be much larger than environmental ones. This situation, being less relevant for experiments or the quantum measurement problem, has not been studied in much detail. Usually, in this case (referred to as the quantum limit of decoherence qulim ), the pointer states are the eigenstates $`|\varphi _n`$ of the system’s Hamiltonian. The decay of the non-diagonal elements $`\varphi _n|\rho |\varphi _m`$ of the density matrix $`\rho `$ is Gaussian, and its rate is linearly proportional to the magnitude of the interaction Hamiltonian qulim . However, as we show below, the quantum limit of decoherence can be more subtle.

Analytical studies of dephasing/decoherence are not always possible. In many cases, they include approximations which can be quite stringent (e.g. Markovian behavior of the bath). In this work, we solve directly hans  the compound “system-plus-bath” time-dependent Schrödinger equation, and the analytical approximations we use are checked against this exact numerical solution.

First, let us study a single central spin $`𝐬`$ ($`s=1/2`$) interacting with a bath of spins $`𝐈_k`$ ($`I_k=1/2`$), $`k=1,\mathrm{}N`$. The corresponding Hamiltonian is

$`=_0+𝒱=_S+_B+𝒱,`$ (1)
where $`_S`$ and $`_B`$ are the Hamiltonians of the system (central spin) and the bath, correspondingly, and $`𝒱`$ is the system-bath interaction. The central spin is subjected to an external field $`2\mathrm{\Delta }`$ applied along the $`x`$-axis, i.e. $`_S=2s_x\mathrm{\Delta }=\sigma _x\mathrm{\Delta }`$ ($`\sigma _{x,y,z}`$ are Pauli’s matrices). We consider an Ising-type interaction $`𝒱=\sigma _z_{k=1}^NJ_k\sigma _k^z`$ (where $`\sigma _k^z`$ is Pauli’s matrix of the spin $`𝐈_k`$), and assume that the Hamiltonian of the bath is zero, i.e. the bath has no internal dynamics. Initially, the system and the bath are in an uncorrelated product state $`|s|b`$; the state of the bath $`|b`$ is a superposition of all possible basis states with random coefficients. This model can be used for description of decoherence in various systems, from electron spins to SQUIDs stamprok ; the initial conditions then correspond to the temperature $`\mathrm{\Delta }TJ_k`$.

In Fig. 1, we show the results for $`N=14`$ bath spins; the coupling constants are randomly distributed from zero to $`J_k^{\text{max}}=0.125`$: $`J_k=`${0.123, 0.06425, 0.079, 0.009, 0.0585, 0.03525, 0.012, 0.00525, 0.0945, 0.049, 0.1105, 0.02575, 0.07625, 0.11225}. Initial state of the central spin is defined by the values $`\sigma _x(0)=0.447`$, $`\sigma _y(0)=0`$, $`\sigma _z(0)=0.894`$; the external field $`\mathrm{\Delta }=4.0`$. The observable $`s_z(t)`$ demonstrates damped oscillations: the central spin precesses around the $`x`$-axis with gradually decreasing amplitude. Only $`\sigma _z(t)`$ is shown: oscillations of $`\sigma _y`$ are identical (only shifted in phase), and $`\sigma _x`$ stays constant with good precision. To characterize the dynamics of decoherence, we calculate the quadratic entropy $`S_e(t)=1Tr\rho ^2(t)`$ decgen , Fig. 1(b). It is more convenient than von Neumann’s entropy, while both characterize the same, how strongly mixed is the state of the system.

Initially the amplitude of oscillations drops, and simultaneously the entropy rises very fast. But, after the entropy has come close to its maximum, and the system has decayed into a mixture of pointer states $`|p_1`$ and $`|p_2`$, the oscillations are clearly seen for a very long time with the amplitude of about 20% of the initial value. I.e., long-lasting quantum oscillations can indeed exist in a decohered system. It means that the pointer states $`|p_1`$ and $`|p_2`$ have non-trivial dynamics.

Precision of the simulations can be checked since this model is exactly solvable. The evolution operator is $`U(t)=\mathrm{cos}\mathrm{\Omega }ti(\sigma _zB+\sigma _x\mathrm{\Delta })\mathrm{\Omega }^1\mathrm{sin}\mathrm{\Omega }t`$, where $`\mathrm{\Omega }=[\mathrm{\Delta }^2+B^2]^{1/2}`$, and $`B=J_kS_k^z`$. Bath spins are static, and the bath state $`|b`$ is a random superposition of a large number ($`2^{14}1.610^4`$) of the eigenstates of the operator $`B`$. Therefore, by the central limit theorem, the trace over the bath spins is equivalent to averaging over the Gaussian random field $`B`$ with zero average and dispersion $`b^2_kJ_k^2`$, and for the case $`\mathrm{\Delta }b`$:

$`\sigma _z(t)`$ $`=`$ $`{\displaystyle \frac{\sigma _z(0)}{2\pi b^2}}Re{\displaystyle 𝑑B\mathrm{exp}\left[2i\mathrm{\Delta }tiB^2t/\mathrm{\Delta }\right]}`$
$`\times `$ $`\mathrm{exp}[(1/2)B^2/b^2]`$
The envelope of the $`\sigma _z(t)`$ oscillations is $`\sigma _z^{\text{(env)}}(t)=\sigma _z(0)\left[1+(2t/\tau _1)^2\right]^{1/4}`$, where $`\tau _1=\mathrm{\Delta }/b^2`$. Thus, initially there is the usual Gaussian damping (quadratic with time), but afterwards it changes to a slow power-law decay $`1/\sqrt{2t/\tau _1}`$. The analytical results for the envelope, and for the oscillations of $`\sigma _z(t)`$ themselves, coincide perfectly with the numerical simulations, see Fig. 1(c) and (d).

The analytical solution makes clear that the rapidness of the system’s dynamics ($`\mathrm{\Delta }J_k`$ in our case) is needed for long-time oscillations in a decohered system. Fast motion of the system eliminates from the evolution operator the decohering terms which are of first order in the system-bath coupling $`B`$, and only second-order terms survive, as in Eq. Quantum oscillations without quantum coherence. In contrast to the conclusion of Ref. qulim , the decoherence remains rather fast, although the spin-bath interaction does not have diagonal elements in the $`\sigma _z`$ basis.

We have considered above a static environment, but slow internal dynamics does not influence our qualitative conclusions. As an example, we add an external field acting on the bath spins, so the Hamiltonian of the bath becomes $`_B=h_x_kS_k^x`$, where $`h_x\mathrm{\Delta }`$ (all other parameters are kept the same). The results of simulations are shown in Fig. 2. Again, the slowly decaying oscillations are seen for long times, long after the entropy came almost to saturation, but the oscillations decay faster than in the case of zero $`h_x`$. Rather interestingly, the decay of the $`\sigma _z(t)`$ oscillations is practically independent of the $`h_x`$ value for $`h_x`$ varying from 0.005 to 1.0 (i.e. by more than two orders of magnitude).

An exact analytical solution is not available for $`h_x0`$. To analyze the dynamics of the system, we use the cumulant expansion for the evolution operator (the Magnus expansion, widely used e.g. in the theory of NMR)waugh : $`U(t)=\mathrm{exp}(i_0t)\mathrm{exp}[iX(t)]`$, where

$$X(t)=_0^t𝑑s𝒱(s)\frac{i}{2}_0^t_0^s𝑑s𝑑u[𝒱(s),𝒱(u)]+\mathrm{}$$
(3)
where $`𝒱(t)=\mathrm{exp}(i_0t)𝒱\mathrm{exp}(i_0t)`$, and only secular terms (growing linearly with time without oscillations) are retained. At short times $`t\tau _2=1/h_x`$, the bath is practically static, and the results are the same as for the case $`h_x=0`$. For long times $`t\tau _2`$, independently of the value of $`h_x`$,

$$U(t)=\mathrm{exp}\left[i\sigma _xt\mathrm{\Delta }i\sigma _xt(B_z^2+B_y^2)/(4\mathrm{\Delta })\right],$$
(4)
where $`B_z=_kJ_kS_k^z`$, and $`B_y=_kJ_kS_k^y`$. To evaluate the evolution operator (4), we consider the bath in the mean-field manner similar to the Mermin model mermin , and replace $`B_{y,z}`$ by Gaussian random fields. The envelope is $`\sigma _z^{\text{(env)}}(t)=\sigma _z(0)\left[1+(b^2t/\mathrm{\Delta })^2\right]^{1/2}`$, i.e. the oscillations decay as $`1/t`$. These conclusions agree reasonably with the results of numerical simulations, see Fig. 2(a): for $`h_x`$ varying from $`510^3`$ to 1.0, the shape of the envelope of $`\sigma _z(t)`$ does not depend on $`h_x`$ and is close to the analytical curve.

Another way of turning on the dynamics of the bath is to couple bath spins with each other by small exchange interactions. We have studied the Ising-type exchange $`_B=_kA_{kl}S_k^xS_l^x`$, with constant $`A_{kl}=A`$ and random $`A_{kl}`$. The results are qualitatively the same as above, what agrees with our considerations and arguments of Ref. qulim , which suggest that the exact dynamics of the bath is unimportant as long as it is slow.

Finally, we study the quantum oscillations in decohered many-spin systems. In this case, we see most clearly the two stages, the first one associated with the decay of the system into a mixture of dynamical pointer states (when the entropy rises close to its maximum), and the second one, when the oscillations of the pointer states decay.

As a simple example, we consider two spins $`𝐬_1`$ and $`𝐬_2`$ with Heisenberg anisotropic coupling $`J`$ between them, so that the system’s Hamiltonian $`_0=2J𝐬_1𝐬_2+J/2`$. Spins are coupled to the static bath via isotropic Heisenberg exchange, $`𝒱=_kJ_k(𝐬_1+𝐬_2)𝐒_k`$. For the simulations presented in Fig. 3, we use $`J=8.0`$, and the values of the coupling constants are the same as above. The initial state of the system is the product $`||`$, i.e. the symmetric superposition of the triplet $`|s=1,s_z=0`$ and the singlet $`|s=0`$ states. We present only the oscillations of the $`z`$-component of the first spin $`\sigma _1^z(t)`$: oscillations of $`\sigma _2^z(t)`$ are just shifted in phase by $`\pi `$, and all other components of the central spins remain practically constant.

The results of the numerical simulations are shown in Fig. 3 (note the time scale in this case). Again, there is an initial sharp decrease of the amplitude of oscillations and an increase of entropy; moreover, the correlations $`C_{12}^{\alpha \beta }=\sigma _1^\alpha \sigma _2^\beta `$ between the two central spins follow the behavior of the entropy supporting the conclusion about practically complete initial decoherence of the system. Nevertheless, again, the oscillations exhibit a long tail, although the entropy and the correlations has saturated. Analytical treatment of the Heisenberg case also supports this picture. The exact evolution operator is

$$U(t)=\mathrm{exp}[iJ𝐬^2t/2i𝐬\underset{k}{}J_k𝐒_k],$$
(5)
and its mean-field evaluation gives $`\sigma _1^{z\text{(env)}}(t)=(1/3)[12(b^2t^21)\mathrm{exp}(b^2t^2/2)]`$, i.e. usual fast Gaussian decay, but the final value is 1/3 rather than zero. This prediction, valid for short times, is in very good agreement with the numerical results, see Fig. 3 (d). However, there is a subsequent decay of oscillations, associated with the next-order terms, quadratic in $`b`$. We did not manage to obtain an analytical expression, but in analogy with the single-spin $`h_x0`$ case, we can assume $`1/t`$ decay, and indeed, the numerical results in Fig. 3(d) are in good agreement with the $`1/t`$ assumption.

Detailed analysis of correlations between the two spins allows us to reconstruct the full density matrix, and to check the structure of pointer states explicitly. The results are shown in Fig. 4. At the first stage, the diagonal elements of the density matrix change: part of the $`|s=1,s_z=0`$ spectral weight is transferred equally to the $`|s=1,s_z=\pm 1`$ states. The direction of the system’s total spin randomizes (although incompletely) due to rotation around the randomly oriented effective field generated by the bath (see Eq. 5). Also, at the first stage the non-diagonal element of the density matrix $`s=0|\rho |s=1,s_z=0`$ decays rapidly, not to zero, but to 0.3. The other non-diagonal elements remain small (not shown). Thus, by the end of the first stage, the system has decayed into a mixture of pointer states belonging to the subspaces of different $`s_z`$. However, these pointer states are not static: there are long-living oscillations within the subspace $`s_z=0`$ (containing two states, one singlet and one triplet). These oscillations do not take place inside decoherence-free subspaces (DFS) dfs ; dfs1 , and are due to non-trival fast quantum dynamics of the pointer states. This picture holds also for other initial conditions, but, in general, the dynamics inside the DFS and the oscillations of pointer states show up simultaneously.

The conclusions presented here are confirmed (results not shown) by calculations for different sets of the system’s parameters, for the environments of different sizes, etc. We have also checked the case when the central spins are coupled differently with the environment, and the picture remains the same. Our qualitative conclusions also hold for larger central systems, such as a ring of 4 coupled spins.

In summary, we have studied numerically the dephasing and decoherence processes which take place in some generic systems of interacting spins 1/2 coupled to a spin bath. We compare the dynamical increase of entropy of the central system and the decay of quantum oscillations. We found that in some cases quantum oscillations take place long after the entropy came close to saturation, i.e. that the quantum coherent oscillations can exist in decohered systems. The oscillations exhibit long tails which decay with time as $`1/t`$ or $`1/\sqrt{t}`$, and are observable long after the system has decayed into a mixture of pointer states. Therefore, the pointer states can exhibit non- trivial fast dynamics. We have shown this also by direct analysis of the density matrix.

The authors would like to thank A. Melikidze for helpful discussions. This work was partially carried out at the Ames Laboratory, which is operated for the U. S. Department of Energy by Iowa State University under Contract No. W-7405-82 and was supported by the Director of the Office of Science, Office of Basic Energy Research of the U. S. Department of Energy. Support from the Dutch “Stichting Nationale Computer Faciliteiten (NCF)” is gratefully acknowledged. This work was partially supported by Russian Basic Research Foundation, grant 01-02-16108.


text:
# 1 Introduction

## 1 Introduction

Light Neutrino may play a relevant role in Hot Dark Matter models within a hot-cold dark matter (HCDM) scenario. Their clustering in Galactic, Local dark halos offer the possibility to overcome the Cosmic Black Body opacity ($`>410^{19}eV`$) (GZK) at highest energy cosmic ray astrophysics. These rare events almost in isotropic spread are probably originated by blazars AGN, QSRs in standard scenario, and they should not come, if originally of hadronic nature, from large distances (above tens Mpcs) because of the electro-magnetic dragging friction of cosmic 2.75 K BBR and of the inter-galactic radio backgrounds (GZK cut-off). Indeed as noted by Greisen, Zatsepin and Kuzmin (K.Greisen 1966, Zat’sepin 1966), proton and nucleon mean free path at E $`>510^{19}EeV`$ is less than 30 $`Mpc`$ and asymptotically nearly ten $`Mpc`$.; also gamma rays at those energies have even shorter interaction length ($`10Mpc`$) due to severe opacity by electron pair production via microwave and radio background interactions (R.J.Protheroe 1997). Nevertheless these powerful sources (AGN, Quasars, GRBs) suspected to be the unique source able to eject such UHECRs, are rare at nearby distances ($`<10÷20Mpc`$, as for nearby $`M87`$ in Virgo cluster); moreover there are not nearby $`AGN`$ in the observed UHECR arrival directions. Strong and coherent galactic(R.J.Protheroe 1997) or extragalactic (Farrar et all. 2000) magnetic fields, able to bend such UHECR (proton, nuclei) directions are not really at hand. The needed coherent lengths and strength are not easily compatible with known cosmic magnetic fields. Finally in this scenario the $`ZeV`$ neutrons born, by photo-pion proton conversions on BBR, may escape the magnetic fields bending and should keep memory of the arrival direction, leading to (unobserved) clustering toward the primary source. Secondaries EeV photons (by neutral pion decays) should also abundantly point and cluster toward the same nearby $`AGN`$ sources (P.Bhattacharjee 2000, Elbert et all. 1995) contrary to $`AGASA`$ data. Another solution of the present GZK puzzle, the Topological defects ($`TD`$), assumes as a source, relic heavy particles of early Universe; they are imagined diffused as a Cold Dark Matter component, in galactic or Local Group Halos. Nevertheless the $`TD`$ fine tuned masses and ad-hoc decays are unable to explain the growing evidences of doublets and triplets clustering in $`AGASA`$ $`UHECR`$ arrival data. On the other side there are growing evidences of self-correlation between UHECR arrival directions with far Compact Blazars at cosmic distance well above GZK cut-off (Tinyakov P.G.et Tkachev 2001). Therefore the solution of UHECR puzzle based on primary Extreme High Energy (EHE) neutrino beams (from AGN) at ZeV $`E_\nu >10^{21}`$ eV and their undisturbed propagation from cosmic distances up to nearby calorimeter (made by relic light $`\nu `$ in dark galactic or local dark halo (Fargion et Salis 1997, Fargion, Mele et Salis 1999, Weiler 1999, Yoshida et all. 1998) is still, the most favorite convincing solution for the GZK puzzle. New complex scenarios for each neutrino mass spectra are then opening and important signature of UHECR Z,WW showering must manifest in observed anisotropy and space-time clustering.

## 2 Relic $`\nu _r`$ neutrino masses and Hot Halo Clustering

If relic neutrinos have a mass larger than their thermal energy (1.9 $`K^0`$) they may cluster in galactic or Local Group halos; at eVs masses the clustering seem very plausible and it may play a role in dark hot cosmology(Fargion 1983). Their scattering with incoming extra-galactic EHE neutrinos determine high energy particle cascades which could contribute or dominate the observed UHECR flux at $`GZK`$ edges. Indeed the possibility that neutrino share a little mass has been reinforced by Super-Kamiokande evidence for atmospheric neutrino anomaly via $`\nu _\mu \nu _\tau `$ oscillation. An additional evidence of neutral lepton flavour mixing has been very recently reported also by Solar neutrino experiment (SNO,Gallex,K2K). Consequently there are at least two main extreme scenario for hot dark halos: either $`\nu _\mu ,\nu _\tau `$ are both extremely light ($`m_{\nu _\mu }m_{\nu _\tau }\sqrt{(\mathrm{\Delta }m)^2}0.05eV`$) and therefore hot dark neutrino halo is very wide, smeared and spread out to local group clustering sizes (increasing the radius but loosing in the neutrino density clustering contrast), or $`\nu _\mu ,\nu _\tau `$ may share degenerated ($`eV`$ masses) split by a very tiny different values. In the latter fine-tuned neutrino mass case ($`m_\nu 0.4eV1.2eV`$) (see Fig,2 and Fig.3) the Z peak $`\nu \overline{\nu }_r`$ interaction (Fargion et Salis 1997, Fargion, Mele et Salis 1999, Weiler 1999, Yoshida et all. 1998) will be the favorite one; in the second case (for heavier non constrained neutrino mass ($`m_\nu >3eV`$)) only a $`\nu \overline{\nu }_rW^+W^{}`$(Fargion et Salis 1997, Fargion, Mele et Salis 1999), and the additional $`\nu \overline{\nu }_rZZ`$ interactions, (see the cross-section in Fig.1)(Fargion et all. 2001) considered here will be the only ones able to solve the GZK puzzle. Indeed the relic neutrino mass within HDM models in galactic halo near $`m_\nu 4eV`$, corresponds to a lower and $`Z`$ resonant incoming energy

$$E_\nu =\left(\frac{4eV}{\sqrt{m_{\nu }^{}{}_{}{}^{2}+p_\nu ^2}}\right)10^{21}eV.$$

This resonant incoming neutrino energy is unable to overcome GZK energies while it is showering mainly a small energy fraction into nucleons ($`p,\overline{p},n,\overline{n}`$), (see $`Tab.1`$ below), at energies $`E_p`$ quite below. (see $`Tab.2`$ below).

$$E_p=2.2\left(\frac{4eV}{\sqrt{m_{\nu }^{}{}_{}{}^{2}+p_\nu ^2}}\right)10^{19}eV.$$

Therefore too heavy ($`>1.5eV`$) neutrino mass are not fit to solve GZK by Z-resonance; on the contrary WW,ZZ showering as well as t-channel showering may naturally keep open the solution. In particular the overlapping of both the Z and the WW, ZZ channels described in fig.1, for $`m_\nu 2.3eV`$ while solving the UHECR above GZK they must pile up (by Z-resonance peak activity) events at $`510^{19}eV`$, leading to a bump in AGASA data. There is indeed a first marginal evidence of such a UHECR bump in AGASA and Yakutsk data that may stand for this interpretation. More detailed data are needed to verify such very exciting possibility. Similar result regarding the fine tuned relic mass at $`0.4eV`$ and $`2.3eV`$, (however ignoring the WW ZZ and t-channels and invoking very hard UHE neutrino spectra) have been independently reported recently (Fodor et all. 1999) . Most of us consider cosmological light relic neutrinos in Standard Model at non relativistic regime neglecting any relic neutrino momentum $`p_\nu `$ term. However, at lightest mass values the momentum may be comparable to the relic mass; moreover the spectra may reflect additional relic neutrino-energy injection which are feeding standard cosmic relic neutrino at energies much above the same neutrino mass. Indeed there may exist, within or beyond Standard Cosmology, a relic neutrino component due to stellar, Super Nova, GRBs, AGN past activities, presently red-shifted into a KeV-eV spectra, piling into a relic neutrino grey-body spectra. Therefore it is worth-full to keep the most general mass and momentum term in the target relic neutrino spectra. In this windy ultra-relativistic neutrino cosmology, eventually leading to a neutrino radiation dominated Universe, the halo size to be considered is nearly coincident with the GZK one defined by the energy loss lenght for UHECR nucleons ($`20Mpcs`$). Therefore the isotropic UHECR behaviour is guaranteed but a puzzle related to uniform source distribution seem to persist. Nevertheless the UHE neutrino- relic neutrino scattering do not follow a flat spectra as shown in figure 2, (as well as any hypothetical $`\nu `$ grey body spectra). This leave open the opportunity to have a relic relativistic neutrino component at eVs energies as well as the observed non uniform UHECR spectra. This case is similar to the case of a very light neutrino mass much below $`0.1`$ eV. As we noticed above, relic neutrino mass above a few eVs in HDM halo are not consistent with naive Z peak; higher energies interactions ruled by WW,(D.Fargion, B.Mele et A.Salis 1999; K.Enqvist et all. 1989) ZZ cross-sections (Fargion 2001) may nevertheless solve the GZK cut-off. In this regime there will be also possible to produce by virtual W exchange, t-channel, $`UHE`$ lepton pairs, by $`\nu _i\overline{\nu }_jl_i\overline{l}_j`$, leading to additional electro-magnetic showers injection.
As we shall see these important and underestimated signal will produce UHE electrons whose final trace are TeVs synchrotron photons able to break the IR-TeV cosmic cut-off. The hadronic tail of the Z or $`W^+W^{}`$ cascade maybe the source of final nucleons $`p,\overline{p},n,\overline{n}`$ able to explain UHECR events observed by Fly’s Eye and AGASA (Y.Uchihori et al. 2000) and other detectors. The same $`\nu \overline{\nu }_r`$ interactions are source of Z and W that decay in rich shower ramification. The average energy deposition for both gauge bosons among the secondary particles is summarized in Table 1A below.

## 3 UHECR channels in Z showers

Although protons (or anti-protons) are the most popular and favorite candidate in order to explain the highest energy air shower observed, one doesn’t have to neglect the signature of final neutron and anti-neutrons as well as electrons and photons. Indeed the UHECR neutrons are produced in Z-WW showering at nearly same rate as the charged nucleons. Above GZK cut-off energies UHE $`n`$,$`\overline{n}`$, share a life lenght comparable with the Hot Galactic Dark Neutrino Halo. Therefore they may be an important component in UHECRs. Moreover prompt UHE electron (positron) interactions with the galactic or extra-galactic magnetic field or soft radiative backgrounds may lead to gamma cascades and from PeVs to TeVs energies.
Gamma photons at energies $`E_\gamma 10^{20}`$ \- $`10^{19}eV`$ may freely propagate through galactic or local halo scales (hundreds of kpc to few Mpc) and could also contribute to the extreme edges of cosmic ray spectrum and clustering (Yoshida et al. 1998, Fargion et al. 2001).
The ratio of the final energy flux of nucleons near the Z peak resonance, $`\mathrm{\Phi }_p`$ over the corresponding electro-magnetic energy flux $`\mathrm{\Phi }_{em}`$ ratio is, as in tab.1 $`e^+e^{},\gamma `$ entrance, nearly $`\frac{1}{8}`$. Moreover if one considers at higher $`E_\nu `$ energies, the opening of WW, ZZ channels and the six pairs $`\nu _e\overline{\nu _\mu }`$, $`\nu _\mu \overline{\nu _\tau }`$, $`\nu _e\overline{\nu _\tau }`$ (and their anti-particle pairs) t-channel interactions leading to highest energy leptons, with no nucleonic relics (as $`p,\overline{p}`$), this additional injection favors the electro-magnetic flux $`\mathrm{\Phi }_{em}`$ over the corresponding nuclear one $`\mathrm{\Phi }_p`$ by a factor $`1.6`$ leading to $`\frac{\mathrm{\Phi }_p}{\mathrm{\Phi }_{em}}\frac{1}{13}`$. This ratio is valid at $`WW,ZZ`$ masses because the overall cross section variability is energy dependent. At center of mass energies above these values, the $`\frac{\mathrm{\Phi }_p}{\mathrm{\Phi }_{em}}`$ decreases more because the dominant role of t-channel (Fig1). We focus here on Z, and WW,ZZ channels showering in hadrons for GZK events. The important role of UHE electron showering into TeV radiation is discussed below.

There is an upper bound density clustering for very light Dirac fermions due to the maximal Fermi degenerancy whose adimensional density contrast is $`\delta \rho m_\nu ^3`$, while one finds (Fargion 1983) that the neutrino free-streaming halo grows only as $`m_\nu ^1`$. Therefore the overall interaction probability grows $`m_\nu ^2`$, favoring heavier non relativistic (eVs) neutrino masses. In this frame above few eV neutrino masses only WW-ZZ channel are operative. Nevertheless the same lightest relic neutrinos may share higher Local Group velocities (thousands $`\frac{Km}{s}`$) or even nearly relativistic speeds and it may therefore compensate the common density bound:

$$n_{\nu _i}=1.910^3\left(\frac{m_i}{0.1eV}\right)^3\left(\frac{v_{\nu _i}}{210^3\frac{Km}{s}}\right)^3$$
(1)
From the cross section side there are three main interaction processes that have to be considered leading to nucleons in the of EHE and relic neutrinos scattering.

channel 1. The $`\nu \nu _rZ`$ annihilation at the Z resonance.

channel 2. $`\nu _\mu \overline{\nu _\mu }W^+W^{}`$ or $`\nu _\mu \overline{\nu _\mu }ZZ`$ leading to hadrons, electrons, photons, through W and Z decay.

channel 3. The $`\nu _e`$ \- $`\overline{\nu _\mu }`$, $`\nu _e`$ \- $`\overline{\nu _\tau }`$, $`\nu _\mu `$ \- $`\overline{\nu _\tau }`$ and antiparticle conjugate interactions of different flavor neutrinos mediated in the $`t`$-channel by the W exchange (i.e. $`\nu _\mu \overline{\nu _{\tau _r}}\mu ^{}\tau ^+`$). These reactions are sources of prompt and secondary UHE electrons as well as photons resulting by hadronic $`\tau `$ decay.

Their cross-section values are plotted in Fig.1. The asymptotic behaviour of these cross section is proportional to $`(\frac{M_W^2}{s})\mathrm{ln}(\frac{s}{M_W^2})`$ for $`sM_Z^2`$.
The nucleon arising from WW and ZZ hadronic decay could provide a reasonable solution to the UHECR events above GZK. We’ll assume that the fraction of pions and nucleons related to the total number of particles from the W boson decay is the almost the same of Z boson. So W hadronic decay ($`P0.68`$) leads on average to about 37 particles, where $`<n_{\pi ^0}>9.19`$, $`<n_{\pi ^\pm }>17`$, and $`<n_{p,\overline{p},n,\overline{n}}>2.7`$. In addition we have to expect by the subsequent decays of $`\pi `$’s (charged and neutral), kaons and resonances ($`\rho `$, $`\omega `$, $`\eta `$) produced, a flux of secondary UHE photons and electrons. As we already pointed out, the particles resulting from the decay are mostly prompt pions. The others are particles whose final decay likely leads to charged and neutral pions as well. As a consequence the electrons and photons come from prompt pion decay. On average it results (Fargion, Grossi et Lucentini 2001; Fargion, Grossi, Lucentini et Troia 2001) that the energy in the bosons decay is not uniformly distributed among the particles. Each charged pion will give an electron (or positron) and three neutrinos, that will have less than one per cent of the initial W boson energy, while each $`\pi ^0`$ decays in two photons, each with 1 per cent of the initial W energy. In the Table 1A above we show all the channels leading from single Z,W and Z pairs as well as t-channel in nuclear and electro-magnetic components.

This interactions, as noted in Table 1A are leading to electro-magnetic showers and are not offering any nuclear secondary.

## 4 The Boosted Z-UHECR spectra

Let us examine the destiny of UHE primary particles (nucleons, electrons and photons) ($`E_e<10^{21}eV`$) produced after hadronic or leptonic W decay. As we already noticed in the introduction, we’ll assume that the nucleons, electrons and photons spectra (coming from W or Z decay) after $`\nu \nu `$ scattering in the halo, follow a power law that in the center of mass system is $`\frac{dN^{}}{dE^{}dt^{}}E^\alpha `$ where $`\alpha 1.5`$. This assumption is based on detailed Monte Carlo simulation of a heavy fourth generation neutrino annihilations(Yu. A.Golubkov 1998; D. Fargion,Yu.A.Golubkov, M.Yu.Khlopov 1999; D. Fargion, R. Konoplich et all. 2000) and with the model of quark - hadron fragmentation spectrum suggested by Hill (C.T.Hill 1983).

In order to determine the shape of the particle spectrum in the laboratory frame, we have to introduce the Lorentz relativistic transformations from the center of mass system to the laboratory system. The number of particles is clearly a relativistic invariant $`dN_{lab}=dN^{}`$, while the relation between the two time intervals is $`dt_{lab}=\gamma dt^{}`$, the energy changes like $`ϵ_{lab}=\gamma ϵ^{}(1+\beta \mathrm{cos}\theta ^{})=ϵ^{}\gamma ^1(1\beta \mathrm{cos}\theta )^1`$, and finally the solid angle in the laboratory frame of reference becomes $`d\mathrm{\Omega }_{lab}=\gamma ^2d\mathrm{\Omega }^{}(1\beta \mathrm{cos}\theta )^2`$. Substituting these relations one obtains

$$\left(\frac{dN}{dϵdtd\mathrm{\Omega }}\right)_{lab}=\frac{dN_{}}{dϵ_{}dt_{}d\mathrm{\Omega }_{}}\gamma ^2(1\beta \mathrm{cos}\theta )^1=$$

$$=\frac{ϵ_{}^\alpha \gamma ^2}{4\pi }(1\beta \mathrm{cos}\theta )^1=\frac{ϵ^\alpha \gamma ^{\alpha 2}}{4\pi }(1\beta \mathrm{cos}\theta )^{\alpha 1}$$
(2)
and integrating on $`\theta `$ (omitting the lab notation) one loses the spectrum dependence on the angle.

The consequent fluence derived by the solid angle integral is:

$$\frac{dN}{dϵdt}ϵ^2=\frac{ϵ^{\alpha +2}\gamma ^{\alpha 2}}{2\beta \alpha }[(1+\beta )^\alpha (1\beta )^\alpha ]\frac{2^{\alpha 1}ϵ^{\alpha +2}\gamma ^{\alpha 2}}{\alpha }$$
(3)
There are two extreme case to be considered: the case where the interaction occurs at Z peak resonance and therefore the center of mass Lorents factor $`\gamma `$ is frozen at a given value (eq.1) and the case (WW,ZZ pair channel) where all energies are allowable and $`\gamma `$ is proportional to $`ϵ^{1/2}`$. Here we focus only on Z peak resonance. The consequent fluence spectra $`\frac{dN}{dϵdt}ϵ^2`$, as above, is proportional to $`ϵ^{\alpha +2}`$. Because $`\alpha `$ is nearly $`1.5`$ all the consequent secondary particles will also show a spectra proportional to $`ϵ^{1/2}`$ following a normalized energies shown in Tab.2, as shown in Fig.(2-6). In the latter case (WW,ZZ pair channel), the relativistic boost reflects on the spectrum of the secondary particles, and the spectra power law becomes $`ϵ^{\alpha /2+1}=ϵ^{0.25}`$. These channels will be studied in details elsewhere. In Fig. 1 we show the spectrum of protons, photons and electrons coming from Z hadronic and leptonic decay assuming a nominal primary CR energy flux $`20eVs^1sr^1cm^2`$, due to the total $`\nu \overline{\nu }`$ scattering at GZK energies as shown in figures 2-6. Let us remind that we assume an interaction probability of $`1\%`$ and a corresponding UHE incoming neutrino energy $`2000eVs^1sr^1cm^2`$ near but below present $`UHE`$ neutrino flux bound from AMANDA and Baikal as well as Goldstone data.

## 5 The UHECRs from Relic $`\nu `$ Masses

The role of each relic neutrino mass is summirized from the convolutions of the UHE neutrino spectra with the relic neutrino mass, its density as well as the cross-sections described above. The case of Z-resonance event with a single neutrino mass has a narrow fine tuned energy mass windows (0.4 eV-1.2 eV) described respectively in Figures 3-4.

We remind again that a heavier neutrino mass $`(2eVs)`$ imply the rise of WW-ZZ channels and a pile up of Z resonance cross-section at lower UHECR spectra. This feature maybe already responsible for the tiny bump in observed events around $`510^{19}eV`$. The lighter neutrino mass possibilities (near 0.1 eV) are comparable with present Super-Kamiokande atmospheric neutrino mass and are leading to the exciting scenario where more non degenerated Z-resonances occur (Fargion 2001). These scenario are summarized in Fig. 5 (for nominal example $`m_{\nu _\tau }`$ = 0.1 eV; $`m_{\nu _\mu }`$ = 0.05 eV). The twin neutrino mass inject a corresponding twin bump at highest energy. Another limiting case of interest takes place when the light neutrino masses are extreme, nearly at atmospheric (SK,K2K) and solar (SNO) neutrino masses. This case is described in two different versions in Fig.6 (assuming comparable neutrino densities) and Fig.7 (keeping care of the lightest neutrino density diluitions). The relic neutrino masses are assumed $`m_{\nu _\tau }`$ = 0.05 eV; $`m_{\nu _\mu }`$ = 0.001 eV). A more complex scenario is also possible when it takes place both a narrow twin bump (Fig5) $`and`$ a wider twin bump (Fig 6-7) because of a small neutrino tau-muon mass splitting overlapping with a wider one due to lightest neutrino electron mass.

## 6 UHECRs Clustering and Anisotropy

The neutrino mass play a role in defining its Hot Dark Halo size and the consequent enhancement of UHECR arrival directions due to our peculiar position in the HDM halo. Indeed for a heavy $`2eV`$ mass case HDM neutrino halo are mainly galactic and/or local, reflecting an isotropic or a diffused amplification toward nearby $`M31`$ HDM halo. In the lighter case the HDM should include the Local Cluster up to Virgo. To each size corresponds also a different role of UHECR arrival time. The larger the HDM size the longer the UHECR random-walk travel time (in extra-galactic random magnetic fields) and the wider the arrival rate lag between doublets or triplets. The smaller is the neutrino halo the earlier the UHE neutron secondaries by Z shower will play a role: indeed at $`E_n=10^{20}eV`$ UHE neutron are flying a Mpc and their directional arrival (or their late decayed proton arrival) are more on-line toward the source. This may explain the high self collimation and auto-correlation of UHECR discovered very recently (Tinyakov et all. 2001). The UHE neutrons Z-Showering fits with the harder spectra observed in clustered events in AGASA (Takeda et all.2001). The same UHECR alignment may explain the quite short (2-3 years)(Takeda 2001) lapse of time observed in AGASA doublets. Indeed the most conservative scenario where UHECR are just primary proton from nearby sources at GZK distances (tens of Mpcs) are no longer acceptable either because the absence of such nearby sources and because of the observed stringent UHECR clustering ($`2^o2.5^o`$) (Takeda et all. 2001) in arrival direction, as well as because of the short ($`3`$ years) characteristic time lag between clustered events. Finally the same growth with energy of UHECR neutron (and anti-neutron) life-lengths (while being marginal or meaning-less in tens Mpcs GZK flight distances) may naturally explain, within a $`Mpc`$ Z Showering Neutrino Halo, the arising harder spectra revealed in doublets-triplet spectra (Aoki et all. 2001).

## 7 The apparent Tinyakov-Glushkov Paradox

The same role of UHE neutron secondaries from Z showering in HDM halo may also solve an emerging puzzle: the correlations of arrival directions of UHECRs found recently (Glushkov et all. 2001) in Yakutsk data at energy $`E=810^{18}eV`$ toward the Super Galactic Plane are to be compared with the compelling evidence of UHECRs events ($`E=310^{19}eV`$ above GZK) clustering toward well defined BL Lacs at cosmic distances (redshift $`z>0.10.2`$) (Tinyakov et Tkachev 2001; Tinyakov et al. 2001). Where is the real UHECR sources location? At Super-galactic disk (50 Mpcs wide, within GZK range) or at cosmic ($`300Mpcs`$) edges? It should be noted that even for the Super Galactic hypothesis (Glushkov et all. 2001) the common proton are unable to justify the high collimation of the UHECR events. Of course both results (or just one of them) maybe a statistical fluctuation. But both studies seem statistically significant (4.6-5 sigma) and they seem in obvious disagreement. There may be still open the possibility of $`two`$ new categories of UHECR sources both of them located at different distances above GZK ones (the harder the most distant BL Lac sources). But it seem quite unnatural the UHECR propagation by direct nucleons where the most distant are the harder. However our Z-Showering scenario offer different solutions: (1) The Relic Neutrino Masses define different Hierarchical Dark Halos and privileged arrival direction correlated to Hot Relic Neutrino Halos. The real sources are at (isotropic) cosmic edges (Tinyakov et Tkachev 2001; Tinyakov et al. 2001), but their crossing along a longer anisotropic relic neutrino cloud enhance the interaction probability in the Super Galactic Plane. (2) The nearest SG sources are weaker while the collimated BL Lacs are harder: both sources need a Neutrino Halo to induce the Z-Showering UHECRs. More data will clarify better the real scenario.

## 8  The TeV Tails from UHE electrons

As it is shown in Table 1A-B and Figures above, the electron (positron) energies by $`\pi ^\pm `$ decays is around $`E_e210^{19}eV`$ for an initial $`E_Z10^{22}eV`$ ( and $`E_\nu 10^{22}eV`$). Such electron pairs while not radiating efficently in extra-galactic magnetic fields will be interacting with the galactic magnetic field ($`B_G10^6G`$) leading to direct TeV photons:

$$E_\gamma ^{sync}\gamma ^2\left(\frac{eB}{2\pi m_e}\right)$$

$$27.2\left(\frac{E_e}{210^{19}eV}\right)^2\left(\frac{m_\nu }{0.4eV}\right)^2\left(\frac{B}{\mu G}\right)TeV.$$
(4)
The same UHE electrons will radiate less efficiently with extra- galactic magnetic field ($`B_G10^9G`$) leading also to direct peak $`27.2`$ GeV photons. The spectrum of these photons is characterized by a power of law $`dN/dEdTE^{(\alpha +1)/2}E^{1.25}`$ where $`\alpha `$ is the power law of the electron spectrum, and it is showed in Figures above. As regards the prompt electrons at higher energy ($`E_e10^{21}eV`$), in particular in the t-channels, their interactions with the extra-galactic field first and galactic magnetic fields later is source of another kind of synchrotron emission around tens of PeV energies $`E_\gamma ^{sync}`$:

$$6.810^{13}\left(\frac{E_e}{10^{21}eV}\right)^2\left(\frac{m_\nu }{0.4eV}\right)^2\left(\frac{B}{nG}\right)eV$$
(5)
$$6.810^{16}\left(\frac{E_e}{10^{21}eV}\right)^2\left(\frac{m_\nu }{0.4eV}\right)^2\left(\frac{B}{\mu G}\right)eV$$
(6)
The corresponding energy loss length instead is (O.E.Kalashev, V.A.Kuzmin, D.V.Semikoz 2000)

$$\left(\frac{1}{E}\frac{dE}{dt}\right)^1=3.8\times \left(\frac{E}{10^{21}}\right)^1\left(\frac{B}{10^9G}\right)^2kpc.$$
(7)
For the first case the interaction lenght is few Kpcs while in the second one in few days light flight. Again one has the same power law characteristic of a synchrotron spectrum with index $`E^{(\alpha +1/2)}E^{1.25}`$. Gammas at $`10^{16}÷10^{17}`$ eV scatters onto low-energy photons from isotropic cosmic background ($`\gamma +BBRe^+e^{}`$) converting their energy in electron pair. The expression of the pair production cross-section is:

$$\sigma (s)=\frac{1}{2}\pi r_0^2(1v^2)[(3v^4)\mathrm{ln}\frac{1+v}{1v}2v(2v^2)]$$
(8)
where $`v=(14m_e^2/s)^{1/2}`$, $`s=2E_\gamma ϵ(1\mathrm{cos}\theta )`$ is the square energy in the center of mass frame, $`ϵ`$ is the target photon energy, $`r_0`$ is the classic electron radius, with a peak cross section value at

$$\frac{4}{137}\times \frac{3}{8\pi }\sigma _T\mathrm{ln}183=1.2\times 10^{26}cm^2$$

Because the corresponding attenuation length due to the interactions with the microwave background is around ten kpc, the extension of the halo plays a fundamental role in order to make this mechanism efficient or not. As is shown in Fig.3-4, the contribution to tens of PeV gamma signals by Z (or W) hadronic decay, could be compatible with actual experimental limits fixed by CASA-MIA detector on such a range of energies. Considering a halo extension $`l_{halo}>100kpc`$, the secondary electron pair creation becomes efficient, leading to a suppression of the tens of PeV signal. So electrons at $`E_e3.510^{16}eV`$ loose again energy through additional synchrotron radiation (O.E.Kalashev, V.A.Kuzmin, D.V.Semikozb 2000), with maximum $`E_\gamma ^{sync}`$ around

$$79\left(\frac{E_e}{10^{21}eV}\right)^4\left(\frac{m_\nu }{0.4eV}\right)^4\left(\frac{B}{\mu G}\right)^3MeV.$$
(9)
Anyway this signal is not able to pollute sensibly the MeV-GeV; the relevant signal pile up at TeVs.

Gamma rays with energies up to 20 TeV have been observed by terrestrial detector only by nearby sources like Mrk 501 (z = 0.033) or very recently by MrK 421 (z = 0.031). More recent evidences of tens TeVs from (three times more) distant blazar 1ES1426+428 (z = 0.129) make even more dramatic the IR-TeV cut-off. This is puzzling because the extra-galactic TeV spectrum should be, in principle, significantly suppressed by the $`\gamma `$-rays interactions with the extra-galactic Infrared background, leading to electron pair production and TeVs cut-off. The recent calibration and determination of the infrared background by DIRBE and FIRAS on COBE have inferred severe constrains on TeV propagation. Indeed, as noticed by Kifune (Kifune 1997), and Protheroe and Meyer (Protheroe et Meyer 2000) we may face a severe infrared background - TeV gamma ray crisis. This crisis imply a distance cut-off, incidentally, comparable to the GZK one. Let us remind also an additional evidence for IR-TeV cut-off is related to the possible discover of tens of TeV counterparts of BATSE GRB970417, observed by Milagrito(R. Atkins et all. 2000), being most GRBs very possibly at cosmic edges, at distances well above the IR-TeV cut-off ones. In this scenario it is also important to remind the possibilities that the Fly’s Eye event has been correlated to TeV pile up events in HEGRA (Horns et al. 1999). The very recent report (private communication 2001) of the absence of the signal few years later at HEGRA may be still consistent with a limited UHE TeV tail activity. To solve the IR-TeV cut-off one may alternatively invoke unbelievable extreme hard intrinsic spectra or exotic explanation as gamma ray superposition of photons or sacrilegious Lorentz invariance violation (G. Amelino-Camelia,et all 1998).

## 9 Conclusion

UHECR above GZK may be naturally born by UHE $`\nu `$ scattering on relic ones. The target cosmic $`\nu `$ may be light and dense as the needed ones in HDM model (few eVs). Then their $`W^+W^{},ZZ`$ pair productions channel (not just the Z resonant peak) would solve the GZK puzzle. At a much lighter, but fine tuned case $`m_\nu 0.4eV`$, $`m_\nu 1.5eV`$ assuming $`E_\nu 10^{22}eV`$, one is able to solve at once the known UHECR data at GZK edge by the dominant Z peak; in this peculiar scenario one may foresee (fig.2-3) a rapid decrease (an order of magnitude in energy fluence) above $`310^{20}eV`$ in future data and a further recover (due to WW,ZZ channels) at higher energies. The characteristic UHECR fluxes will reflect the averaged neutrino-neutrino interactions shown in Fig.2-7. Their imprint could confirm the neutrino masses value and relic density. At a more extreme lighter neutrino mass, occurring for $`m_\nu m_{\nu _SK}0.05eV`$, the minimal $`m_{\nu _\tau },m_{\nu _\mu }`$ small mass differences might be reflected, in a spectacular way, into UHECR modulation quite above the GZK edges. The ”twin” lightest masses (Fig.5-6-7) call for either gravitational $`\nu `$ clustering above the expected one or the presence of relativistic diffused background. Possible neutrino gray body spectra, out of thermal equilibrium, at higher energies may also arise from non standard early Universe. The UHECR acceleration is not yet solved, but their propagation from far cosmic volumes is finally allowed. The role of UHE neutrons in Z-showering, their directional flight leading to clustering in self collimated data is possibly emerging by harder spectra. Peculiar secondaries of TeVs tails may be precursor and afterglows signal correlated to past or future UHECRs pointing toward the same far sources. The IR-TeV solution may be just be a necessary corollary of the Z-Showering GZK solution (Fargion, Grossi et Lucentini 2001). The time and space directional may be a new fundamental test of present Z-Showering model. The discover of UHE neutrino at GZK energies might be testify on ground by UHE $`\tau `$ air-shower, born by direct $`10^{19}`$eV UHE $`\nu `$ crossing small Earth crust depth, flashing from the horizontal edges to mountain,balloon and satellite detectors (Fargion 2000). The new generation UHECR data within next decade, may also offer the probe of lightest elementary particle masses, their relic densities, their spatial map distribution and energies and the most ancient and evasive shadows of earliest $`\nu `$ cosmic relic backgrounds.

DISCUSSION

JIM BEALL: For photo-pion production in the AGN sources, where do the pions come from? This is a somewhat rhetorical question.

D.FARGION: In the UHE $`\nu `$-$`\nu `$ scattering to the overcome GZK cut-off we consider two places where pions occurs:

1. Near the AGN source: the UHE photon (or nuclei) hitting BBR (2,7K) photons or local thermal photons lead to: $`p+\gamma \mathrm{}^+\pi ^++n`$
$`p+\gamma \mathrm{}^+\pi ^0+p`$
or
$`p+\gamma N\pi +p`$
$`p+\gamma N\pi +n`$
These UHE pions decays to UHE $`\nu `$. The UHE $`\nu ^+`$ relic $`\nu `$ produce Z (UHE parent of UHECR). The Z decay into 2.7 nucleons and $`30`$ pions (neutral + charged).
The neutral pions decay into UHE $`\gamma `$, the charged ones decay into $`\mu `$ $`e`$, and $`\nu _e`$,$`\nu _\mu `$.

WOLFGANG KUNDT: What are your objections to a very near population of CR sources, d$`Kpc^2`$. (I Like to think of slingshot-acceleration of CR’s by neutron-star magnetospheres).

D.FARGION: I am sorry because I was first understanding your ”near objects” as Farrar-Bierman-Piran solution of AGN near (M-87 - like $``$ few Mpc) + $`B_G`$. They should be also UHE neutron, $`\overline{n}`$ ($`p+\gamma \pi ^++n`$ by photopion), source that should point directly to M-87 (this was not observed). About ”galactic” (Kpc radius) I think they should follow the galactic disk distribution with strong quadrupole modulation in UHECR (unobserved). For HALO source tiny dipole anisotropy might be hidden but present data seem not correlated with known PSR distribution.

S.COLAFRANCESCO: Which is the contribution on $`\mathrm{\Omega }_0`$ by the relic neutrinos you need ($`m_\nu 0.1÷10eV`$)in the dark halo ($`R_{halo}1÷2Mpc`$) in order to solve the GZK problem?

D.FARGION: The amount of dark matter I consider in cosmology is linearly dependent on $`m_\nu `$ mass (for simplicity here I mean an unique flavour mass dominant component)

$$\mathrm{\Omega }_\nu 10^2\left(\frac{m_\nu }{eV}\right)$$

The amount of dark matter in galactic or local group halo has a density which may reach a density contrast:

$$\frac{\delta \varrho }{\varrho }10^2÷10^3\left(\frac{m_\nu }{eV}\right)^\beta $$

model dependent $`0\beta 2`$. The bounds we found up to day are fixing:

$`m_\nu >0.1eV`$(extended extragalactic halo); $`2eV<m_\nu <5eV`$ (galactic halo).