text:
# A Tip of the TOE

## 1 Introduction

The matter of sand has ever since fascinated mankind. Already twenty thousand years ago attempts to classify the grainy zoo have been made, . In the Cretacious Age many advances with respect to sand have been undertaken, especially it has been proved by  that chalk and sand are—at a fundamental level—two manifestations of the same thing. Unfortunately, since the end of the Cretatious Age, sand research activities suffered from a considerable lack of financial resources, a time period which is called the Big Sand Crisis (BSC). This led to the the well-known revolutionary movement “Every Scientist Needs a Camel”, . After all, at least a few new concepts have been invented, among them the so-called wheel, which then, however, has been discarded due to its impracticability, . (However, there are some unteachable fanatics who still try to demonstrate the usefulness of wheels on sand in the annual Paris-Dakar competition.)

Some centuries after the BSC a completely new impulse came from ancient Lei Li Ga, , who was the first to put special interest in the dynamics of sand.

The so far highly disordered efforts in sand have been fibre bundled by the famous philosoph Goe-The, who formulated the central question of today’s complex theory of sand, :

Herauszufinden, was den Sand im Innersten zusammenband.2
22footnotetext: to be translated maybe as: To find out, what sand bound together in its inner.

Starting from this point, many interesting theories about sand dynamics have been developed. One of the most ambitious (and ambiguous) models was the well-known sandstring theory developed by the old Egyptian priest Edua ’Wit X . With the postulate of the non-existence of infinitely small grains of sand, several severe problems of ordinary sand-theory (e.g. incurable sand-divergencies, causality problems etc.) seemed to be repaired. However, over the millennia, this model did not bear any fruits in sand. (After all, there have been some interesting side results for weavers, in particular new techniques for knotting.) This led to a complete restart of research in the field of sand. Fortunately, at least some interesting questions have been solved independently of sandstring theory, , , .

Modern, post-string sand theories have also their origin in Egypt. In particular, a comprehensive quantum theory of the dynamics of sand dunes has been developed, , . Piles of stones (so called “pyramides”) have been constructed for a macroscopic test of the dynamic theory, especially for the prominent dune-tunnel effect. However, the pyramid-models emerged as too roughly textured, so there could not be proven anything with these (although they are still very impressive). Indeed, experimental clues for the correctness of the Egyptian Theory of Quantum Sand Dynamics (ETQSD) has been found only recently by Swedish Inger Zeil, . Of course, to people living in North-Africa, his results have always been obvious. However, as they aroused some interest outside the Big Desert, he eventually will receive the prize without any bells.

Presently, research in sand dynamics came to an interlocutory end with a stone’s3 33footnotetext: Unfortunately, the authorship of this work is still unclear. field inequalities, , which form the basis of our approach.

The paper is organized as follows. In section two, the main aspects of the standard sand-wind duality are recapitulated. In section three a stone’s field inequalities (specialized for the heat kernel of the earth) are presented. From these fundamental considerations, section four will lead us in a natural way to a completely new and attractive theory of sandstorms. In section five follows a brief conclusion and outlook with respect to some recent aspects of sand theory will be given.

## 2 Sand-Wind Duality

This section is devoted to a brief recapitulation of the sand-wind duality of grain dynamics. Since sand is quantized (’grained’), one has to use the well-known cat-equation, which reads (in natural units)

$$(W+V_H)|\mathrm{sand}=i\frac{}{t}|\mathrm{sand}.$$
(1)
$`W`$ is the wind operator, $`V_H`$ stands for the heat field operator. We use the idempotency of sand dunes (two sand dunes thrown upon each other give again a sand dune, since the superfluous grains merely drain in the sea of sand, as introduced by the great wizard Pamdirac) and $`WV_H=V_HW=0`$ (heat and wind never occur at the same time). Now, multiplying this equation by its complex conjugate and dressing the cat with its bra leads immediately to

$$\mathrm{sand}|W|\mathrm{sand}=\mathrm{sand}|(\frac{}{t})^2|\mathrm{sand}.$$
(2)
Due to the quadratic time derivative on the right hand side of this equation, by interpreting the operators as states and the states as operators we can immediately read off the duality between sand and wind. Most interesting consequences of this phenomenon follow from the fact that a (sand theoretical) distinction between the sand-field and the wind-field is impossible. This is the basis of our modern understanding of sandstorms.

In classical sand theory, it has been believed that the genesis of sandstorms lies merely in the wind field, interacting with the grains of sand, pushing them up and around. Modern understanding of sandstorms using the sand-wind duality of grain dynamics, however, has taught that the wind is a mere manifestation of the sand field dynamics. Thus, the usual theory of sandstorms is simply that the interaction of the sand-wind-field with itself is sufficient to create a sandstorm. Indeed, it has been shown by  that the movement of the grains of sand functions as source term for the wind field, not the other way as suggested by classical theory: It’s the grains generating the wind.

So far, this is somwhat awkward but, after all, well understood. Unfortunately, combinig the principles of grain dynamics with the Relative Generality Theory of a stone, we will find that this simple, non-general model is not sufficient to a full understanding of sandstorms.

## 3 The Heatfield

In the preceding section we have deliberately assumed that the reader is familiar with some of the basic notions describing the coupling of wind to sand. However, in order to fully exploit the origin of the breakdown of sand-wind-duality in sandstorms, we shall need a more thorough understanding of the generalistic physics of sand with respect to the heat field, i.e. of a stone’s theory and its grainization. An excellent recent introduction can be found in  – in chapter 42, of course.

Let’s begin with a brief historical survey. As it is well known, the experiments of Michael’s son (an ethiopian sprinter) and a cheetah  have given a first evidence that nothing moves faster than sand and that all grains (in a sandstorm) have the same speed, at least within the experimental accuracy of that time. This led to a relatively special view of matters. However, it still remained unclear how such a theory could be combined with the interaction of wind, sand and heat in the spirit of the old ton theory . Only a few millennia later, this problem has been solved in an ingeniously simple manner by a stone.

Obviously, the sand in the desert is rather hot, in fact even hotter than the air surrounding it. Up to a stone’s era many scientists desperately tried to understand this phenomenon. A stone turned it into a new paradigm instead. In brief, his idea can be stated as follows: The heat moves the grains, but the grains are the source for the heatfield at the same time. (If we could switch off the dynamics and wait long enough then the air would be as hot as the sand.) More quantitatively, a stone’s inequalities read

$`T_g^{\mu \nu }`$ $``$ $`T_w^{\mu \nu },`$ (3)
$`G^{\mu \nu }`$ $``$ $`T_\text{?}^{\mu \nu }`$ (4)
from which one immediately derives his famous saying “E=m ceh ceh”. (The ancient letter $`=`$ has no analogue nowadays.) Note that the sand is still treated as a classical (in the sense of continuous) field (of characteristic zero) here. The coupling of the spin-2-heatfield to wind, which is spin-1, i.e. a vector-field, is described by the force $`F^\mu =h^{\mu \nu }W_\nu `$ acting on the wind quanta.

Due to its many nontrivial predictions, the success of a stone’s theory was overwhelming. Among these predictions were the delay of the time shown by a sandglass in an external heatfield and, most striking, the deflection of light in heatfields, often referred to as “Fata Morgana”. (By the way, this effect is also the origin of the common misbelief that dromedars have two hunches.)

Even more so, it turned out that the grainization of this theory was no great deal 444footnotetext: since a spin-2 grain, i.e. the heat quanta, appears without saying in the grainized theory. A typical effect that can only be understood in the framework of the full theory of grain dynamics, as it is a typical grain effect of the sand-wind-interaction, is the movement of sandpiles . We have already given a detailed account of this application in the preceding section. We have also shown therein that the resulting sand-wind duality leads to our present understanding of the genesis of sandstorms.

But we also mentioned, that this model has some drawbacks to which we shall turn our attention now.

First of all and honestly speaking, the above described mechanism for the genesis of sandstorms is only of limited theoretical value. In fact, no one has ever been able to prove that such extreme states of the sand-wind-field like sandstorms really do exist. (Actually, it is not even clear, whether this theory describes anything realistic at all.) A modern researcher can hardly take a “problem”of this kind serious, of course, but unfortunately it is not the end of the story:

Evidently, there is only a finite amount of energy in the heatfield of the earth. It came as a big surprise, when it has been discovered that certain kinds of sandstorm, so called grain ray bursts (GRB) had energies seemingly exceeding this upper bound, if a stone’s (grainized) theory was valid in this regime.

For a few hundred years this discovery aroused a lot of confusion and many scientists even started to doubt the spherical shape of the earth, until Re-Ez  remembered the forgotten singularity theorems of and , , for the heatfield . The former, , has also shown that such singularities will “radiate” large amounts of sand.

Assuming that the GRBs stem from the creation of sand due to the presence of a singularity of the heatfield he overcame all problems and eventually caused the (long overlooked) second revolution in the field of sand.

In fact, it soon became clear that there do exist many visible singularities of the heatfield. For instance, the sun is nothing but such a singularity, arising periodically and then moving at the sky. Due to the emittance of sand5 55footnotetext: Amusingly, this mechanism resulted in the ancient misbelief that the sun is made out of yellow sand., it disappears again after a certain time, depending on the state of the heatfield, i.e. the temperature. (In summer it will last longer than in winter.)

Quite recently, even the existence of a very massive heat-kernel in the center of earth has been established.

Most importantly and even more surprisingly, however, it turned out that many ordinary (non-GBR) sandstorms can be traced back to such singularities. This motivates our conjecture that all sandstorms are caused by such singularities.

Finally, the interpretation of the genesis of sandstorms as caused by singularities also resolved another puzzle:

Since sand is fermionic (it obviously respects the exclusion principle and, of course, the grains have spin $`\frac{1}{2}`$) and interacts only via heat-interactions, the system should be stable, seemingly in contradiction to the existence of sandstorms. However, sandstorms are only a local instability in the superfield, just like oases . (For this reason, the energy stored in a sandstorm does not exceed the maximal available energy.) Globally, the system is in equilibrium! This fact might be of some importance.

## 4 Sandstorms Revisited

So let’s assume that in the center of each sandstorm a singularity of the heatfield is sitting.666footnotetext: The objection that such a singularity has never been observed in experiments with non-GRBs can be ignored: People entering a sandstorm voluntarily should not be taken serious. As it is well known, in the eye of a (sand)storm the wind vanishes, and thus it is commonly believed that $`WV_H=V_HW=0`$ holds even there. But that is wrong, since $`V_H`$ is infinite in the center of a sandstorm. Accordingly, the sand-wind-duality breaks down in this region, and “that is the poodle’s core ”.

Before we can explore this fact and its consequences in detail we should describe the mechanism of the emergence of these singularities, however. This mechanism is quite similar to the mechanism discovered by Msw the Elder , when he tried to understand the so-called solar-new-tiro puzzle. (The problem why the sun emits so few sand, which is solved by a resonance mechanism. At the same time the Msw-resonance explains why the sun emerges only periodically and why it is red around the times of its appearance, respectively disappearance.)

Thus, we immediately infer that small fluctuations in the heatfield will eventually lead, by a similar resonance, to the emergence of a singularity if we apply the following well-established result :

Theorem 1. Under the above assumptions, there exists a point $`x𝐑^3`$ in which the degree of singularity of the heatfield-distribution $`\mathrm{sand}|h^{\mu \nu }|\mathrm{sand}`$ is strictly larger than $`C_{storm}`$.

Thus, the genesis of a sandstorm will be inevitable. Moreover, it is now evident, that the standard sand-wind duality is broken, viz.

$$\underset{rx}{lim}\mathrm{sand}|V_HW|\mathrm{sand}(r)0.$$
(5)
Hence, the grains in a sandstorm are driven by the wind! This is the most nontrivial and surprising result of this paper. It can and should be verified empirically.

## 5 Conclusion and Outlook

In this paper we presented in a highly suggestive and self-explaining way a new idea about the emergence of sandstorms. We have shown that the movement of the grainized sand is caused by the wind field. This result seems to be in consistency with the predictions of the original (but wrong) classical sand-theory. Note, however, that this “accordance” is not even qualitative: There are no similarities of the classical theory with our completely new and astonishing result (5).

Of course, this paper has only sketched some basic features of the new sand theory. A complete review is in preparation. There furthergoing questions about noncommutative sand, the super-connection between quicksand and worme holes and some quite interesting (although avantgardistic) new features in sand theory will be presented . Moreover, we will overcome some of the flaws this paper is suffering from, e.g. a renormalization of the infantilities due to the self-citation  will be done. This works the following way: Assuming that the problem is already removed at the first order. Then it follows from an inductive proof that no infantilities will arise at any higher order due to recursive citations. Now, curing the problem at the first loop level is rather simple. This completes the proof of renormalizability of this work.

We should also admit that this work suffers from a considerable lack of computer assistance. Unfortunately, the current generation of the abacus (used here at the Pyramid College) seems unable to perform simulations of such complexity. Not to speak of resolving a heart-tree in the fog.

Acknowledgements
I would like to thank my colleagues Mario Paschke and Volkmar Putz who made me aware of the reference . Without their support and friendship I would not exist. Furthermore I have to thank my namesake P.P. Longstockings, from whom much inspiration has been drawn in making this paper: “No one really knows…”

Figure 1 (The so-called $`\mathrm{\Omega }^{}`$-facsimile).


format: text

text:
# SLOW DYNAMICS AND LOCAL QUASI-EQUILIBRIUM. RELAXATION IN SUPERCOOLED COLLOIDAL SYSTEMS

## I Introduction

The dynamics of slow relaxation systems exhibits peculiar characteristics which make it essentially different from the dynamics of systems relaxing in much shorter time scales angell -nowak . The aging behaviour of the correlation functions struick6  and the violation of the fluctuation-dissipation theorem kurchan  and of the Stokes-Einstein relation bonn6  are among those significant features which have attracted the interest of many researchers in the last few years. Relaxation phenomena in glasses, polymers, colloids and granular matter provide innumerable situations demanding new theoretical developments whose implementation constitutes a challenge for nonequilibrium statistical mechanics theories.

The origin of that peculiar behaviour lies in the fact that during its evolution, the system rests permanently away from equilibrium. This feature would justify why results obtained when the system is close to equilibrium are not necessary valid when that condition is not fulfilled. A simple example illustrating this point is supplied from the slow dynamics of a two-level system. During the transition, the system does not thermally equilibrate with the bath and as a consequence equipartition law is not valid. It was shown in origin6  that coarsening the level of description of the path that the system follows in the configuration space going from diffusion to activation regimes leads to violation of the fluctuation-dissipation theorem. The fast decay of the diffusion modes makes the system always close to equilibrium whereas the activation process occurring at longer time scales does not guarantee validity of the theorem.

In this paper we will present an example illustrating the peculiar dynamics of slow relaxation systems. Using a Fokker-Planck description we derive a generalized Stokes-Einstein relation for supercooled colloidal liquids showing that diffusion coefficient and inverse of viscosity are not proportional to each other at low frequencies, when the system becomes activated. At higher frequencies, the particles undergo diffusion in the solvent and that relation holds.

The paper will be distributed as follows. In Section II, we will briefly review the generalization of the Onsager fluctuation theory to nonequilibrium *aging* states. Section III will be devoted to the formulation of a mean field theory describing the dynamics of high concentrated suspensions. This theory will be used in Section IV to interpret experimental results supporting the violation of the Stokes-Einstein relation in supercooled colloidal liquids. In the conclusions section we will summarize our main results.

## II Slow relaxation dynamics and local quasi-equilibrium

When the relaxation of a system is slow enough, its evolution takes place through a sequence of states characterized by a set of thermodynamical quantities in which the intensive parameters are, in general, different from those of the bath. It is then said that the system is at quasi-equilibrium.

A statistical mechanics description of those systems involves the non-stationary conditional probability density $`P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)`$, with $`\underset{¯}{\alpha }=(\alpha _1,\mathrm{},\alpha _n)`$ representing the state vector of the system, and $`\underset{¯}{\alpha }_0`$ the initial state at time $`t_0`$. The local quasi-equilibrium states are characterized by the probability density $`P_{qe}(\underset{¯}{\alpha },t)`$ satisfying

$$\underset{t\mathrm{}}{lim}P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)=P_{qe}(\underset{¯}{\alpha },t).$$
(1)
When relaxation takes place in a short time scale, $`P_{qe}(\underset{¯}{\alpha },t)`$ reduces to the local equilibrium probability density in which the intensive parameters are those of the bath.

The evolution of the conditional probability density is governed by a Fokker-Planck equation. To obtain it, we will first formulate the continuity equation

$$\frac{}{t}P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)=\frac{}{\underset{¯}{\alpha }}\left[P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)\underset{¯}{v}_{\underset{¯}{\alpha }}(\underset{¯}{\alpha }_0,\underset{¯}{\alpha };t_0,t)\right],$$
(2)
in which $`P\underset{¯}{v}_{\underset{¯}{\alpha }}`$ represents the phase space probability current with $`\underset{¯}{v}_{\underset{¯}{\alpha }}`$ the stream velocity in $`\underset{¯}{\alpha }`$-space. That velocity can be obtained from a mesoscopic version of the nonequilibrium thermodynamics for slow relaxation systems we have proposed Slow6 .

The existence of a local quasi-equilibrium state in $`\underset{¯}{\alpha }`$-space enables one to formulate the Gibbs equation

$$T(t)\delta s_{qe}(t)=\delta e_{qe}(t)\mu _{qe}(t)\delta P(t)𝑑\underset{¯}{\alpha },$$
(3)
where $`s_{qe}(t)`$ is the entropy and $`e_{qe}(t)`$ the mean internal energy per unit mass at quasi-equilibrium . We will assume that the temperature $`T(t)`$ of the system at local quasi-equilibrium is only a function of time whereas the chemical potential $`\mu _{qe}(t)`$ can be a function of $`\underset{¯}{\alpha }`$.

The entropy of the system can be expressed in terms of the probability density by means of the Gibbs entropy postulate kampenSuper  in the form

$$s(t)=\frac{k}{m}P(t_0|t)\mathrm{ln}\frac{P(t_0|t)}{P_{qe}(t)}d\underset{¯}{\alpha }+s_{qe}(t),$$
(4)
where $`k`$ is the Boltzmann constant and $`m`$ the molecular mass.

The form of the current $`P\underset{¯}{v}_{\underset{¯}{\alpha }}`$ can be inferred from the entropy production by relating it, as done in nonequilibrium thermodynamics degroot6 , with the corresponding thermodynamic force. The entropy production is obtained by combining the rate of change of $`s_{qe}(t)`$ obtained from Eq. (3) with the time derivative of Eq. (4). After using Eq. (2) and integrating by parts assuming that $`P(t_0|t)\underset{¯}{v}_{\underset{¯}{\alpha }}(t_0,t)`$ vanishes at the boundary, one may identify the non-equilibrium chemical potential $`\mu (t_0;t)=\mu _{qe}(t)+\frac{kT(t)}{m}\mathrm{ln}\frac{P(t_0|t)}{P_{qe}(t)}`$ and the corresponding thermodynamic force $`\frac{\mu (t)}{\underset{¯}{\alpha }}`$, which is conjugated to the probability current $`P\underset{¯}{v}_{\underset{¯}{\alpha }}`$. Following the scheme of nonequilibrium thermodynamics one simply obtains

$$\underset{¯}{v}_{\underset{¯}{\alpha }}(t_0;t)=\underset{¯}{\underset{¯}{B}}(t_0;t)\left[\underset{¯}{X}(\underset{¯}{\alpha },t)+\frac{kT(t)}{m}\frac{}{\underset{¯}{\alpha }}\mathrm{ln}P(t_0|t)\right],$$
(5)
where we have defined the generalized force degroot6

$$\underset{¯}{X}(\underset{¯}{\alpha },t)\frac{}{\underset{¯}{\alpha }}\left[\frac{kT(t)}{m}\mathrm{ln}P_{qe}(t)\right].$$
(6)
The time dependent transport coefficients $`\underset{¯}{\underset{¯}{B}}(t_0;t)`$ (related to the Onsager coefficients) incorporate memory effects through its dependence on time Fp-nonmarkovian6 , oliveiraEURO , and the local quasi-equilibrium probability density is given by

$$P_{qe}(t)=P_ee^{\frac{m}{kT(t)}[\mu _{qe}e_{qe}]},$$
(7)
where $`P_e`$ is the equilibrium distribution function. Eq. (7) can be obtained by using the statistical definition of the entropy in the Gibbs equation (3), Slow6 . Now, by substituting Eq. (5) into (2) the resulting Fokker-Planck equation is

$$\frac{P(t_0|t)}{t}=\frac{}{\underset{¯}{\alpha }}\underset{¯}{\underset{¯}{B}}(t_0;t)\left[\underset{¯}{X}(\underset{¯}{\alpha },t)P(t_0|t)+\frac{kT(t)}{m}\frac{P(t_0|t)}{\underset{¯}{\alpha }}\right],$$
(8)
Using the Fokker-Planck equation (8), we can calculate the evolution equation for the equal-time correlation function of $`\underset{¯}{\alpha }`$-variables which for sufficiently long times leads to the following expression for the temperature of the system

$$T(t)\frac{m}{k}\underset{¯}{X}(t)\underset{¯}{\alpha }(t)_{qe},$$
(9)
with $`\mathrm{}_{qe}`$ denoting average at local quasi-equilibrium. This result reduces to the one corresponding to fast relaxation processes degroot-pp88  in which the average is performed at local equilibrium.

In Ref. Slow6 , we have shown that the temperatures of the system and of the bath $`T_B(t)`$ are related by $`T(t)=AT_B(t)`$, expressing lack of thermal equilibration. Equilibrium is reached when $`A=1`$ in which case, for linear thermodynamic forces, one recovers the expression for equipartition law. The presence of $`T(t)`$ in the last term of Eq. (8), gives rise to a modified version of the fluctuation-dissipation theorem.

## III Fokker-Planck description of concentrated colloidal suspensions

In this section, we will study diffusion in concentrated colloidal suspensions by means of an effective medium approach elaborated on the grounds of the theory we have introduced in the previous section. We will analyze the motion of a test colloidal particle through the suspension by using the Fokker-Planck equation describing the evolution of the two-time probability density $`P(\underset{¯}{\alpha }_0,t_0|\underset{¯}{\alpha },t)`$. In this case, the phase space vector is simply $`\underset{¯}{\alpha }=(\stackrel{}{u},\stackrel{}{r})`$, where $`\stackrel{}{u}`$ and $`\stackrel{}{r}`$ represent the velocity and position of the test particle.

To derive the Fokker-Planck equation, we will first formulate the Gibbs equation

$$T\rho \delta s=\rho \delta e\rho ^1p\delta \rho m\mu \delta P𝑑\stackrel{}{u},$$
(10)
which incorporates the effects of interactions among particles through the term containing the excess of osmotic pressure $`p`$. Here $`e`$ is the internal energy, $`m`$ the mass of the particle and $`\rho =mP𝑑\stackrel{}{u}`$ the mass density. In this case, the quasi-equilibrium probability density is given by Slow6

$$P_{qe}=e^{\frac{m}{2kT}\left[\mu _{qe}\frac{1}{2}u^2\rho ^1p\right]}.$$
(11)
Here, the interaction of the test particle with the other colloidal particles is represented by means of the term $`\rho ^1p`$ which can be expressed in terms of a virial expansion in $`\rho `$, whereas $`\mu _{qe}`$ constitutes the ideal chemical potential. Introducing the fugacity $`z=Pa`$, with the activity coefficient given by $`ae^{\frac{m}{kT}\rho ^1p}`$, Eq. (4) can be expressed as

$$s=\frac{k_B}{\rho }P\mathrm{ln}\frac{z}{z_{qe}}d\stackrel{}{u}+s_{qe},$$
(12)
where $`z_{qe}=e^{\frac{m}{kT}[\mu _{qe}\frac{1}{2}u^2]}`$ is the fugacity at quasi-equilibrium. Notice that in the limit of infinite dilution $`\rho ^1p0`$, thus implying $`a=1`$, which corresponds to the ideal case.

The existence of interactions among the test and the other colloidal particles, gives rise to a force term in the continuity equation for $`P`$. In particular, according to Eq. (11), the force entering into this term is given by $`(\rho ^1p)`$ (see for example Ref. mayorgas ). Now, by taking into account the definition of the activity coefficient $`a`$, the continuity equation for $`P`$ becomes

$$\frac{P}{t}+P\stackrel{}{u}\left[\frac{kT(t)}{m}\mathrm{ln}a\right]\frac{P}{\stackrel{}{u}}=\frac{}{\stackrel{}{u}}P\stackrel{}{v}_\stackrel{}{u},$$
(13)
where $`P\stackrel{}{v}_\stackrel{}{u}`$ is the probability current, with $`\stackrel{}{v}_\stackrel{}{u}`$ the stream velocity in $`\stackrel{}{u}`$-space.

In accordance with the general formalism of nonequilibrium thermodynamics, the expression of the probability current $`P\stackrel{}{v}_\stackrel{}{u}`$ can be obtained from the entropy production which only contains dissipative terms and is simply given by

$$\sigma =\frac{m}{T}P\stackrel{}{v}_\stackrel{}{u}\frac{\mu }{\stackrel{}{u}}𝑑\stackrel{}{u}.$$
(14)
Eq. (14) has been obtained by taking the time derivative of Eq. (12), using Eqs. (13), (11) and (10), and the balance equation for the energy $`e(t)`$, neglecting viscous dissipation degroot6 . In order to obtain Eq. (14), we have also identified the nonequilibrium chemical potential

$$\mu (t_0;t)=\mu _{qe}(t)+\frac{kT(t)}{m}\mathrm{ln}\frac{z(t_0;t)}{z_{qe}(t)}.$$
(15)
By taking into account the definition (6), the resulting generalized Fokker-Planck equation for the two-time probability density is

$$\frac{P}{t}+\stackrel{}{u}P\left[\frac{kT(t)}{m}\mathrm{ln}a\right]\frac{P}{\stackrel{}{u}}=\frac{}{\stackrel{}{u}}\beta (t)\left[\stackrel{}{u}P+\frac{kT}{m}\frac{P}{\stackrel{}{u}}\right],$$
(16)
where $`\beta `$ characterizes the dissipation of the kinetic energy of the test particle. This coefficient is, in general, function of time, position and density; its dependence on time is a consequence of the non-Markovian nature of the stochastic process Fp-nonmarkovian6 , oliveiraEURO , tokuyama .

From Eq. (16), we will now derive the evolution equation for the velocity field by taking the time derivative of the momentum $`\rho \stackrel{}{v}=m\stackrel{}{u}P𝑑\stackrel{}{u}`$ and using Eq. (16). After integrating by parts one obtains

$$\rho \frac{d\stackrel{}{v}}{dt}=\beta \rho \stackrel{}{v}\frac{kT(t)}{m}\left[1+\frac{\mathrm{ln}a}{\mathrm{ln}\rho }\right]\rho ,$$
(17)
where we have used the fact that, in the absence of an externally imposed flow, the second moment $`\stackrel{}{\stackrel{}{\mathrm{P}}}=m(\stackrel{}{u}\stackrel{}{v})(\stackrel{}{u}\stackrel{}{v})P𝑑\stackrel{}{u}`$ can be approximated by $`\stackrel{}{\stackrel{}{\mathrm{P}}}=\frac{kT}{m}\rho \stackrel{}{\stackrel{}{1}}`$, nosotrosPRE6 .

At times such that $`t\beta ^1`$, the particle enters the diffusion regime in which one may neglect the time derivative in Eq. (17), obtaining then the constitutive relation for the mass diffusion current $`\rho \stackrel{}{v}`$ which substituted into the mass continuity equation $`\frac{\rho }{t}=(\rho \stackrel{}{v})`$, yields the Smoluchowski equation

$$\frac{\rho }{t}=\left(D\rho \right).$$
(18)
From Eq. (17) we may identify the diffusion coefficient

$$D=\frac{kT(t)}{m\beta }\left[1+\frac{\mathrm{ln}a}{\mathrm{ln}\rho }\right].$$
(19)
The form of Eq. (19) coincides with the expression given in Ref. degroot6 . Notice, however, that it contains the temperature of the system at quasi-equilibrium instead of the temperature of the bath.

### III.1 The generalized Stokes-Einstein relation in supercooled colloidal systems

It has been shown in bonn6  that in supercooled colloidal liquids at short frequencies the Stokes-Einstein relation is not fulfilled. Measurements of the viscosity and the diffusion coefficients obtained in Refs. mortensen6 , mason6 , indicate the existence of a frequency domain in which those quantities are not proportional (see figure 1). It is our purpose to offer here an explanation for those experimental results.

In mortensen6  it was shown that the colloidal system can be assumed isotropic and homogeneous at the considered volume fractions, ranging from the freezing to the glass transition values $`\varphi _g`$. Consequently, the transport coefficients can be assumed to be independent of the wave vector $`\stackrel{}{k}`$. We will first introduce the effective viscosity

$$\eta (\omega )=\frac{kT}{6\pi R}D^1,$$
(20)
where $`R`$ is the radius of the particle.

The dependence of the viscosity on the frequency can be inferred through the expression of the penetration length $`\lambda `$ of the diffusive modes. In the dilute case, $`\lambda `$ is given through $`\lambda ^1=\frac{1}{R}\left(\tau \omega \right)^{\frac{1}{2}}`$, with $`\tau =\frac{R^2}{6D_0}`$ the characteristic diffusion time boon-yip  and $`D_0=\frac{kT_B}{6\pi R\eta _0}`$ the diffusivity in the pure solvent, with $`\eta _0`$ its viscosity. At higher concentrations, it can be assumed to be

$$\lambda ^1\frac{1}{\mathrm{\Gamma }R}\left(\tau \omega \right)^\gamma ,$$
(21)
where the exponent $`\gamma `$ and the scaling function $`\mathrm{\Gamma }`$ may, in general, depend on the volume fraction. The quantity $`\mathrm{\Gamma }R`$ can be interpreted as an effective radius.

The penetration length gives an estimation for the average size of a cage formed by the surrounding particles cohen6 . At high frequencies, the test particle performs a free Brownian motion inside the cage and the viscosity is that of the pure solvent. At short frequencies, collisions of the test particle with the other particles modify the viscosity. Thus, we will assume the following expression for the effective viscosity

$$\eta (\omega )=\eta _0\left(1+R\lambda ^1\right),$$
(22)
valid up to first order in $`R\lambda ^1`$ reflecting the fact that cage diffusion is important when $`R\lambda `$, cohen-cage .

By taking into account the general relation between system and bath temperatures obtained from our formalism, $`T=AT_B`$, Slow6 , Eq. (20) can be written as

$$D\eta =\frac{kT_B}{6\pi R}A,$$
(23)
which constitutes a generalization of the Stokes-Einstein relation to systems exhibiting slow relaxation dynamics.

From Eq. (9) and from its definition, the quantity $`A`$ is essentially proportional to the velocity moments of the particle. In particular, when $`\underset{¯}{X}(t)`$ is a linear function, $`A`$ is proportional to the velocity correlation function. It is then plausible to assume the following expansion ngai , tokuyamaslow

$$A\left[1+b\left(\tau \omega \right)^ϵ\right],$$
(24)
where $`b`$ and the exponent $`ϵ`$ may, in general, be functions of the volume fraction $`\varphi `$. Comparison of Eq. (24) with the corresponding relation obtained in bonn6  by fitting the experimental data yields $`\mathrm{ln}bϵ\left(40\frac{\varphi }{\varphi _g}37\right)+2.659`$, and $`ϵ=0.77`$.

To contrast the theoretical expression of the effective viscosity (22) with the experimental results reported in bonn6 , we will take the logarithm of Eq. (22) and expand the result in terms of the variable $`\mathrm{log}\tau \omega `$, around $`\tau \omega =1`$. Up to first order we obtain

$$\mathrm{log}\eta \mathrm{log}\left[\eta _0\left(1+\frac{1}{\mathrm{\Gamma }}\right)\right]+\frac{\gamma }{1+\mathrm{\Gamma }}\mathrm{log}\tau \omega .$$
(25)
A direct comparison of Eq. (25) with the fitting of experimental data given by the expression $`\mathrm{log}\eta c_0+s(\varphi )\mathrm{log}\tau \omega `$, with $`s(\varphi )6.3913\varphi `$ and $`c_0`$ a constant, gives the dependence of $`\gamma `$ and $`\mathrm{\Gamma }`$ on the volume fraction.

The expression of the diffusion coefficient at the frequencies we are considering can now be derived by using Eq. (23). In figure 2, we represent the diffusion coefficient given through Eq. (23) as a function of the reduced frequency. At higher concentrations ($`\varphi =0.56,\mathrm{\hspace{0.17em}0.55}`$), a closer agreement with the experimental values of $`D`$ given in bonn6  has been obtained with $`ϵ=0.77`$ (solid lines). At lower concentrations ($`\varphi =0.53,\mathrm{\hspace{0.17em}0.52}`$), a closer agreement with experiments has been obtained with $`ϵ=0.67`$ (dotted lines). This dependence of the fitting on the value of the exponent $`ϵ`$ suggests that it is also a function of the volume fraction.

## IV Conclusions

In this paper, we have proposed a general formalism to analyze the dynamics of systems which relax in long time scales of the order of the observation time or even longer. The fact that the system always evolves through transient states never reaching equilibration with the bath is the origin of a peculiar behaviour different from that of systems whose propagating modes relax rapidly ending up in a quiescent state. The dynamics is non-Markovian and the correlation functions exhibit aging effects.

We have formulated a generalization of Onsager’s theory to nonequilibrium aging states to derive a Fokker-Planck equation which captures the main characteristics of the dynamics such as non-stationarity through the two-time probability density, non-Markovianity through the time-dependence of the transport coefficients and lack of thermal equilibration between system and bath. This equation is the main result of our analysis from which the behaviour of the correlation functions follows.

We have applied the theory to study the relaxation in supercooled colloidal liquids for which violation of the Stokes-Einstein relation at low frequencies has recently been found experimentally. The relation between diffusion coefficient and viscosity we have obtained explains those experiments. The key ingredients in this interpretation are the existence of a quasi-equilibrium state with a temperature different from that of the bath which is proportional to the velocity moments and the fact that the viscosity depends on the penetration length of the diffusion modes resulting from the underlying activated process taking place at sufficiently long times.

The formalism proposed thus offers an interesting theoretical framework for the study of glassy behaviour in colloidal liquids.

## V Acknowledgments

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

## VI Acknowledgments


text:
# One dimensional quantum walk with unitary noise

## I Introduction

Random walk models describe a great variety of diffusion phenomena in physical systems. Such phenomena include the diffusion of particles in a fluid (Brownian motion), the motion of vacancies in a crystal and of atoms on a crystalline surface. Related models are also used to describe the spatial structure of systems such as polymer chains.

A random walker on a lattice hops at each time step from its present site to one of its nearest neighbors. The hopping direction is picked randomly, with, for example, the probability to hop to each of $`c`$ nearest neighbors given by $`1/c`$. Random walk models on discrete lattices as well as in the continuum have been studied extensively.

Consider a random walker on a one dimensional lattice. Denote the probability to find it at site $`n`$ at time $`t`$ by $`P_t(n)`$. The time evolution of $`P_t(n)`$ is described by the recursion equation

$$P_{t+1}(n)=\frac{1}{2}[P_t(n1)+P_t(n+1)],n=0,\pm 1,\pm 2,\mathrm{}$$
(1)
For a walker starting from the origin at $`t=0`$, the probabilities, $`P_t(n)`$, of the walker to be at site $`n`$ at time $`t`$, are given by the components of the Pascal triangle, namely

$$P_t(n)=\{\begin{array}{cc}\frac{t!}{2^t\left(\frac{tn}{2}\right)!\left(\frac{t+n}{2}\right)!}:\hfill & n=t,t+2,\mathrm{},t\hfill \\ 0:\hfill & n=t+1,t+3,\mathrm{},t1,\hfill \end{array}$$
(2)
and $`P_t(n)=0`$ for $`|n|>t`$. Note that at even times only the even sites can be occupied, while at odd times only the odd sites can be occupied. A continuum description of the random walk can be obtained when the lattice constant $`\mathrm{\Delta }n0`$ and the time step $`\mathrm{\Delta }t0`$, such that $`D=\mathrm{\Delta }n^2/(2\mathrm{\Delta }t)`$ converges to a finite value. The parameter $`D`$ is called the diffusion coefficient and it quantifies the rate in which the random walker moves. In this limit the random walk can be described by the diffusion equation

$$\frac{dP_t(n)}{dt}D\frac{d^2P_t(n)}{dn^2}=0,$$
(3)
where $`n`$ and $`t`$ are continuous variables. The solution of this equation, approached by Eq. (2) in this limit, is the Gaussian probability distribution

$$P_t(n)=\frac{1}{\sqrt{4\pi Dt}}e^{\frac{n^2}{4Dt}}$$
(4)
with $`D=1/2`$. The standard deviation of this distribution, $`\sigma (t)=\sqrt{2Dt}`$, thus takes the form $`\sigma (t)=\sqrt{t}`$.

Recently, quantum analogues of the classical random walk model have been studied on the one dimensional lattice AharonovY93 ; Meyer96 ; Ambainis01 ; Mackay01 ; Travaglione2002  as well as on more general graphs AharonovD01 ; Childs01 . Consider the quantum walk on the one dimensional lattice. Unlike the classical random walker that occupies a single site at a time, the quantum walker can be in an extended state. This state is a superposition of all the basis states $`|n`$, $`n=0,\pm 1,\pm 2,\mathrm{}`$, where $`|n`$ is the state in which the walker is located at site $`n`$. In addition to the spatial degree of freedom the walker has a chirality qubit, which can be in a superposition of the states $`|R`$ and $`|L`$, and determines the direction of the next hop. Each move of the walker consists of two unitary operations. The first one is a unitary transform, taken here to be the Hadamard transform, on the chirality qubit (expressed in the standard basis). The second one is the actual move, in which a walker with chirality $`|R`$ moves to the right and a walker with chirality $`|L`$ moves to the left. This definition resembles the classical walk in the sense that a walker in a basis state $`|n|s`$ ($`s=R`$ or $`L`$) at time $`t`$, will have equal probabilities to be found in the $`n+1`$ and $`n1`$ sites if a measurement is performed at time $`t+1`$. However, the coherent motion of the state gives rise to constructive and destructive interference that strongly modifies the emerging probability distribution. As a result, a quantum walk starting from the origin at time $`t=0`$ moves away faster than the classical random walk. The standard deviation $`\sigma (t)`$ of the probability distribution $`P_t(n)`$ for the quantum walk increases with time according to $`\sigma (t)t`$, while for the classical random walk $`\sigma (t)\sqrt{t}`$.

One of the main difficulties in the experimental implementation of quantum algorithms is the sensitivity of the quantum systems to noise Bernstein97 ; Preskill ; Nielsen00  and decoherence Zurek01 . The problem of decoherence can be tackled by using quantum error correction methods Shor95 ; Shor96a ; Steane96 ; Steane96a ; Knill97  as well as decoherence free sub-spaces Zanardi97 ; Lidar98 ; Bacon99 ; Kempe01 . However, these methods require significant overhead, making the quantum circuits more complicated. It is thus useful to examine the effect of decoherence on various quantum algorithms, implemented on unprotected quantum circuits Azuma2002 . The quantum walk model is particularly suitable for this task since its behavior is qualitatively different from the classical random walk. The effect of decoherence on the quantum walk has been studied recently Brun2002 ; Kendon2002 ; Brun2003 ; Kendon2003 ; Kendon2003b . It was found that the quantum walk is highly sensitive to decoherence, namely, even weak decoherence gives rise to classical-like behavior in the long time limit. Decoherence is a fundamental problem because it is a result of the unavoidable interaction of the system with the environment, that makes the two of them entangled. This interaction can be described by a variety of noise models, bringing the system into a mixed state, by essentially tracing out the environment degrees of freedom.

A related problem is the effect of unitary noise, which appears as a result of fluctuations and drifts in the generating Hamiltonian of any given unitary operation Shapira2003 . This noise tends to reduce the performance of quantum devices, but does not cause an entanglement with the environment, namely the system remains in a pure state.

In this paper we analyze the effect of unitary noise on the discrete one-dimensional quantum walk model, using computer simulations. We find that even a tiny noise level will eventually induce a crossover from the quantum walk into a classical-like behavior, characterized by $`\sigma (t)\sqrt{t}`$. The crossover time $`T`$ is calculated numerically. It is found that $`T\alpha ^2`$, where $`\alpha `$ is the standard deviation of the noise distribution.

The paper is organized as follows. The quantum walk model on the one-dimensional lattice is described in Sec. II. The noisy quantum walk is introduced in Sec. III. Simulations and results are presented in Sec. IV, followed by a summary in Sec. V.

## II The quantum walk model

Consider a quantum walker on a one-dimensional lattice. The state of the system at time $`t`$ is

$$|\psi (t)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\underset{s=L}{\overset{R}{}}a_{n,s}(t)|n|s,$$
(5)
where the amplitudes $`a_{n,s}`$ are complex numbers that satisfy $`_{n,s}|a_{n,s}|^2=1`$. The state vector $`|n`$, where $`n=0,\pm 1,\pm 2,\mathrm{}`$, represents the position of the walker on the lattice. The state vector $`|s=|R`$ or $`|L`$ represents the chirality degree of freedom. The chirality consists of a single qubit, and its state determines the coefficients $`a_{n,s}(t+1)`$. Operations in the chirality space in the quantum model, replace the randomized decision that appears in the classical walk model and determines the hopping direction. When measuring the walker position at a certain time $`t`$, the probability to find the walker at site $`n`$ is given by:

$$P_t(n)=|a_{n,R}(t)|^2+|a_{n,L}(t)|^2.$$
(6)
The time evolution of the quantum walk is expressed by

$$|\psi (t+1)=\widehat{Q}_0|\psi (t),$$
(7)
where $`|\psi (t)`$ is the quantum state at time $`t`$. The quantum walker’s step is defined by

$$\widehat{Q}_0=\widehat{T}\widehat{U}_0,$$
(8)
where $`\widehat{U}_0`$ is a unitary operator that applies only in the chirality space and takes the role of the “coin” in the classical random walk. Here we focus the case of the Hadamard walk, in which $`\widehat{U}_0=\widehat{I}\widehat{w}_0`$, where

$$\widehat{w}_0=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)$$
(9)
is the Hadamard operator. The states of the chirality qubit are expressed in the standard basis

$$|R=\left(\begin{array}{c}1\\ 0\end{array}\right);|L=\left(\begin{array}{c}0\\ 1\end{array}\right),$$
(10)
and $`\widehat{I}`$ is the spatial identity operator. The computational basis state vectors are transformed by the Hadamard operator according to

$`\widehat{U}_0|n|R`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}|n(|R+|L)`$
$`\widehat{U}_0|n|L`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}|n(|R|L).`$ (11)
The translation operator $`\widehat{T}`$ then performs the walker’s move according to the chirality state such that:

$`\widehat{T}|n|R`$ $`=`$ $`|n+1|R`$
$`\widehat{T}|n|L`$ $`=`$ $`|n1|L.`$ (12)
These operators yield the following recursion equations for the amplitudes of the quantum state

$`a_{n,R}(t+1)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}[a_{n1,R}(t)+a_{n1,L}(t)]`$
$`a_{n,L}(t+1)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}[a_{n+1,R}(t)a_{n+1,L}(t)].`$ (13)
For a given initial state, the recursion equations provide the probability of finding the walker at any site $`n`$ at time $`t`$. Here we focus on the case that at $`t=0`$ the walker is located at the origin, namely, $`a_{n,s}(0)=0`$ for all $`n0`$ (and $`s=L,R`$).

## III The noisy quantum walk

The quantum walk, like other quantum computing systems, may be affected by unitary noise. This noise is due to fluctuations and drifts in the parameters of the quantum Hamiltonian of the system. The perturbed Hamiltonian is still Hermitian and therefore generates time evolution operators that are unitary, but now include a stochastic part. Formally, one can write a noisy unitary operator as

$$\widehat{U}=\widehat{U}_0e^{i\widehat{A}},$$
(14)
where $`\widehat{A}`$ is a stochastic Hermitian operator determined by the perturbation and $`\widehat{U}_0`$ is the time evolution operator of the original quantum process without the perturbation.

The Hadamard walk model consists of Hadamard gates that act in the chirality space at each time step of the walker. Here we consider the effect of unitary noise in the Hadamard operator on the quantum walker. To this end we describe a move of the noisy walker at time $`t`$ as:

$$|\psi (t+1)=\widehat{Q}(t)|\psi (t).$$
(15)
where

$$\widehat{Q}(t)=\widehat{T}(\widehat{I}\widehat{w}_0e^{i\widehat{a}(t)}).$$
(16)
The operator $`\widehat{a}(t)`$ is a stochastic and Hermitian operator that applies in the single-qubit chirality space at time $`t`$. It can be expanded in the basis of the Pauli operators such that

$$\widehat{a}(t)=\alpha _1(t)\widehat{\sigma }_1+\alpha _2(t)\widehat{\sigma }_2+\alpha _3(t)\widehat{\sigma }_3,$$
(17)
where $`\alpha _k(t)`$, $`k=1,2,3`$ are real stochastic variables, and the $`\sigma `$’s are the Pauli operators, represented in the spin basis as:

$$\widehat{\sigma }_1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right);\widehat{\sigma }_2=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right);\widehat{\sigma }_3=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$
(18)
Although, in principle, the $`2\times 2`$ identity operator should also be included in the expansion (17), it is omitted because it can only change the overall phase in the quantum walk process, and does not have any effect on the walker’s measurement probabilities.

We will now analyze the properties of the noisy Hadamard walk according to the characteristics of the real stochastic variables $`\alpha _k(t)`$, $`k=1,2`$ and $`3`$. In the analysis we assume that there is no correlation between different noise components $`k`$ and $`k^{}k`$ as well as between different walker steps at times $`t`$ and $`t^{}t`$, namely

$$\alpha _k(t)\alpha _k^{}(t^{})=\delta _{k,k^{}}\delta _{t,t^{}}\alpha ^2,$$
(19)
where $`\delta _{k,k^{}}`$ is the Kronecker delta function. Furthermore, we focus on the isotropic case in which $`\alpha _1`$, $`\alpha _2`$ and $`\alpha _3`$ are taken from the same distribution $`p(\alpha )`$, which is unbiased, namely

$$\alpha _k(t)=0,k=1,2,3,$$
(20)
with standard deviation $`\alpha `$.

To analyze the time evolution of the Hadamard walk in the presence of unbiased and isotropic unitary noise, we have performed computer simulations of the model. In these simulations the Hadamard gate is replaced by a noisy one, so the walker step at time $`t`$ is given by Eq. (16). The Pauli coefficients $`\alpha _1(t)`$, $`\alpha _2(t)`$ and $`\alpha _3(t)`$ are taken from a Gaussian distribution with zero average and a standard deviation $`\alpha `$. The magnitude of the isotropic noise is characterized by the values of $`\alpha `$. To obtain proper statistics the results were averaged over a sufficient number of runs. For each value of $`\alpha `$, we have applied $`10,000`$ steps of the noisy Hadamard walk for at least $`200`$ runs in the case of weak noise ($`\alpha <0.07`$) and more than $`4,000`$ runs in case of strong noise ($`\alpha 0.07`$). The distribution $`P_t(n)`$ of the walker’s position was then examined and its moments are calculated.

The probability to measure the noisy Hadamard walker in site $`n`$ at time $`t`$ is given by:

$$P_t(n)_\alpha =|a_{n,R}(t)|^2_\alpha +|a_{n,L}(t)|^2_\alpha ,$$
(21)
where $`a_{n,R}(t)`$ and $`a_{n,L}(t)`$ are amplitudes of the walker’s state $`|\psi (t)`$ \[see Eq. (5)\]. Here $`_\alpha `$ denotes the averaging over the noise, taken from a Gaussian distribution with a zero average and standard deviation $`\alpha `$.

Here we consider only initial conditions for which at $`t=0`$, the walker is located at the origin, namely $`a_{n,s}(0)=0`$ for all $`n0`$. For these initial states it is guaranteed that $`P_t(n)=0`$ for all $`|n|>t`$. Therefore the first and the second moments of the spatial distribution at time $`t`$ (which are both averaged over the noise) are given by:

$$\overline{n(t)}_\alpha =\underset{n=t}{\overset{t}{}}nP_t(n)_\alpha $$
(22)
and

$$\overline{n^2(t)}_\alpha =\underset{n=t}{\overset{t}{}}n^2P_t(n)_\alpha .$$
(23)
The standard deviation $`\sigma _\alpha (t)`$ of the spatial distribution of the noisy Hadamard walk at time $`t`$ is given by

$$\sigma _\alpha ^2(t)\overline{n^2(t)}_\alpha \overline{n(t)}_{\alpha }^{}{}_{}{}^{2}.$$
(24)
## IV Simulations and results

To analyze the effect of noise on the Hadamard walk we have performed direct computer simulations of the system, using Eq. (16). The initial condition used in the simulations is

$$|\psi (0)=\frac{1}{\sqrt{2}}(|0|R+i|0|L).$$
(25)
For this state the probability distribution $`P_t(n)`$ of the Hadamard walk without noise turns out to be symmetric around the origin Travaglione2002 , namely $`P_t(n)=P_t(n)`$. In Fig. 1 we present the probability distribution $`P_t(n)`$, $`n=t,t+2,\mathrm{},t`$, for $`t=250`$ (a), $`1000`$ (b), and $`10,000`$ (c) steps. The distribution exhibits two main peaks, on opposite sides of the origin, that move away from each other as time evolves. In the central part it forms a plateau, in which $`P_t(n)1/t`$. It also exhibits wild oscillations under the envelope, that decays towards the origin and rises towards the two peaks. This distribution is clearly very different from the Gaussian distribution obtained for the classical random walk. The formation of $`P_t(n)`$ can be understood as a result of constructive interference in the front of the propagating wave function, while the center is dominated by destructive interference.

When unitary noise is added, it perturbs the structure produced by constructive and destructive interference. The distribution $`P_t(n)`$ for a weak noise level ($`\alpha =0.025`$) is shown in Fig. 2 for $`t=250`$ (a), $`1000`$ (b), and $`10,000`$ (c) steps. The effects of the competition between the behavior of the original Hadamard walk and the unitary noise are clearly shown. First, alongside the oscillating wings on the edges, an incoherent component is generated around the center \[Fig. 2(a)\]. Gradually, the central peak increases absorbing weight from the two oscillating wings \[Fig. 2(b)\], which later disappear completely, \[Fig. 2(c)\]. However, the probability distribution does not converge to a Gaussian even at much longer times, since at any time step the quantum effects continue to be dominant. The incoherent component becomes more dominant as time evolves, due to the accumulation of noise effects. However, the coherent features of the quantum walk appear in every time step and the probability distribution does not seem to converge to a Gaussian.

In Fig. 3 we present the probability distribution $`P_t(n)`$ for the case of strong noise ($`\alpha =0.8`$) at $`t=1000`$. At this noise level the quantum effects are completely suppressed by the noise. The oscillations are smoothed and the probability distribution converges towards a Gaussian distribution with a zero mean, that expands with time. It thus approaches the behavior of the classical random walk.

In Fig. 4 we present the time dependence of the standard deviation $`\sigma _\alpha (t)`$ of $`P_t(n)`$ for the noisy Hadamard walk with symmetric, unbiased noise. The standard deviation $`\sigma _\alpha (t)`$, given by Eq. (24), is shown for different noise levels $`\alpha `$. The top curve is for the noiseless quantum walk, the next four curves, from top to bottom, are for $`\alpha =0.025`$, $`0.05`$, $`0.1`$ and $`0.2`$, and the last curve at the bottom is for the classical random walk. The initial state used in these simulations is

$$|\psi (0)=|0|R.$$
(26)
We observe that as the noise level $`\alpha `$ increases, the standard deviation curves for the noisy quantum walk move down from the top curve of the noiseless quantum walk and approach the bottom curve describing the classical random walk.

The standard deviation $`\sigma _0(t)`$ for the noiseless quantum walk as a function of time $`t`$ takes the form Ambainis01

$$\sigma _0(t)=qt.$$
(27)
From our numerical simulations we obtain that the coefficient $`q=0.4505\pm 0.0005`$, which is in perfect agreement with the analytical results of Ref. Ambainis01 . The standard deviation for the classical random walk is given by

$$\sigma _{\text{classical}}(t)=\sqrt{t}.$$
(28)
As long as the noise level is not too high, the noisy Hadamard walk shows the typical behavior of a quantum walk for short times. However, beyond some crossover time the motion of the noisy quantum walk acquires diffusive features in the sense that $`\sigma _\alpha (t)`$ behaves like the classical random walk. However, the shape of $`P_t(n)`$ does not approach a Gaussian and maintains a far reaching tail.

In order to evaluate the crossover time, and examine its dependence on $`\alpha `$, we will use the short-time and long-time limits of $`\sigma _\alpha (t)`$. For short times $`\sigma _\alpha (t)`$ increases linearly with $`t`$, according to Eq. (27). Then, gradually the accumulation of the noise modifies the curve until in the long time limit it forms a square-root shape. The long-time tail of $`\sigma _\alpha (t)`$ is well fitted by $`\sigma _\alpha (t)=K(\alpha )\sqrt{t}+C(\alpha )`$, and the functions $`K(\alpha )`$ and $`C(\alpha )`$ can be obtained. It turns out that for long enough times $`C(\alpha )`$ can be neglected, giving rise to

$$\sigma _\alpha (t)=K(\alpha )\sqrt{t}.$$
(29)
The function $`K(\alpha )`$ is then obtained by linear regression of the points $`[\sqrt{t},\sigma _\alpha (t)]`$, for each value of $`\alpha `$. The function $`K(\alpha )`$ for a relevant range of $`\alpha `$ values is shown in Fig. 5. It turns out to be a monotonically decreasing function, that in the limit of small $`\alpha `$ takes the form $`K(\alpha )1/\alpha `$. In the limit of very large $`\alpha `$, $`K(\alpha )1`$, and the noisy quantum walk coincides with the classical random walk, for which $`K=1`$ (dashed line). The monotonically decreasing behavior of $`K(\alpha )`$ means that as the noise level increases, the broadening of the distribution $`P_t(n)`$ in the long time limit becomes slower. This is due to an earlier suppression of the quantum coherence, that broadens the distribution much faster than the noise.

We will now consider the crossover time $`T_2(\alpha )`$, which is based on the behavior of $`\sigma _\alpha (t)`$. It is defined as the time at which the long time asymptote of the noisy Hadamard walk \[given by Eq. (29)\] intersects with the short time asymptote, namely the ideal noiseless linear curve \[Eq. (27)\]. This crossover time is given by

$$K(\alpha )\sqrt{T_2(\alpha )}=qT_2(\alpha ),$$
(30)
and therefore

$$T_2(\alpha )=\left[\frac{K(\alpha )}{q}\right]^2.$$
(31)
The crossover time $`T_2(\alpha )`$, obtained from Eq. (31), is shown in Fig. 6 as a function of $`\alpha `$, on a log-log scale. It is found that

$$T_2=c_2\alpha ^\eta ,$$
(32)
where $`c_2=0.62\pm 0.01`$ and $`\eta =2.05\pm 0.08`$. The crossover time can be interpreted as the time required for the noise to scramble the phases and thus eliminate the structure of the constructive and destructive interference that characterize the Hadamard walk. This time can be estimated by

$$\alpha \sqrt{T_2}1$$
(33)
namely,

$$T_2\frac{1}{\alpha ^2},$$
(34)
in agreement with the simulation results. Unlike the classical random walk that produces a symmetric probability distribution for a particle that starts from the origin, the Hadamard walk exhibits an inherent asymmetry. For example, the initial state of Eq. (26) produces a distribution $`P_t(n)`$ with much more weight on the right vs. the left hand side. Moreover, the average position of the walker moves to the right at a constant speed, namely

$$\overline{n(t)}=vt.$$
(35)
where $`v=0.293\pm 0.003`$, in perfect agreement with the analytical results of Ref. Ambainis01 . In Fig. 7 we show the evolution in time of the average position $`\overline{n(t)}_\alpha `$ \[given by Eq. (22)\] for the noisy Hadamard walk with different values of the noise $`\alpha `$. The dashed linear line at the top is the result for the noiseless Hadamard walk, given by Eq. (35). Below it, from top to bottom are the simulation results for $`\alpha =0.025,0.03,0.04,0.07`$ and $`0.1`$. For short times, the average $`\overline{n(t)}_\alpha `$ of the noisy Hadamard walk tends to follow the straight line of Eq. (35), up to some crossover time, $`T_1(\alpha )`$, at which it saturates and approaches a constant value:

$$\overline{n(t)}_\alpha n_\alpha .$$
(36)
The asymptotic value, $`n_\alpha `$, decreases as the noise level $`\alpha `$ is increased. This is due to the fact that at higher noise levels, the buildup of the asymmetric pattern, which is a quantum effect, is suppressed more quickly. We define the crossover time $`T_1(\alpha )`$ as the time $`t`$ at which the line of Eq. (35), describing the noiseless Hadamard walk, intersects the asymptotic horizontal line $`\overline{n(t)}_\alpha =n_\alpha `$, namely $`vT_1=n_\alpha `$. Therefore

$$T_1=\frac{n_\alpha }{v}.$$
(37)
In Fig. 8 we present the crossover time $`T_1(\alpha )`$ vs. $`\alpha `$ on a log-log scale. It is well fitted by a power-law function of the form

$$T_1=c_1\alpha ^\rho ,$$
(38)
where $`c_1=0.20\pm 0.01`$ and $`\rho =2.06\pm 0.08`$, which is consistent with the expected asymptotic result of $`\rho =2`$, based on an argument similar to Eq. (33). We thus find that both $`T_1(\alpha )`$ and $`T_2(\alpha )`$ exhibit similar dependence on the noise level $`\alpha `$, up to a constant factor.

We will now discuss the connection between our results for the effect of unitary noise on the quantum walk and earlier results on decoherence effects Brun2002 ; Kendon2002 ; Brun2003 ; Kendon2003 ; Kendon2003b . Decoherence can be described by an additional quantum operation that applies on the chirality qubit at each step of the quantum walk. In a particular choice of the noise model, this operation is given by the completely positive map that consists of the operators Brun2003

$`\widehat{A}_0=\sqrt{p}|RR|`$
$`\widehat{A}_1=\sqrt{p}|LL|`$
$`\widehat{A}_2=\sqrt{1p}\widehat{I},`$ (39)
where $`\widehat{I}`$ is the identity operator. This set of operators satisfies

$$\underset{n=0}{\overset{2}{}}\widehat{A}_n^{}\widehat{A}_n=\widehat{I}.$$
(40)
It can be interpreted as a measurement of the chirality qubit that is performed with probability $`p`$. Analytical calculations show that in the long time limit, the first moment of the decohered quantum walk, starting in the state $`|0|R`$, saturates and approaches the value:

$$\overline{n(t)}_pn_p=\frac{(1p)^2}{p(2p)}.$$
(41)
For weak decoherence, or $`p1`$, $`n_p1/p`$, which resembles our result $`n_\alpha 1/\alpha ^2`$, if $`p`$ is replaced by $`\alpha ^2`$. This is a sensible connection since $`\alpha `$ represents some matrix elements that multiply the amplitudes while $`p`$ is a probability.

The results of Ref. Brun2003  for the second moment of the spatial distribution of the decohered quantum walk can be expressed by

$$\sigma _p(t)=K(p)\sqrt{t},$$
(42)
where

$$K(p)=\sqrt{1+\frac{2(1p)^2}{p(2p)}}.$$
(43)
Like $`K(\alpha )`$, this is a monotonically decreasing function, that in the limit of small $`p`$ takes the form $`K(p)1/\sqrt{p}`$, while $`K(p)1`$ for $`p1`$. Therefore, associating $`p`$ with $`\alpha ^2`$, one obtains the same scaling behavior of the noise effects in the two noise models for both the first and second moments of the distribution.

## V Summary

We have studied the effect of unitary noise on the discrete one dimensional quantum walk. We have shown that when the quantum walk is exposed to unitary noise it exhibits a crossover from the quantum behavior for short times to classical-like behavior for long times. For times shorter than the crossover time, the standard deviation $`\sigma _\alpha (t)`$ of the spatial distribution of the noisy Hadamard walk, scales linearly with $`t`$, like the noiseless Hadamard walk. Beyond the crossover time, $`\sigma _\alpha (t)`$ scales like $`\sqrt{t}`$, namely it acquires diffusive behavior, like the classical random walk. The crossover time was also characterized using the average position of the random walker, namely the first moment, $`\overline{n(t)}_\alpha `$, of the spatial distribution, which scales like $`t`$ for short times and saturates to a constant value for long times. In both cases the crossover time was found to scale as $`\alpha ^2`$ where $`\alpha `$ is the standard deviation of the random noise.

###### Acknowledgements.

OB thanks the Center for Quantum Computer Technology and the Department of Physics at the University of Queensland for hospitality during a sabbatical leave when this work was initiated. The work at the Hebrew University was supported by the EU Grant No. IST-1999-11234.


text:
# Numerical Simulations of Intermittent Transport in Scrape-Off Layer Plasmas

## Abstract

Two-dimensional fluid simulations of interchange turbulence for geometry and parameters relevant for the scrape-off layer of confined plasmas are presented. We observe bursty ejection of particles and heat from the bulk plasma in the form of blobs. These structures propagate far into the scrape-off layer where they are lost due to transport along open magnetic field lines. From single-point recordings it is shown that the blobs have asymmetric conditional wave forms and lead to positively skewed and flat probability distribution functions. The radial propagation velocity may reach one tenth of the sound speed. These results are in excellent agreement with recent experimental measurements.

Recently, several experimental investigations have revealed a strongly intermittent nature of particle and heat transport in the scrape-off layer (SOL) of magnetized plasmas ct:boedo ; ct:antar ; ct:terry . There are strong indications that this is caused by localized structures in the form of plasma “blobs” propagating radially far into the SOL. It has been suggested that this is due to the dipolar vorticity field caused by vertical guiding center motions in plasmas in non-uniform magnetic fields ct:krash , reflecting the compressibility of the diamagnetic current. An outstanding challenge is to give a self-consistent description of the emergence and evolution of such structures, which also capture their statistical properties. Here an attempt towards this goal is presented, yielding favorable agreement with experimental measurements. This is achieved by focusing on the collective two-dimensional dynamics perpendicular to the magnetic field while using a simplified description of particle and heat losses to limiters or end plates along open field lines.

Some of the most prominent features of experimental single-point measurements are the presence of asymmetric conditional wave forms as well as strongly skewed and flat probability distribution functions (PDF’s) of the density and temperature signals ct:boedo ; ct:antar ; ct:terry . Simple interpretations as well as advanced imaging techniques give a picture of field-aligned blobs or filaments propagating out of the bulk plasma with radial velocities up to one tenth of the sound speed. These highly non-linear thermal structures have amplitudes which significantly exceed the back-ground levels. The associated intermittent transport may have severe consequences for magnetic confinement experiments by producing large heat bursts on plasma facing components. A change of sign in the asymmetry of the fluctuation time series PDF close to the last closed flux surface (LCFS) indicates that the structures are generated by a “flapping” of the edge pressure gradient, ejecting blobs of excess particles and heat out of the bulk plasma ct:boedo . The blob structures thus seems to be generated close to the LCFS and subsequently propagate far into the SOL, where they are subject to strong damping due to parallel losses. This separation of driving and damping regions in configuration space has been discarded in several previous studies of SOL turbulence ct:sarazin ; ct:benkadda .

In this Letter we present a novel model for interchange turbulence in slab geometry and numerical solutions in qualitative agreement with experimental measurements. The model geometry comprises distinct production and loss regions, corresponding to the edge and SOL of magnetized plasmas. The separation of these two regions defines an effective LCFS, though we do not include magnetic shear in our model. In the edge region, strong pressure gradients maintain a state of turbulent convection. It is demonstrated that a self-regulation mechanism involving differential rotation leads to a repetitive expulsion of hot plasma into the SOL, resulting in asymmetric conditional wave forms, strongly non-Gaussian probability distributions and significant cross-field transport by localized objects.

Assuming cold ions and neglecting electron inertia effects, a three-field model may be derived for quasi-neutral electrostatic perturbations of the full particle density $`n`$, electric potential $`\varphi `$ and electron temperature $`T`$. Using the Bohm normalization and slab coordinates with $`\widehat{𝐳}`$ along the magnetic field we obtain ct:model ,

$$\left(\frac{}{t}+\widehat{𝐳}\times \varphi \right)\mathrm{\Omega }𝒞(p)=\nu _\mathrm{\Omega }^2\mathrm{\Omega }\sigma _\mathrm{\Omega }\mathrm{\Omega },$$

$$\frac{dn}{dt}+n𝒞(\varphi )𝒞(nT)=\nu _n^2n\sigma _n(n1)+S_n,$$

$$\frac{dT}{dt}+\frac{2T}{3}𝒞(\varphi )\frac{7T}{3}𝒞(T)\frac{2T^2}{3n}𝒞(n)=\nu _T^2T\sigma _T(T1)+S_T,$$

where time is normalized by the ion gyration period, $`1/\omega _{ci}`$, and spatial scales by the hybrid gyration radius, $`\rho _s=c_s/\omega _{ci}`$. The particle density $`n`$ and temperature $`T`$ are normalized to fixed characteristic values at the outer wall. We further define the two-dimensional advective derivative, the magnetic field curvature operator and the toroidal magnetic field, respectively, by

$$\frac{d}{dt}=\frac{}{t}+\frac{1}{B}\widehat{𝐳}\times \varphi ,𝒞=\zeta \frac{}{y},B=\frac{1}{1+ϵ+\zeta x}.$$

The vorticity is given by $`\mathrm{\Omega }=_{}^2\varphi `$, the inverse aspect ratio $`ϵ=a/R_0`$ and $`\zeta =\rho _s/R_0`$ where $`a`$ and $`R_0`$ are the minor and major radius of the device, respectively. The terms on the right hand side of the model equations describe external sources $`S`$, parallel losses along open field lines through the damping rates $`\sigma `$, and collisional diffusion with coefficients $`\nu `$. The geometry and boundary conditions are sketched in Fig. 1.

In the absence of external forcing and dissipative processes the model equations non-linearly conserves the global energy

$$E=𝑑𝐱\left[\frac{1}{2}\left(_{}\varphi \right)^2+\frac{3}{2}nT\right],$$

where the integral extends over the whole plasma layer under consideration. Thus, the curvature terms due to magnetic field inhomogeneity correctly yield a conservative energy transfer from the confined heat to the convective motions. We further define the kinetic energy of the fluctuating and mean components of the flows,

$$K=𝑑𝐱\frac{1}{2}\left(_{}\stackrel{~}{\varphi }\right),U=𝑑𝐱\frac{1}{2}v_0^2,$$
(1)
where the zero index denotes an average over the periodic direction $`y`$ and the spatial fluctuation about this mean is indicated by a tilde. The linearly damped mean flows, $`v_0=\varphi _0/x`$, does not yield any radial convective transport and hence form a benign path for fluctuation energy. The energy transfer rates from thermal energy to the fluctuating motions, and from the fluctuating to the mean flows, are given respectively by

$$F_p=𝑑𝐱nT𝒞(\varphi ),F_v=𝑑𝐱\stackrel{~}{v}_x\stackrel{~}{v}_y\frac{v_0}{x}.$$
(2)
Note that $`F_p`$ is essentially a measure of the domain integrated convective thermal energy transport, while $`F_v`$ shows that structures tilted such as to transport positive poloidal momentum up the gradient of a sheared flow will sustain the flow against collisional dissipation ct:garcia ; ct:gb ; ct:naulin .

In the following we present results from numerical simulations of the interchange model using parameters relevant for SOL plasmas. The dimensions of the simulation domain is $`L_x=2L_y=400`$ and the LCFS is located at $`x_{\text{LCFS}}=100`$. The parameters used for the simulation presented here are $`ϵ=0.25`$, $`\zeta =10^3`$, and the collisional diffusion $`\nu =5\times 10^3`$ is taken to be the same for all fields. The parallel loss rate of temperature is assumed to be five times stronger than that on the density and vorticity, $`\sigma _n=\sigma _\mathrm{\Omega }=\sigma _T/5=10^3/2\pi q`$, since primarily hot electrons are lost through the end sheaths ct:s . The damping rates $`\sigma _n`$ and $`\sigma _\mathrm{\Omega }`$ correspond to losses along one connection length $`2\pi R_0q`$ with the acoustic speed $`c_s`$, where $`q=3`$ is the safety factor at the edge. Finally, the radial line-integral of the sources $`S_n`$ and $`S_T`$ equals $`0.2`$, and the shape of the sources and parallel loss coefficients shown in Fig. 1 are given by $`\delta =16`$ and $`\xi =2`$. For the numerical solution we have employed an Arakawa scheme for the advective non-linearities and a third order stiffly stable scheme for the time integration ct:model . The spatial resolution is 512 and 256 grid points in the radial and poloidal directions, respectively, and the total time span of the simulation is $`4\times 10^6`$.

In Fig. 2 we show the typical evolution of the particle confinement $`P`$ and heat confinement $`H`$ in the edge and SOL regions, defined by

$$P_{\text{edge}}=_0^{x_{\text{LCFS}}}𝑑xn_0(x,t),P_{\text{SOL}}=_{x_{\text{LCFS}}}^{L_x}𝑑xn_0(x,t),$$

and similarly for the heat confinement $`H`$. From the figure we observe that plasma and heat gradually builds up in the edge at the same time as it is decaying in the SOL region. This is repetitively interrupted by rapid changes in which plasma and heat is lost from the edge to the SOL. More than 20% of the edge plasma may be lost during individual bursts. Also note from the figure that the normalized heat confinement is much less than the particle confinement due to the larger loss rate in the SOL region of the former. Further insight is revealed by Fig. 3 which shows the evolution of the kinetic energy contained by the mean and fluctuating motions, confer Eq. (1), as well as the collective energy transfer terms defined in Eq. (2). From the figure we observe that the convective energy and thermal transport appears as bursts during which particles and heat are lost from the edge into the SOL region. As discussed in Refs. ct:garcia, ; ct:gb, ; ct:naulin, , this global dynamics is caused by a self-regulation mechanism in which kinetic energy is transfered from the fluctuating to the mean components of the flows, and subsequently damped by collisional dissipation. As the thermal confinement is allowed to vary, this results in characteristic sawtooth oscillations ct:garcia . Note, however, the clear demonstration in Figs. 2 and 3 that thermal energy is ejected in a bursty manner from the edge and into the SOL region, where it is eventually lost by transport along open field lines.

The time-averaged profiles of the plasma density and the mean poloidal flows are shown in Fig. 4. Parallel losses in the SOL region result in average profiles peaked inside the LCFS, and weakly decaying throughout the SOL. Also shown in Fig. 4 are typical instantaneous profiles during a quiet period ($`t_{\text{quiet}}`$) and during a burst ($`t_{\text{burst}}`$). We observe significant deviations from the average profiles, with a more peaked density profile in the edge region during quiet phases. During bursts there is a substantial increase of the density profile in the SOL region due to the convective plasma transport. The temperature profiles have a similar structure but with lower amplitudes due to the larger loss rate. The self-sustained poloidal flow profiles are strongly sheared in the edge region, and have larger amplitudes during the the strong fluctuation period. This figure clearly indicates that fluctuations are driven in the strong pressure gradients in the edge region.

The statistics of single-point recordings at different radial positions $`P_i`$ indicated in Fig. 1 completely agree with experimental measurements. In Fig. 5 we present the probability distribution functions of the density signals taken from a long-run simulation containing more than a hundred burst events. The flat and strongly skewed distributions indicate the high probability of large positive fluctuations corresponding to blobs of excess plasma. The skewness and flatness factors take values up to 3 and 15, respectively, except for the outermost radial position ($`P_7`$) where very few structures arrive. At the other points the PDF’s have similar structure with a pronounced exponential tail towards large values, a characteristic feature of turbulent convection in the presence of sheared flows ct:garcia ; ct:gb .

The conditionally averaged temporal wave forms calculated from the same signals, using the trigger condition $`n>4n_{\text{rms}}`$ at each individual point, are presented in Fig. 6. An asymmetric wave form with a sharp rise and a relatively slow decay, as observed in experimental measurements, is clearly seen ct:boedo ; ct:antar ; ct:terry . The maximum density excursions significantly exceed the background level, and decay rapidly as the structures propagate through the SOL. The number of realizations for the conditional averaging decreases gradually from 125 at the innermost probe to 2 at the outermost one. Using a negative amplitude for the conditional averaging results in very few realizations, again showing the presence of blob-like structures. From two-dimensional animations and simple statistical correlation measurements we further find that the radial propagation velocity of these structures is typically about one tenth of the sound speed but with a large statistical variance. This is in excellent agreement with experimental measurements ct:boedo ; ct:antar ; ct:terry .

In Fig. 7 we show the spatial structure of the density, temperature, vorticity and electrostatic potential during a quiet period and during a burst. These correspond to the same times as the instantaneous profiles shown in Fig. 4. In the quiet period there are only weak spatial fluctuations with the plasma and heat well confined within the LCFS. During bursts we observe strong structures in all fields, which have propagated far into the SOL. Notice again the much weaker perturbation in the temperature field as compared to the density due to the difference in the parallel loss rates. Moreover, while blob-like structures are observed for the density and temperature fields, the vorticity displays a roughly dipolar structure as expected from theory and experimentally measured ct:krash ; ct:boedo .

In this Letter we have proposed a new model for interchange turbulence and demonstrated that its numerical solutions are in good agreement with that reported from experimental investigations of SOL turbulence and transient transport events ct:boedo ; ct:antar ; ct:terry . An important feature of our model is the spatial separation between forcing and damping regions which has not been accounted for in previous studies ct:sarazin ; ct:benkadda . While our model does not in detail describe the transition region between closed and open magnetic field lines, its solution clearly demonstrates a sound mechanism for the origin and nature of intermittent transport in the SOL of magnetized plasmas. Our work verifies the present experimental working hypothesis in terms of field-aligned blob-like structures propagating far into the scrape-off layer, which was recently questioned by Sarazin et alct:sarazin . The blob-like transport results in strongly skewed and flat PDF’s and asymmetric conditional wave forms. This is caused by a “flapping” of the edge pressure profile, ejecting blobs of excess particles and heat out of the bulk plasma. The excellent qualitative agreement between simulation data and experimental observations gives strong confidence that the two-dimensional structure in potential, vorticity, density and temperature of the blobs reveals the actual spatial shape of these blobs as occurring in experiments.

This work was supported by the Danish Center for Scientific Computing through grants no. CPU-1101-08 and CPU-1002-17. O. E. Garcia has been supported by financial subvention from the Research Council of Norway.

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7


text:
# Site determination and thermally assisted tunneling in homogenous nucleation

(September 8, 2003)

## Abstract

A combined low-temperature scanning tunneling microscopy and density functional theory study on the binding and diffusion of copper monomers, dimers, and trimers adsorbed on Cu(111) is presented. Whereas atoms in trimers are found in fcc sites only, monomers as well as atoms in dimers can occupy the stable fcc as well as the metastable hcp site. In fact the dimer fcc-hcp configuration was found to be only 1.3 meV less favorable with respect to the fcc-fcc configuration. This enables a confined intra-cell dimer motion, which at temperatures below 5 K is dominated by thermally assisted tunneling.

In recent years important progress has been made in understanding epitaxial growth processes on the atomic scale, from both the experimental and the theoretical side. Early experiments employed field ion microscopy (FIM), whereas in recent years scanning tunneling microscopy (STM) proved to be a perfect technique for this task.

Interestingly, the first steps in the homoepitaxial growth on fcc (111) noble-metal surfaces still belong to the unsolved problems. Compared to the more open metal surfaces, the behavior of single adatoms is more complicated on (111) surfaces, because adatoms can reside in two different binding sites, the so-called fcc and hcp sites. Whereas the fcc site continues the ABC layer stacking of the bulk the hcp site breaks this sequence. As the energy difference between the two sites and the related diffusion barriers can be very small, the site occupation preference is expected to depend critically on the size of the metal structures, i. e., it is different for monomers, dimers and trimers as was demonstrated experimentally for Ir/Ir(111) Wang90 ; Busse03 .

The experimental determination of the fcc/hcp site preference of individual adatoms is a difficult task, and only two experimental studies exist. Wang and Ehrlich determined single adatoms of Ir/Ir(111) to be located on hcp sites WangEhrlich89 ; WangEhrlich9192 , whereas in the case of Pt/Pt(111) a clear fcc site preference was found Goelzhaeuser96 . A general rule predicting the site preference was later developed for metals having a partly filled $`d`$-band Piveteau92 ; Papadia96  that can be applied in both cases. Whereas no experimental work can be found in the literature for systems without a partly filled $`d`$-shell, several theoretical studies exist. First-principles theoretical work by Stumpf and Scheffler StumpfScheffler  and Feibelman Feibelman92  predicted an hcp site preference for individual Al atoms on Al(111). Other theoretical studies on copper, silver, and gold do not find any clear hcp or fcc site preference for a single adatom adsorbed on the fcc(111) surface Boisvert95 ; Kuerpick01 . Clearly, very subtle contributions of both the electronic structure and the local relaxation have to be taken into account GiesenIbach , and therefore no general rule is available to predict the preferred binding site. Although in homoepitaxy Ir/Ir(111) is the only known experimental example of an hcp site preference so far, several similar systems may well exist.

In this work we determined the adsorption site of Cu metal adatoms to be the fcc site, by means of a unique technique of STM diffusion studies in combination with atomic manipulation, applicable also for systems that cannot be studied by FIM. Moreover a particular localized diffusion behavior of Cu dimers has been observed experimentally, with the properties as predicted earlier in a theoretical work by Bogicevic et al. Bogicevic98 . Despite the large mass of Cu atoms, this localized diffusion is dominated below 5 K by thermally assisted tunneling of Cu adatoms. Copper as a noble metal is characterized by an almost complete $`d`$-band occupation, making the existing, qualitative theoretical analysis of site preference Piveteau92 ; Papadia96  indeterminate. Instead the experimental results are complemented by an extensive set of large-scale ab-initio density functional theory (DFT) calculations of the site preference on Cu(111).

Our experiments were performed with a low-temperature STM operated at $`421`$ K. The sample was cleaned by sputtering and annealing cycles. Bias voltages refer to the sample voltage with respect to the tip. To measure possible tip effects we determined the hopping rate for different tunneling parameters and varied the tip-adatom interaction time by taking image series with fixed parameters but for different time intervals. Using atomic manipulation the monomers and dimers were positioned far away from defect sites. Great care was taken to calibrate the temperature measurements precisely.

Our calculations are within the generalized gradient approximation (GGA) using the plane-wave ultrasoft-pseudopotential code DACAPO DacapoLB . We used slab supercell geometries with a vacuum separation $`13`$ Å, and with a five- (four-)layer slab thickness and a 3-by-3 (4-by-5) atom in-surface extension to map the monomer (dimer) adsorption and dynamics. The two lowermost layers were kept frozen in a bulk-lattice structure. The adsorbates and surface atoms were relaxed to 0.05 eV/Å using an energy cutoff of 30 Ry, with dipole correction enabled. Our DFT calculations should determine the local relaxation, binding, and vibrations extremely accurately. However, ab-initio DFT can encounter difficulties in resolving the energy variation between different adsorption sites if they are in the range of only a few meV PtPuzzle01 . Consequently, we have tested the robustness of our results and used a carefully controlled sampling of $`k`$-points: a standard choice of $`2\times 2`$ with a 0.1 eV Fermi smearing tested against $`4\times 4`$ $`k`$-point calculations with a 0.025 eV smearing. In addition, we have tested the monomer site-preference result against calculations using a fully relaxed four-layer 5-by-4 unit cell.

The diffusion paths of individual adatoms were mapped by taking a series of consecutive STM images on Cu/Cu(111). Five series of more than 100 images each were taken at different temperatures between 11 and 12 K, allowing the determination of the diffusion barrier as $`E_B=(37\pm 5)`$ meV Repp00 , in good agreement with other experiments Wulfhekel96 ; Knorr02  as well as theory Stoltze94 ; Kuerpick01 . The mapping of the relative movement of one adatom in between subsequent images clearly reveals the hexagonal lattice of the substrate surface, showing that only one adsorption site is observed. As pointed out by Wang and Ehrlich WangEhrlich9192 , the conventional technique of mapping diffusion paths is not applicable for a site determination if only one of the sites is occupied. On the other hand, lateral manipulation on the atomic scale EiglerSchweizer90 ; Meyer00  can in general be used as a tool to access also metastable sites by artificial occupation of non-equilibrium configurations Heinrich02 .

For this type of extension of the mapping technique, individual adatoms were repeatedly repositioned using lateral manipulation at a temperature of 7 K. At this temperature the monomers can be positioned in both the fcc and the hcp sites, making the site determination possible. If the temperature is increased to $`10`$ K the individual adatoms regain their equilibrium position, which can now be identified as the fcc site (see Fig. 1).

In order to have a point of origin for this extended site-mapping technique, the positions were determined with respect to a monatomic copper row built up atom by atom. The row ran parallel to one of the close-packed \[$`1\overline{1}0`$\] directions. The lateral distances between the individual adatoms and the row were $`|\mathrm{\Delta }x|=n\sqrt{3/8}a_0`$, with integer values of $`n`$ for the fcc positions and fractional values of $`n`$ for the hcp positions. In addition, this shows that the atoms in the row reside in fcc positions in agreement with our theoretical calculations TrimerRowDFT .

The diffusion experiments, in which no hcp monomer occupation could be identified, revealed that the hcp-to-fcc hopping rate is at least 75 times higher than the fcc-to-hcp hopping rate, providing a lower bound for the energy difference of the two sites of $`|\mathrm{\Delta }E|kT\mathrm{ln}(75)=4`$ meV. Alternatively, the hopping from hcp to fcc sites, which sets in at a temperature that is $`\mathrm{\Delta }T1.5`$ K lower than the onset for fcc-to-hcp hopping, allows a rough estimate of the energy difference for the two sites as $`\mathrm{\Delta }EE_B\mathrm{\Delta }T/T`$. Taking the experimental uncertainties into account, this provides an upper bound of $`\mathrm{\Delta }E8`$ meV, subject to the assumption that the fcc and hcp diffusion prefactors differ by no more than a factor of ten (in fact our theory calculations yield identical prefactors of $`\nu _0=1\times 10^{12}`$s-1 monovibration ).

This result is consistent with a previous DFT study Bogicevic00  and congruent with our DFT calculations that yield an fcc preference of about 6 meV. An identical energy difference of 6 meV arises in a corresponding DFT study of monomer adsorption on a prerelaxed frozen surface, or when we increase the $`k`$-point sampling while lowering the smearing. A similar fcc preference arises when, instead, we use a four-layer 5-by-4 unit cell. These tests Bogicevic00  suggest that the monomer calculations are converged with respect to the choice of DFT parameters. The small monomer-adsorption preference reflects subtle differences in the interplay between the adsorbate-induced density reponse and the exact location of subsurface atoms, as we shall discuss further in a forthcoming publication.

In experiments, dimers do not exhibit any translational diffusion up to a temperature of 21 K. Nonetheless they already diffuse locally at a temperature as low as 7 K, as can be seen from the unstable appearance of the dimer in the STM image shown in Fig. 2(a). As predicted by theory Bogicevic98  and verified by the site determination demonstrated in Fig. 1, the local diffusion takes place in a cell centered around one atom in the topmost layer of the substrate surface. This cell of local diffusion consists of three hcp and three fcc sites around the on-top position of the substrate atom (see scheme in Fig. 1). In the following we will use the notation $`ff`$, $`hh`$, $`fh`$, $`fb`$, $`hb`$, and $`bb`$ for fcc-fcc, hcp-hcp, fcc-hcp, fcc-bridge, hcp-bridge, and bridge-bridge dimer occupation, respectively.

As can be seen in Fig. 2(b) a dimer can be stabilized by long-range interactions Repp00 ; HyldgaardPersson00 ; Fichthorn00  with another monomer positioned close to the dimer. Moreover, reducing the temperature to below 7 K, the dimer intra-cell diffusion becomes so slow that it is possible to trace the individual diffusion steps by recording a series of STM images. The stable $`ff`$ and the metastable $`fh`$ configuration can be observed (see Fig. 2(c) ). The relative probability of observation at 5 K ($`P_{fh}/P_{ff}1/20`$) yields an energy difference of only $`\mathrm{\Delta }E=kT\mathrm{ln}(P_{fh}/P_{ff})=(1.3\pm 0.5)`$ meV. The $`hh`$ configuration could not be found.

An Arrhenius plot of the dimer $`ff`$ to $`fh`$ hopping rate (Fig. 3) documents that thermally assisted tunneling Goran ; LauhonHo ; Heinrich02  enables the hopping of the dimers at very low temperatures ($`T5`$ K), despite the large mass $`m_{\mathrm{Cu}}`$ of a Cu atom. In a simple model with a parabolic barrier at the transition state $`fb`$, characterized by an imaginary frequency $`\omega _{fb}`$ through $`E=E_B+1/2m_{\mathrm{Cu}}\omega _{fb}^2x^2`$, the transmission probability due to tunneling for an energy $`ϵ`$ is known to be approximately $`\mathrm{exp}(2i\pi (E_Bϵ)/\mathrm{}\omega _{fb})`$ HillWheeler . In the case of thermally assisted tunneling this probability has to be multiplied by the probability for thermal activation to the energy $`ϵ`$ of $`\mathrm{exp}(ϵ/kT)`$. This results in a critical temperature $`T_C=\mathrm{}\omega _{fb}/2i\pi k`$, which marks the limit between thermally activated hopping and thermally assisted tunneling. Above this temperature we obtain the conventional Arrhenius behavior with the activation energy $`E_B`$, whereas for the low-temperature regime, only thermally assisted tunneling to the ground state is relevant, and the slope in the Arrhenius plot is determined by the energy difference $`\mathrm{\Delta }E`$ of the two ground states. Only in the temperature region just below $`T_C`$ do all intermediate processes come into play (see fits in Fig. 3). The fit to the data yields $`\nu _0=810^{11\pm 0.5}`$ s-1 and $`E_B=18\pm 3`$ meV. Assuming the barrier to be cosine-shaped, $`\omega _{fb}`$ and thus $`T_C`$ can be estimated from the distance $`d=1.47`$ Å between the fcc and the hcp sites and the barrier height $`E_B`$ as $`T_C\sqrt{E_B/2m_{\mathrm{Cu}}}\mathrm{}/kd=6.0`$ K, in agreement with experiment. The surprising observation of the tunneling of a heavy Cu atom can be understood when considering that the atom not only has a very low barrier but also an exceptionally small distance to overcome. In contrast to the recently reported tunneling of carbon adatoms Heinrich02 , we do not see any influence of descrete vibrational levels of the dimer. This indicates a very strong vibronic coupling of the dimer to the substrate, clearly ruling out the possibility of a Cu2/Cu(111) quantum rotor Bogicevic98 .

From calculations using $`2\times 2`$ $`k`$-point sampling (Fig. 4, dashed lines), we predict an intra-cell dimer-energy variation that is in qualitative agreement with the original Al-dimer study Bogicevic98 . As in these previous theoretical calculations and consistent with the experiments, we find that intra-cell rotation remains possible even when low temperatures inhibit the net inter-cell diffusion, for which we determine a concerted-sliding barrier $`E_{\mathrm{B},\mathrm{dif}}\text{190\hspace{0.17em}meV}`$. Using the more accurate $`4\times 4`$ $`k`$-point sampling we predict an energy variation (Fig. 4, solid lines), in even better agreement with experimental observations. The $`4\times 4`$ $`k`$-point calculations establish a near degeneracy between the $`ff`$, $`fh`$, and $`hh`$ dimer configurations and a significantly reduced effective rotation barrier, i. e. the energy of the $`fb`$ configuration. The accuracy improvement obtained by using the $`4\times 4`$ $`k`$-point sampling in the dimer-configuration survey is selective and correlated with the dimer-induced surface distortion as also previously found in the Al-addimer study Bogicevic98 . We attribute the experimental observation of dimer confinement to the rotation process, $`fffbfh`$ (and back), which has a lower barrier than the $`fffbfhhbhh`$ rotation and the direct $`ffbbhh`$ transition.

Trimers were found to be always of fcc-fcc-fcc type and stable in the temperature range of $`521`$ K.

In summary, we have shown that individual copper adatoms as well as adatoms in dimers, trimers, and larger structures prefer fcc over hcp sites on Cu(111). The experimentally observed difference in binding energy for the monomer is only about 6 meV. Moreover, in the case of copper dimers, we were able to observe a confined intra-cell metal dimer motion for a wide temperature range of $`421`$ K, which at low temperatures is enabled by thermally assisted tunneling. Even though the small energy differences inherent in this system may attain the limits of current DFT accuracy, a detailed DFT analysis yields monomer and dimer energetics that are in very good agreement with the experiments (see Table 1).

We thank Prof. G. Ehrlich and Prof. M. Scheffler for their very valuable comments, and acknowledge partial funding by the EU TMR project “AMMIST”, the EU IST-FET project “NICE”, the Deutsche Forschungsgemeinschaft Project No. RI 472/3-2, the Chalmers UNICC (for supercomputing), and the Swedish Foundation for Strategic Research via ATOMICS (PH).