text:
# The probability distribution of the conductance in anisotropic systems

## Abstract

We investigate the probability distribution $`p(g)`$ of the conductance $`g`$ in anisotropic two–dimensional systems. The scaling procedure applicable to mapping the conductance distributions of localized anisotropic systems to the corresponding isotropic one can be extended to systems at the critical point of the metal–to–insulator transition. Instead of the squares used for isotropic systems, one should use rectangles for the anisotropic ones. At the critical point, the ratio of the side lengths must be equal to the squre root of the ratio of the critical values of the quasi–one–dimensional scaling functions. For localized systems, the ratio of the side lengths must be equal to the ratio of the localization lengths.

The presence of disorder may allow a system to make a transition from metallic to insulating behaviour by varying the Fermi energy in an energy range where both extended and localized states are found, separated by a mobility edge. Characterizing this transition, one can employ transport properties, such as the conductance, or properties of the system’s eigenstates, such as the correlation length for extended, metallic states or the localization length $`\xi `$ for insulating states. At the mobility edge, a determination of the complete probability distribution $`p(g)`$ of the conductance $`g`$ (in units of $`e^2/h`$) is needed. The critical point of the transition from metallic states to Anderson localized ones is of particular interest. The distributions are well–known to be normal and log–normal off the mobility edge towards the extended and the localized regime respectively, whereas the exact form of the critical distribution is still under investigation. For example, contrary to expectations, the critical distribution seems to vary even within the same universality class, depending on the boundary conditions perpendicular to the direction in which transport occurs. Also, questions about the exact form of the large–$`g`$ tail ($`g>1`$) remain unanswered. Where calculations in $`2+\epsilon `$ dimensions indicate higher cumulants to diverge with system size, leading to a power law tail, numerical calculations in three dimensions and analytical results for quasi–one–dimensional wires show an exponential decay.

Anisotropic systems have recently been the focus of particular attention. It is generally accepted that anisotropy does not change the universality class and that isotropic results can be recovered by performing a proper scaling of the anisotropic results. For anisotropic systems in a localized state, it is reasonable to assume that scaling the dimensions of the system by the corresponding localization lengths will make the system effectively isotropic. This procedure has been applied successfully to the scaling function $`\mathrm{\Lambda }=\lambda _M/M`$, which is a function of $`\xi /M`$, where $`\lambda _M`$ denotes the finite size localization length of a quasi–one–dimensional strip of finite width $`M`$ and $`\lambda _M\xi `$ as $`M\mathrm{}`$. It was also shown that the same scaling procedure works for the probability distribution in such a system.

In order to test the approach for critical states one must either face the numerical challenge of large three–dimensional systems or take into account additional interactions (beyond the disorder potential) such as spin–orbit coupling. Another possibility is the introduction of external magnetic fields as e.g. in integer quantum Hall systems or tight binding models with random magnetic flux. However, as a result of Anderson localization, extended states do not exist in two–dimensional systems of non–interacting electrons in a magnetic field, except at a singular energy near the center of each of the Landau subbands. At these critical energies $`E_c`$ the localization length $`\xi `$ diverges with a critical exponent $`\nu `$: $`\xi |EE_c|^\nu `$.

Because significant finite size effects have to be expected, we decided to concentrate our research on two–dimensional systems although the exact form of the critical distribution of the conductance depends on the dimensionality of the system. The investigation of the self–averaging quantity $`\xi `$ in integer quantum Hall systems yielded very encouraging results, supporting the expectation that quantities of anisotropic systems can indeed be mapped to isotropic values by a simple rescaling scheme.

In this paper, we show a method of mapping the probability distributions of the conductance of anisotropic two–dimensional systems with a magnetic field perpendicular to the plane or with spin–orbit coupling to the probability distribution of the conductance for the corresponding isotropic system at the critical point, using a tight binding model. It turns out that the ratio of the squares of the side lengths $`L_x`$, $`L_y`$ of the anisotropic system should be chosen equal to the ratio of the critical values $`\mathrm{\Lambda }_x^c`$, $`\mathrm{\Lambda }_y^c`$ of the quasi–one–dimensional scaling functions:

$$\frac{L_x^2}{L_y^2}=\frac{\mathrm{\Lambda }_x^c}{\mathrm{\Lambda }_y^c}$$
(1)
In the following, we first describe the models and the numerical method we employed. Then we present and discuss our numerical results and finally summarize the conclusions of this work.

The tight binding model uses the Hamiltonian

$$=\underset{n,\tau }{}|n\tau \epsilon _nn\tau \left|+\underset{n,\tau ,n^{},\tau ^{}}{}\right|n\tau V_{n,n^{}}n^{}\tau ^{}|$$
(2)
where $`n,n^{}`$ denotes the lattice site. Without spin–orbit interaction the “variables” $`\tau ,\tau ^{}`$ take on only one value and the hopping integrals $`V_{n,n^{}}`$ are scalar, otherwise they are $`2\times 2`$ matrices and the spin variables take on the values $`1`$ or $`1`$. In either case the site energies $`\epsilon _n`$ are independent of $`\tau `$ and we take into account interactions only between neighbouring lattice sites.

An external magnetic field enters the Hamiltonian via its vector potential $`𝐀`$ ($`\times 𝐀=𝐁`$), which appears in the phases of the hopping integrals:

$$V_{n,n^{}}=t_{n,n^{}}^0e^{2\pi i(e/h)_{𝐫_n}^{𝐫_n^{}}𝐀(𝐫)d𝐫}$$
(3)
The integral connects the lattice sites $`n`$ (at $`𝐫_n`$) and $`n^{}`$ (at $`𝐫_n^{}`$) in a straight line. For the sytems under consideration, where the magnetic induction $`𝐁`$ is perpendicular to the plane of the two–dimensional lattice, the gauge for the vector potential can be chosen such that the phases vanish in the direction perpendicular to $`𝐀`$ and are integer multiples of some number $`2\pi \alpha `$ in the direction parallel $`𝐀`$. The value of the parameter $`\alpha `$ then completely characterizes the influences of the magnetic field on the system. For rational $`\alpha `$, the denominator determines the number of bands in the density of states of the system without disorder.

The Evangelou–Ziman model incorporates spin–orbit coupling by using the following hopping integrals:

$$V_{n,n^{}}^{\tau ,\tau ^{}}=t_{n,n^{}}^0\left[\delta _{\tau ,\tau ^{}}+\mu i\underset{\nu }{}\sigma _{\tau ,\tau ^{}}^\nu t_{n,n^{}}^\nu \right]$$
(4)
where $`\nu =x,y,z`$ and $`\sigma ^\nu `$ are the Pauli matrices. The parameter $`\mu `$ characterizes the strength of the spin–orbit interaction.

Both systems may be made anisotropic by chosing the value of $`t^0`$ to be different in the two directions within the plane. Otherwise this parameter is a constant, independent of lattice site $`n`$. We bring disorder to the system by chosing all the site energies independently from a rectangular distribution of width $`W`$ centered at $`0`$, so that $`W`$ is a measure of the strength of the disorder. The parameters $`t^\nu `$ are also randomly selected from a uniform distribution on $`[1/2,1/2]`$. The energy scale is set by the larger of the two values for $`t^0`$, which is taken to be unity.

We calculate the conductance from the Landauer formula

$$g=\text{Tr}\left(\mathrm{t}^{}\mathrm{t}\right)$$
(5)
where $`\mathrm{t}`$ is the transmission matrix. We suppose two semi–infinte leads are attached to opposite sides of the sample. Then $`\mathrm{t}`$ determines the transmission of an electron through the sample. The numerical procedure is based on the algorithms published by Ando and Pendry.

The critical conductance distributions we calculated for isotropic quantum Hall systems at different disorder strengths show that finite size effects become stronger the weaker the disorder. Where systems with $`W=4.0`$ and $`W=2.0`$ show a basically size–independent critical distribution of the conductance for squares of $`64\times 64`$ lattice sites, at $`W=0.5`$ finite size effects are still somewhat noticeable up to systems with $`192\times 192`$ lattice sites. The anisotropic quantum Hall systems we investigated are characterized by $`\alpha =\frac{1}{8}`$. The anisotropies we chose were $`t_x^0/t_y^0=0.5`$ at $`W=0.5`$ and $`t_x^0/t_y^0=0.8`$ at $`W=0.1`$. The latter was chosen mainly because we already had the data for the quasi–one–dimensional scaling function. As the disorder is even weaker than in the first case, finite size effects are even stronger, and even at $`240\times 240`$ lattice sites the conductance distribution is far from the one we expect from our calculations of isotropic systems. Therefore, we will not be able to show that our procedure maps the two anisotropic conductance distributions to the critical distribution of isotropic systems for this extreme case. We will however be able to prove the somewhat weaker claim that our method transforms the two anisotropic distributions so that both have the same shape. In a square system, one expects that the distribution in the difficult hopping direction shows a more localized character than the one in the easy hopping direction. In an isotropic system, the distribution obviously cannot depend on the direction of transport. By making the system rectangular rather than square, i.e. shorter in the difficult hopping direction, it should be possible to obtain distributions in the two directions that are the same, thus making the anisotropic system effectively behave isotropically.

The task now is “How do we choose the correct ratio of side lengths of the rectangle?” From the research on localized systems we know that in those cases, the ratio should be equal to the ratio of the localization lengths:

$$\left(\frac{L_x}{\xi _x}=\frac{L_y}{\xi _y}\right)\frac{L_x}{L_y}=\frac{\xi _x}{\xi _y}$$
(6)
for localized systems, as these are obviously the appropriate length scales in their respective directions. This is of no use for critical systems as both localization lengths diverge at the transition. A closely corresponding non–diverging quantity is, however, available in the scaling function $`\mathrm{\Lambda }_M=\lambda _M/M`$, which has a finite critical value, independent of the system width $`M`$. The finite size localization lengths $`\lambda _{M,x}`$ and $`\lambda _{M,y}`$ have $`\xi _x`$ and $`\xi _y`$ respectively as their large–$`M`$ limits, and for big enough systems we can approximate Eq. 6 by

$$\frac{L_x}{L_y}=\frac{\lambda _{M,x}}{\lambda _{M,y}}$$
(7)
for localized systems at “large enough” $`M`$. The meaning of $`M`$ in this context would be that of the system width perpendicular to the direction in which the localization length is measured, i.e. $`M=L_y`$ in $`\lambda _{M,x}`$ and $`M=L_x`$ in $`\lambda _{M,y}`$. Now, by multiplying both sides of Eq. 7 by $`L_x/L_y`$ we have

$$\frac{L_x^2}{L_y^2}=\frac{L_x}{\lambda _{M,y}}\frac{\lambda _{M,x}}{L_y}=\frac{\mathrm{\Lambda }_{M,x}}{\mathrm{\Lambda }_{M,y}}$$
(8)
for “large” localized systems. Now $`\mathrm{\Lambda }_M`$ is a continuous function of $`E`$ and for “large enough” systems at $`E_c`$ should have reached its critical value. Therefore, we arrive at the conclusion that

$$\frac{L_x}{L_y}=\sqrt{\frac{\mathrm{\Lambda }_x^c}{\mathrm{\Lambda }_y^c}}$$
(9)
should be the correct ratio for critical systems in order to make them behave like isotropic ones. Noting the relationship

$$\sqrt{\mathrm{\Lambda }_x^c\mathrm{\Lambda }_y^c}=\mathrm{\Lambda }_{\text{iso}}^c$$
(10)
we can write Eq. 9 in the alternate form

$$\frac{L_x}{L_y}=\frac{\mathrm{\Lambda }_x^c}{\mathrm{\Lambda }_{\text{iso}}^c}=\frac{\mathrm{\Lambda }_{\text{iso}}^c}{\mathrm{\Lambda }_y^c}.$$
(11)
We tested this prediction on the system with $`t_x^0/t_y^0=0.5`$ and $`W=0.5`$, where the ratio $`\sqrt{\mathrm{\Lambda }_x^c/\mathrm{\Lambda }_y^c}`$ is roughly $`1.5`$. The result is shown in Fig. 1 together with the critical distribution for the isotropic system. The agreement is very good. For the other system with $`t_x^0/t_y^0=0.8`$ and $`W=0.8`$, we have to deal with stronger finite size effects and cannot expect to approach the form of the critical distribution we see in Figs. 1 for reasonable system sizes. Instead we merely show in Fig. 2 how the critical distributions change with the ratio of side lengths. The best value for the ratio according to Eq. 9 would be roughly $`1.23`$. Fig. 2 shows results for ratios of $`1.0`$, $`1.25`$ and $`1.5`$. The averages of $`\mathrm{ln}(g)`$ for the easy hopping direction decrease with increasing ratio from $`4.09`$ for the square to $`4.40`$ and $`5.30`$, while the averages for the difficult hopping direction increase from $`4.94`$ for the square to $`4.18`$ and $`4.03`$. Similarly, the standard deviations increase for the easy hopping direction from $`1.90`$ for the square to $`2.07`$ and $`2.30`$, while they decrease for the difficult hopping direction from $`2.24`$ for the square to $`2.01`$ and $`1.91`$. The values for a ratio of $`1.25`$ are not equal but reasonably close, so that for larger systems, where a ratio of $`1.23`$ might be practicable, we expect a better agreement of the two probability distributions.

Taking the best–ratio rectangle as the “undeformed” base, we can also see from Fig. 2 that similar “deformations” have similar effects in the two directions, that is, reducing the ratio by a factor $`\gamma `$ in one direction will cause the ensemble average $`<\mathrm{ln}(g)>`$ in that direction to increase, and the standard deviation to decrease, while the trend is opposite in the perpendicular direction. However, by reducing the ratio by the same factor $`\gamma `$ in the other direction, will result roughly in the same distributions as before, but with the one associated with the easy direction before now assigned to the difficult direction and vice versa.

That the same procedure also works for sytems with spin–orbit coupling is shown in Fig. 3, where we plot the conductance distribution for an isotropic system with $`\mu =1.0`$ at $`E=0.1`$ and $`W_c=6.7`$ together with that of two anisotropic systems, one with $`\mu =1.0`$ and $`t_x^0/t_y^0=0.1`$ at $`E=0.1`$ and $`W_c=1.6`$, the other with $`\mu =1.0`$ and $`t_x^0/t_y^0=0.2`$ at $`E=0.1`$ and $`W_c=2.6`$. The ratio of sidelengths, according to Eq. 9, should be $`23.0`$ for the latter. We chose $`40\times 40`$ lattice sites for the isotropic system, $`10\times 230`$ lattice sites for the strongly anisotropic one and $`20\times 185`$ lattice sites for the anisotropic system with the weaker anisotropy. Again the agreement is very good.

We have shown that the scaling procedure applicable to mapping the conductance distributions of localized anisotropic systems to the corresponding isotropic one can be extended in a straightforward manner to systems at the critical point of the Anderson localization–delocalization transition in both unitary and symplectic two–dimensional systems. Instead of the squares used for isotropic systems, one should use rectangles for the anisotropic ones, with a ratio of side lengths equal to the square root of the ratio of the critical values of the quasi–one–dimensional scaling function.

Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract No. W–7405–Eng–82. This work was supported by the Director for Energy Research, Office of Basic Science. P.M. would like to thank IITAP and Ames Laboratory at Iowa State University for their hospitality and support and the Slovak Grant Agency for financial support.


format: text

text:
# Ferromagnetic instabilities in atomically-thin lithium and sodium wires

\[

## Abstract

Using density functional theory the ground state structural, electronic, and magnetic properties of monatomic lithium and sodium chains with low average density are investigated. A metallic, zigzag ground state structure is predicted but, most interestingly, stable equilibria for chains under tension are predicted to be ferromagnetic, which can be traced to exchange effects arising from occupation of the second subband as a function of the interatomic distance.

\]

Substantial variations in structural, electronic, and magnetic properties of solids have been predicted to occur with decreasing dimensionality, and with the advent of new techniques in atomic manipulation and nanofabrication, it is now possible to realize quasi-one dimensional systems experimentally. Using photoemission spectroscopy, both Luttinger-liquidlike properties  and Peierls instabilities  have already been observed in linear monatomic chains formed on insulating substrates. Despite confirmation of such unusual structural and electronic properties, the existence of magnetic order in quasi-one-dimensional systems remains somewhat controversial. Generalizing an earlier prediction of Paul, the Lieb-Mattis theorem formally prohibits ferromagnetic ground states for strictly one dimensional systems in the absence of spin- or momentum-dependent forces. Although strictly valid for single-band systems, this theorem is no longer applicable when two or more transverse subbands are occupied, and in fact we report here the results of first-principles calculations indicating a ferromagnetic ground state for lithium chains, as well as stable magnetic equilibria for quasi-one dimensional sodium chains under tension. The unexpected development of magnetic order in these classic nearly-free electron systems can be traced to exchange effects arising from the changing occupation of a second subband as a function of the interatomic distance. Accordingly, at high  and now at low densities, lithium and sodium depart radically from the simple-metallic behavior that typifies their normal density electronic structure.

To examine lithium and sodium chains over a wide range of interatomic distances, we use a plane wave implementation of density functional theory  within the local spin-density approximation, the Vienna ab initio Simulation Package (VASP)  Our Li and Na ultrasoft pseudopotentials treat only the single outer $`s`$ electron as valence; they also include non-linear core corrections, and both reproduce previously calculated structural and electronic properties of bulk Li and Na. We use a supercell containing four atoms aligned along only one direction (defining the chain) and include 20 Å of vacuum in the orthogonal directions to minimize interactions between neighboring chains. For several different values of the chain-length per atom, $`x`$, we optimize both the interatomic distance, $`d`$, and the bond-angle, $`\alpha `$; the resulting binding energies appear in Fig. 1. As can be seen, the energy has a pronounced minimum in both cases (around $`x=1.7`$ Å for sodium and $`x=1.55`$ Å for lithium), at which point $`d=3.47`$ Å and $`\alpha =61^\mathrm{o}`$ for sodium and $`d=2.85`$ Å and $`\alpha =63^\mathrm{o}`$ for lithium. Therefore we conclude that linear light-alkali chains are unstable to a lateral buckling into a planar zigzag or polymer-like configuration. Out-of-plane deformations do not result in lower energies. Previous experimental work by Whitman et al. observed that monatomic cesium chains on GaAs and InAs substrates form a zig-zag structure, which then remained stable under tension. Recent first-principles calculations have also suggested a similar zigzag structure in monatomic gold chains ; a similar zigzag geometry was proposed to explain expanded liquid Rb and Cs.

Although both alkali chains satisfy the necessary condition for insulating behavior with two monovalent atoms per unit cell, band overlap nevertheless prevents the systems from being insulating. When perfectly linear, the electronic structure of the each chain consists of a single half-filled $`s`$-like band, with two degenerate, unoccupied $`p`$-like transverse subbands lying above. The electronic structure of the zigzag chain, however, exhibits two partially occupied subbands, one predominately $`s`$-like and one $`p`$-like. Formation of the zigzag structure lifts the two-fold degeneracy of the second subband (corresponding to the two transverse $`p`$-like bands of the linear chain), one of which then drops below the Fermi energy, stabilizing the zigzag structure over the linear chain. This symmetry-breaking also increases the coordination of each atom from two to four nearest-neighbors, lowering Madelung-type terms and further stabilizing the zigzag arrangement. Additionally, we note that distances between neighbors in both chains are longer than the bond length of the molecule (2.99 Å for sodium and 2.71 Å for lithium) but smaller than the shortest interatomic distance in their equilibrium BCC solid phases (3.66 Å and 2.92 Å respectively), reflecting an expected increase in bond length with coordination number.

Along with the energy minimum at a length of $`x=1.7`$, the energy of the sodium chain has two noticeable kinks (at $`x=1.9`$ and 2.2 Å). Similar features are found in the lithium chain, though only one kink is observed (at $`x=1.95`$ Å) upon stretching (Fig. 1(b)). In both chains, we record complementary discontinuous changes in both $`d`$ and $`\alpha `$, each associated with major atomic rearrangements (as will be illustrated in Fig. 4 below). The additional metastable equilibria observed in sodium exhibit surprising electronic and magnetic properties. As $`x`$ increases, the second subband becomes progressively less occupied, and this has significant consequences for the emerging magnetic properties, as we now show. (Although in the following our discussion will focus primarily on the sodium chains, our calculations indicate that equivalent conclusions are obtained for lithium.) The band structure corresponding to sodium chains for $`x=2`$ Å is shown in Fig. 2 and, remarkably, the bands for electrons with opposite spins are shifted with respect to one another, indicating that the sodium chains have a non-zero net spin polarization.

Interestingly the lithium chain already has a net spin polarization at its equilibrium stabilizing length, thereby accounting for the absence of the first kink in Fig. 1(b). The spin-split $`p`$-like band, corresponding to antisymmetric orbitals oriented transverse to the direction of the chain, allows for electrons of a single spin character only and is thus solely responsible for the magnetic polarization of the chain. The occupation number corresponding to this subband is so low ($`0.15`$ $`e`$/atom) that the exchange energy actually overcomes the excess of kinetic and correlation energies, resulting in the ferromagnetic state. We further note that the density of states (DOS) is observed to vary roughly as $`1/\sqrt{Eϵ_m^\sigma }`$ (where $`ϵ_m^\sigma `$ is the energy of subband $`m`$ with spin $`\sigma `$ at the $`\mathrm{\Gamma }`$ point), as is expected for free electrons in one dimension. This behavior arises from van Hove singularities at the bottom of each subband, which are especially prominent in one dimension. We note in passing that although our calculations carried out within the generalized gradient approximation (GGA)  result in a slightly larger interatomic distance $`d`$, they do not change the overall notably “polymeric” structure of the chains. The closing of the second subband and corresponding transition to the nonmagnetic state is found to occur at the same chain length, but the magnetic instability remains for a larger range of chain lengths within the GGA.

The valence charge and spin polarization for sodium chain with $`x=2.0`$ Å appear together in Fig. 3. Notice that the charge is mainly confined to the region between the ions, and is almost uniform along the longitudinal direction of the chain. There does not appear to be any directional bonding between the ions: covalent bonding will therefore not account for the stability of the zigzag structure. The minima on either side of the charge density maxima correspond to the location of the ions, where the valence charge density is small because of the pseudopotential approximation. Because of the $`p`$-symmetry of the second subband, the spin density also exhibits interesting character: in Fig 3(b) two separate maxima can be seen on both sides with a minimum in the center. The spin polarization in this valley has opposite sign as compared to the net magnetic moment. The magnetic polarization of the chain is therefore highly anisotropic and, as previously indicated, the chain is ferromagnetic along $`x`$ but has an antiferromagnetic character along the transverse direction of the chain. Similar charge and spin densities can be observed for lithium zigzag chains.

When the chain length is extended to $`x=2.2`$ Å , the second subband rises above the Fermi energy and the gap between spin-split subbands closes; the ferromagnetic order vanishes, resulting in the final kink observed in Fig. 1(a) for the sodium chain. We may also note that for the paramagnetic state (having only one occupied subband) the energy is well fit to a polynomial function of $`x`$ (as in Fig. 1), and the minimum energy obtained from this fit, corresponding to the cohesive energy of the chain with only the first subband occupied, is higher than the value calculated for the ground state of the chain. Thus filling the second band decreases the total energy and stabilizes the zigzag chain. The net magnetic moment per atom and the optimal geometry of the zigzag chains as a function of $`x`$ for both alkali chains is shown in Fig. 4. The maximum value of the net magnetic moment of the chain (per conduction electron) is $`0.18\mu _B`$ ($`0.16\mu _B`$ for lithium), and it decreases with increasing $`x`$; and for $`x>2.2`$ Å, only one subband is occupied and the chain is once again paramagnetic. The discontinuities in $`d`$, corresponding to the kinks observed for the energy and major atomic rearrangements, can be seen to reflect the onset of magnetic transitions.

According to our results, the spontaneous magnetization in these zigzag chains originates from the exchange between electronic states in the second subband; these results are therefore not in violation of the Lieb-Mattis theorem, which holds strictly for systems in which only a single subband is occupied. Recent calculations have shown that sodium quantum wires, described within the jellium model, undergo spontaneous magnetization for certain radii. In particular, they have predicted that the ground state is fully spin polarized when only the lowest subband is occupied. The effects of band structure and associated orbital hybridization, absent from the jellium model, clearly have an important influence on the stability of the chain. When the electron-ion interaction and atomic structure are explicitly considered, we have found that the ground state of the chain is paramagnetic when only the second subband is occupied. At larger chain lengths (i.e., $`x>2.2`$ Å ) the first subband is completely filled and although the Fermi wave vector makes contact with the edge of the Brillouin zone, we do not observe a Peierls transition. Instead, we find an intermediate state where the chain remains metallic and paramagnetic until $`x=3.5`$ Å , above which dimerization drives the metal-insulator transition.

In summary, we observe that one-dimensional lithium and sodium chains are unstable to zigzag polymer-like structures. The transition to the zigzag structures increases the coordination number and lowers the energy of a minimally-occupied $`p`$-like subband (possessing a large van Hove peak). The spin-splitting of a $`p`$-like subband moves the van Hove peak to lower energies in lithium, and a ferromagnetic ground state is predicted; for sodium, this ferromagnetic ordering is stable only when the chains are under tension. Our calculations clearly establish an intimate connection between the existance of magnetic order and zigzag arrangements. But the fact that magnetic order is found at all in these structures raises the important possibility of the occurance of more general ground states, such as, for example, non-collinear magnetic order, or even static spin waves as originally proposed for the alkalis by Overhauser . Accordingly, future extension of this work along these lines would be of considerable interest. The results presented here for lithium and sodium may not necessarily be extended to chains composed of heavier elements. For the sake of comparison, monoatomic gold chains also show the same zigzag structures and the energy shows two minima with the intermediate maximum associated to the closing of a conduction band. However, we have found a very small magnetic instability in gold chains, $`\mu _{\mathrm{max}}0.02\mu _B`$, and we believe that is related to the presence of $`d`$ bands near the Fermi energy. The possibility of magnetic chain- or stripe-like structure is also predicted to arise in an important system which might be seen to underlie the alkali chains considered here, namely the two-dimensional interacting electron gas at low densities . Future experiments on inert substrates (e.g., cleaved rare-gas crystals), chosen so that they not interfere with the $`p`$-like subbands, would clearly be of very significant interest.

This work was supported by the National Science Foundation (DMR-9988576). This work made use of the Cornell Center for Materials Research Shared Experimental Facilities, supported through the National Science Foundation Materials Research Science and Engeneering Centers program (DMR-0079992). A.B. would like to acknowledge financial support from the Spanish Ministerio de Educación y Cultura under Fulbright Grant No. FU-99-30656084.


text:
# Do Hard Spheres have Natural Boundaries?

## 1 Introduction

Recently it has been shown  that the susceptibility of the two dimensional Ising model at $`H=0`$ has a natural boundary in the complex temperature plane. Furthermore the first 323 terms in the high and low temperature series expansions have been generated ,. This discovery of a natural boundary gives new insights into phase transitions. In this note I explore possible consequences for the phenomenon of hard sphere freezing.

The freezing transition in the classical hard sphere system has been observed in both molecular dynamics and Monte Carlo computations since the 50’s , and 60’s ,  and by now it is well established  that there is a first order transition with the coexistence of a fluid phase with a dimensionless volume fraction of $`\eta _f=0.494`$ and a solid phase with a volume fraction of $`\eta _s=.545`$ where $`\eta `$ is defined as $`\eta =\rho v_{hs}=\rho \frac{4}{3}\pi r_{hs}^3`$ with $`r_{hs}`$ the radius of the hard spheres and $`\rho `$ is the number density $`N/V.`$ A non equilibrium glass transition occurs at $`\eta _g=0.58.`$ By comparison the volume fraction of the closest packed configuration (the FCC lattice) is $`\eta _{cp}=\frac{\pi }{3\sqrt{2}}=0.7405`$ which was conjectured by Kepler in 1611 and recently proven by Hales . There is also a volume fraction $`\eta _{rcp}=0.638`$ at which random close packing occurs.

Because the volume fractions $`\eta _f`$ and $`\eta _s`$ lie well below the densities $`\eta _g`$ and $`\eta _{rcp}`$ where possible non ergodic behavior might occur the freezing transition is expected to be described by the partition function of the canonical ensemble. Yet , somewhat surprisingly, after more than three decades of effort the existence of this transition is still considered controversial. This is very well stated in a recent paper  “ Theory can explain the transition only “afterwards”, that is, properties of the coexisting phases should be inserted into the theories in advance, and the coexisting densities follow by minimizing the free energy functional”.

The theory which is usually used to describe the transition (see for example ref. and references contained therein) has remained virtually unchanged for over 30 years ,. Typically the high density phase is described by a cell (or free volume) approximation  which has the virtue that it is the exact limiting form of the free energy as $`\eta \eta _{cp}.`$ The fluid phase is described by solving a liquid integral equation  based on the Percus–Yevik approximation or the scaled particle approximation  to the fluid free energy. A more modern approach uses density functional theory . These high and low density expressions are then extrapolated out of their regions of validity and in this extrapolated region the densities of coexisting phases are found by a Maxwell construction. This procedure is guaranteed to give a first order transition but cannot actually be called predictive.

On the other hand this description of a first order transition has long ago been called into question by the droplet models of condensation of Fisher  and Langer  which predict an essential singularity at the coexistence density.

These two points of view of first order transitions are not in harmonious agreement and the state of knowledge is perhaps well described by ref.$`\mathrm{}`$ from a fundamental viewpoint they \[the existing theories of freezing\] are still unsatisfactory. The splitting of a hard sphere system into a fluid and a solid branch is, evidently, hidden in the partition function of the system. A thorough theory of hard sphere freezing should therefore identify this property that leads to the observed symmetry breaking”

The purpose of this note is to use the intuition obtained from the recent Ising model computations - to reexamine the hard sphere freezing transition and to propose that the phase transition is due to the formation of a natural boundary in the free energy in the complex density plane.

## 2 Natural Boundaries

Yang and Lee  pioneered the study of phase transitions by means of looking at the zeroes of the partition function on a finite size lattice. These zeros will in general not be on the positive real axis (fugacity, temperature or density) but may approach the axis in the limit as $`V\mathrm{}.`$ We may classify several of the possible limiting behavior of zeroes as follows:

1) The zeroes can lie on a curve in the complex plane and pinch the real axis at a point (say $`T_c`$ for example). In the limit $`V\mathrm{}`$ the curve of zeroes in the complex plane becomes a branch cut and the point of pinching becomes a branch point. This happens in the complex temperature plane for all integrable models and forms the basis of the theory of second order phase transitions. In this picture the thermodynamic functions are singular only at $`T_c,`$ the point of phase transition.

2) The zeroes can lie in some area of the complex plane, pinch the real axis in some line segment and cut the plane into two disconnected pieces each of which may be analytically continued into the region occupied by the zeroes in the finite size system. In this picture the thermodynamic functions such as the free energy and the pressure are analytic at the coexistence density $`\eta _f`$ and the virial expansion will converge for densities greater than $`\eta _f.`$ This picture forms the basis of the theory of first order phase transitions and is exactly the opposite of what happens in a second order transition.

3) The zeroes may pinch the real axis at either a point or a line segment and form a natural boundary in the complex plane. This is the phenomena discovered by Orrick, Nickel, Guttmann and Perk - in the susceptibility of the Ising model.

The first mechanism of phase transition is explicitly seen in every integrable model whose partition function has been computed. However, the computations of ref. - strongly lead to the conjecture that for non integrable models with a second order phase transition there will always be a natural boundary accompanying the singularity at $`T_c`$ which gives the usual critical exponents. In this sense it may be expected that natural boundaries are generic in second order transitions and that it is only the very special symmetries of integrablity which turn the curve of zeroes into a branch cut in the case of integrable models.

For first order transitions such as hard spheres at the freezing density the situation is much less clear. First of all it may be debated whether a freezing transition in which the system undergoes a change in symmetry is the same as the vapor/fluid condensation where there is no symmetry change. Even for the liquid/vapor transition the only models which have ever been solved are mean field models with long range attractive forces such as the Kac potential model of condensation . The pressure of these models is in fact analytic at the coexistence curve. However, for the Ising model at sufficiently low temperatures the free energy is not analytic at zero magnetic field  which in the language of the lattice gas means that there in a singularity at the coexistence boundary. This singularity is of the infinitely differentiable type discussed by Fisher  and Langer  but the additional prediction that the singularity is an isolated essential singularity has not been addressed.

Furthermore no model in a finite number of dimensions with only repulsive forces has ever been solved exactly which exhibits a freezing transition.

There is thus no compelling evidence to support the assumption that the scenario of case 2 must be correct for hard spheres and thus it is possible to raise the suggestion that case three holds for first as well as second order transitions. The difference being that for first order transitions the behavior at the coexistence density is presumably infinitely differentiable (but of course not analytic) instead of having an algebraic singularity.

This conjecture can be proven wrong if it can be demonstrated that there is a lower bound on the radius of convergence of the virial expansion which is greater than the fluid coexistence density of $`\eta _f=0.494.`$ However, the existing bound Lebowitz and Penrose  $`\eta _{pl}=.01809,`$ which comes from Groeneveld’s bound  on the fugacity expansion of the pressure in the grand canonical ensemble , is very far less than this density so the conjecture cannot be ruled out by the existing bounds.

It is of course widely appreciated that the bound $`\eta _{lp}`$ is much less than $`\eta _f.`$ However, because of the excellent agreement of the Padé approximate extrapolation of the first eight virial coefficients with the numerical simulations it is often said that “$`\mathrm{}`$ it (the lower bound of Lebowitz and Penrose) seems to be far below the true radius of convergence” . Indeed the nearest singularity in the Padè approximations is almost exactly at $`\eta _{cp}`$ and no sign of any kind is seen of the freezing transition in the eight term virial expansion data of ref. .

It is exactly at this point that the intuition obtained from the natural boundary in the Ising model is extremely useful. In ref.  the existence in the susceptibility of a dense set of singularities lying on the curve of zeroes of the partition function is demonstrated. But in refs. - where an expansion of the susceptibility is made to 323 terms only a small number (up to 8) of these dense set of singularities can be resolved out of the series expansion. We therefore conclude that unlike the power law singularities of second order transitions which are very strong and can be located with series as short as 8 or 10 terms that the singularities coming from natural boundaries require series of hundreds or thousands of terms to see. It is thus completely plausible and possible that the Padé approximate to the 8 term virial expansion is incapable of locating a natural boundary singularity at $`\eta _f`$ even if it does exist. Indeed it may be argued that since the kissing number (maximum number of spheres which can touch a given sphere) is 12 in three dimensions that the virial expansion cannot possibly include the effects of the geometry of hard spheres until at least the 12th coefficient has been computed. This is reminiscent of the early studies  on the Ising susceptibility where, based on the first 8 terms of the series expansion, a simple algebraic expression was conjectured which in fact omits almost all of the interesting structure of the true answer.

To illustrate the difficulty in seeing natural boundaries in a virial expansion consider the familiar example of the lacunary series

$$\underset{k=1}{\overset{\mathrm{}}{}}a_kz^{n_k}\mathrm{with}\mathrm{lim}_k\mathrm{}k/n_k=0,$$
(1)
which is well known to have a natural boundary. Typical examples of $`n_k`$ are $`k^2`$ and $`2^k.`$ The rapid growth of the powers of $`z`$ indicate that the natural boundary will be invisible in a Pade analysis of terms up to order $`z^8.`$ Note that this example exhibits the general feature of natural boundaries is that there is a sense in which the coefficients in the power series have oscillations which in a Padé analysis will lead to the appearance of more and more singularities in the complex plane as the number of terms in the series increases. We also note that if the coefficients $`a_k`$ fall off sufficiently rapidly (say as $`1/k!`$ for $`n_k=k^2`$ or $`2^k`$) that the singularity at $`z=1`$ will be infinitely differentiable. This is the behavior found by Isakov for the Ising model  and we therefore conclude that this result is just as compatible with a natural boundary as it is with the essential singularity of Fisher  and Langer .

It has, of course, always been possible to assert that the convergence of the virial expansion is determined by terms which have not yet been seen and without some evidence to the contrary it is just as reasonable to ignore this possibility. What has changed with the discovery of the natural boundary in the Ising model susceptibility is that not only do we now have an example where a natural boundary does occur but the reason which it occurs seems to be connected with the fact that at $`H0`$ the Ising model is not integrable. Therefore it would seem that natural boundaries are features of generic models (such as hard spheres) instead of the very specialized models such as the Ising model at $`H=0`$

Finally it is perhaps useful to comment on the question of possible sign changes in the virial coefficients. The mechanism of natural boundary formation can happen if all the virial coefficients are positive as the example of the lacunary series (1) shows. However it has been suggested as long ago as 1957 by Temperley  that “ the possibility of violent oscillations \[in the signs of the virial coefficients\] cannot altogether be ruled out.” Furthermore Ree and Hoover  showed in dimension 8 and 9 that that the fourth virial coefficient is indeed negative but they speculated that in three dimensions the signs would not change before the ninth or tenth coefficient. Indeed the first eight coefficients for hard spheres and discs  are now known to be positive as is the Pade estimate for the ninth coefficient. Nevertheless Temperley’s suggestion of violent oscillations in sign is an appealing mechanism because random oscillations in sign almost always lead to natural boundaries.

## 3 Conclusions

The possibility of a natural boundary in the pressure as a function of density as an explanation of the freezing transition of hard spheres is unorthodox but it does in fact have a certain intrinsic charm. Not only does it put first and second order transitions in some sense on an equal footing but it gives some hope that there is some predictability to first order transitions. If the mechanism of first order transitions is in fact that there are 2 (or more) branches of the free energy which do not know about each other then it is extremely difficult to ever say if you have in fact complexly computed the phase diagram. There may always be some other phase with a lower free energy which you just never thought of. In this situation one will always be in the position described by ref.  that “theory can only explain the transition afterwards”. In the end this reduces statistical mechanics to a descriptive rather than predictive science. Indeed this may account for the fact that many physicists prefer to study second order transitions and are happy to leave first order transitions to chemists and biologists. If, on the other hand, the natural boundary mechanism proves to be correct then the existence of the phase boundary can at least in principle be determined by studying the properties of the fluid phase alone. In many ways this is a much more attractive alternative. At the very least sufficient analysis should be done on the hard sphere system to demonstrate that it does not happen.

Acknowledgments

The author wishes to thank A. Guttmann, W. Orrick and J.H.H. Perk for discussions about their work on the Ising susceptibility. He also is indebted to R. Kamien , T. Lubensky and A. Sokal for extremely illuminating discussions about hard spheres. This work is supported in part by NSF grant DMR0073058.


text:
# Turbulent self-organized criticality

## Abstract

In the prototype sandpile model of self-organized criticality time series obtained by decomposing avalanches into waves of toppling show intermittent fluctuations. The $`q`$-th moments of wave size differences possess local multiscaling and global simple scaling regimes analogous to those holding for velocity structure functions in fluid turbulence. The correspondence involves identity of a basic scaling relation and of the form of relevant probability distributions. The sandpile provides a qualitative analog of many features of turbulent phenomena.

PACS numbers: 05.65.+b, 05.45.Tp, 45.70.Ht

Intermittency is a characterizing feature of turbulence . Intervals of low activity and of irregular duration are separated by sudden bursts in a wide range of strengths. Random quantities like activity variations across fixed space or time distances, are distributed according to non-Gaussian laws, whose moments display also a peculiar multiscaling. At first sight, some qualitative features of turbulent signals seem to be recognizable in the distant context of self-organized critical dynamics . In this field suitable models describe the bursts of activity by which a system, subjected to an external inflow of energy or matter, relaxes under the control of local, nonlinear threshold mechanisms. The bursts are called avalanches, and their time series show indeed some similarities to those of turbulent fluctuating quantities. Of course, these similarities are only apparent. Indeed, as a rule, different avalanches are just uncorrelated, and the probability distribution functions (pdf’s) of avalanche quantities do not show multiscaling.

Quite recently it was realized that the two-dimensional (2D) Bak-Tang-Wiesenfeld (BTW) sandpile, the prototype model of self organized criticality , obeys a peculiar form of multiscaling for the pdf’s of several avalanche measures . This result came after many years of work assuming allmost unanimously the simple finite size scaling form. The unexpected multiscaling makes the 2D BTW one special among models in the same class, and raises the issue of understanding its origin and the possible analogies with other complex scaling phenomena.

In this Letter we show that one can indeed establish a surprisingly close correspondence between BTW sandpile dynamics and fully developed turbulence. This correspondence includes an identity of the basic relations connecting the multiscaling exponents of the relevant pdf’s.

The BTW model is defined on a square lattice ($`L\times L`$) box . On site $`i`$ $`z_i=0,1,2,\mathrm{}`$ is the number of “grains”. If $`z_i<4`$, $`i`$, the configuration is stable. Grain addition to a stable configuration is made by selecting at random a site $`k`$ where $`z_kz_k+1`$. If then $`z_k4`$, toppling occurs, i.e. $`z_kz_k4`$, while each nearest neighbor, $`l`$, of site $`k`$ gets one grain ($`z_lz_l+1`$). If $`k`$ is at the border grains are dissipated. The toppling of site $`k`$ may cause instabilities in the neighbors, leading to further topplings at the next microscopic time step, and so on. Thus, an avalanche made by a total number $`s0`$ of topplings occurs before a new stable configuration is reached and a new grain is added. After many additions the system reaches a stationary critical state in which avalanche properties are sampled.

A key notion in the approach to the BTW is that of wave decomposition of avalanches . The first wave is obtained as the set of all topplings which can take place as long as the site of addition is prevented from a possible second toppling. The second wave is constituted by the topplings occurring after the second toppling of the addition site takes place and before a third one is allowed, and so on. The total number of topplings in an avalanche is the sum of those of all its waves. The global sample of all waves resulting from the decomposition of a sequence of several avalanches has exactly known scaling properties, well confirmed numerically . The wave size has a pdf satisfying finite size scaling $`P_w(s,L)s^{\tau _w}f_w(s/L^{D_w})`$, where, in $`d=2`$, $`\tau _w=1`$, $`D_w=2`$, and $`f_w`$ is a suitable scaling function. Most recently it has been shown that the time series of successive waves sizes $`\{s_m\}`$, $`m=1,2,..`$, has long range autocorrelation in the 2D BTW sandpile . This is not the case for other similar models, like the Manna one , in which the size of a wave is just uncorrelated with that of the preceeding one . The long time autocorrelation of waves is responsible for the fact that the scaling of the BTW avalanche size pdf, $`P_{av}`$, is different from that of $`P_w`$ . It remains a major challenge to identify which features of the BTW wave time series are at the root of the multifractal scaling pattern detected for avalanches . Here we find a key to this issue by establishing a correspondence between $`q`$-th order wave size difference moments of the sandpile, and velocity structure functions of a turbulent fluid flow.

In fully developed turbulence one measures a velocity field $`u(x)`$, with structure functions $`G_q(r)=|u(x+r)u(x)|^q^{1/q}`$ . $`G_q`$ displays a power law $`r`$-dependence $`G_qr^{\zeta _q}`$ in the inertial range $`1<r/\eta <\xi (l/\eta )R^{3/4}`$, where $`\xi `$ is a correlation length, $`R`$ is the Reynolds number, and $`l`$ and $`\eta `$ are the integral and the dissipation scales, respectively. The local dissipation rate $`ϵ(u/x)^2`$ has anomalous fluctuations. Indeed, the values of $`u/x`$ have a broad pdf, with moments obeying multiscaling in $`\xi `$. By matching the scaling of $`|u/x|`$ with that of $`|u(x+\eta )u(x)|`$ in the inertial range, one also expects $`G_q(r\eta )\xi ^{\alpha _q}`$, where $`\alpha _q`$ is another $`q`$-dependent exponent . On the other hand, as one considers $`u`$-increments for $`r/\eta \xi `$, the corresponding pdf’s, due to the lack of correlation, become $`r`$-independent and give information only on the global sample of velocity variations over all $`r`$, i.e. $`G_q(r)|uu|^q^{1/q}`$, for $`r/\eta \xi `$. The pdf’s of $`u(x+r)u(x)`$ become also Gaussian while multiscaling crosses over to a global simple scaling, so that $`G_q(r)\xi ^\alpha `$, for $`r/\eta \xi `$, independent of $`q`$ .

The scaling of $`G_q`$ is often detected by measuring the velocity field $`u`$ at different times at a fixed position in space . The data of this time series are then converted by interpreting time intervals as space intervals, once the average convection velocities are known. This suggests to regard the time series $`\{s_k\}`$ of BTW wave sizes in a similar perspective: also in this case we deal with a long range correlated signal, and it is legitimate to ask whether turbulence-like mechanisms in this signal are responsible for some form of multiscaling. Of course, a conversion of time into space dependence would not make sense in this case. Thus, in the analogy $`t`$ will definitely replace $`r`$. A structure function for the sandpile can be defined as

$$F_q(t)=|s_{k+t}s_k|^q^{1/q},$$
(1)
In turbulence the difference in velocities is a measure of how a perturbation of the average flow grows by propagating from one point to another. In sandpile dynamics the difference in size between two successive waves can also be considered as a measure of how a perturbation evolves in time. The instantaneous entity of the perturbation is represented by the number of topplings occurring in the corresponding wave. In analogy with turbulence, a scaling form one can expect for $`F_q`$ is :

$$F_q(t)\xi ^{\alpha _q}t^{\zeta _q}f_q(t/\xi ),$$
(2)
where $`f_q`$ are scaling functions, and $`\alpha _q`$ and $`\zeta _q`$ are exponents possibly varying with $`q`$, and specified here by the same symbols as the corresponding turbulence ones. $`\xi `$ in the BTW case is a correlation time for waves, which should scale $`L^z`$. It plays the role of the correlation length in turbulence, and is denoted by the same symbol here. The long time correlation is a peculiar feature of 2D BTW wave series, essential for establishing the correspondence with turbulence. Avalanche time series in general show finite correlation time for models of self organized criticality, including the 2D BTW sandpile .

If assumed for $`G_q`$, the form (2) embodies the various scaling regimes described above for turbulence . Indeed, we expect $`f_q(y\mathrm{})y^{\zeta _q}`$ in such a way that the $`t`$ dependence when $`t\xi `$ is correctly absorbed. In this limit one gets then $`F_q\xi ^{\alpha _q+\zeta _q}`$. We know further that, for $`t\xi `$, $`F_q`$ should reproduce the moment scaling of the global wave distribution $`P_w`$. Indeed, wave sizes are totally uncorrelated for $`t\xi `$, so that $`F_q(t\xi )|ss_{P_w}|^q_{P_w}^{1/q}`$, with $`_{P_w}P_w(s,L)ds`$. Since $`P_w`$ has a simple, constant gap ($`D_w=2`$) scaling form, $`F_q(t\xi )L^2`$. Thus, the exponents in Eq. (2) must satisfy the identity

$$\alpha _q+\zeta _q=2/z$$
(3)
for all $`q`$. A scaling relation analogous to 3 is expected to hold for the corresponding turbulence exponents . For $`1<t<\xi `$ one expects $`f_qconst.`$, which leads to a power law $`t`$-dependence, $`F_q(t)t^{\zeta _q}`$, in the analog of the inertial range. Finally, for $`t1`$, which corresponds to $`r\eta `$, one recovers the moments of the equivalent of $`|u/x|`$ in turbulence: $`F_q(1)\xi ^{\alpha _q}`$. The interval between successive waves, taken here as the unit of time, corresponds to $`\eta `$ in the turbulent case.

By analyzing samples of up to $`10^7`$ waves for $`L=128,256,512`$ and $`1024`$, we could verify that the scaling ansatz 2 and relation 3 are very well satisfied. Fig. 1 reports collapse plots of the quantity $`y^{\zeta _q}f_q(y)`$ for various $`q`$’s. The collapse quality increases with $`q`$ and $`z0.7`$ is an optimal choice for all cases. Since $`f_q`$ is approximately constant for $`y<1`$, from the low $`y`$ part of the plots one also obtains estimates of $`\zeta _q`$. The $`t`$-dependence of $`F_q`$ in the inertial range shows some degree of extended self-similarity, a well known property of structure functions in turbulence : log-log plots of $`F_q(t)^q`$ vs $`F_1(t)`$ (Fig. 2) are rather straight in the range corresponding to the inertial regime. In spite of the relatively narrow interval of variation of the moments in this range, this allows rather precise determinations of $`\zeta _q/\zeta _1`$. These, combined with an estimate of $`\zeta _1`$ from the plots in Fig. 1 lead to the values reported in Table I. The exponents $`\alpha _q`$ reported there are determined by log-log plots of $`F_q(1)`$ versus $`L`$, which are definitely linear and show clear $`q`$-dependent slopes. The values of $`z\zeta _q`$ and $`z\alpha _q`$ reported in Table I satisfy Eq. (3) very well.

A key fluctuating quantity in the description of turbulent flow is the dissipation rate averaged over spatial volumes of linear size $`r`$, $`ϵ_r`$ . In order to get further insight into the correspondence with turbulence, we studied an analogous wave quantity, given by the average over an interval $`t`$ of what should correspond to the local dissipation rate:

$$ϵ_t=\underset{k=1}{\overset{t}{}}(s_{k+1}s_k)^2/t$$
(4)
In analogy with turbulence we expect $`ϵ_t^q/ϵ_t^qt^{\mu _q}`$, in the inertial range. Fig 3 shows plots of this quantity which clearly manifest multiscaling and are remarkably similar to those one obtains for $`ϵ_r`$ in turbulence . The plots of course do not include the plateaus one finds in the turbulent dissipation range ($`r<\eta `$), which is not covered by our correspondence. Convergence to zero for all $`q`$ can also be perceived for large $`t`$. These results further support the identification of $`s_{k+1}s_k`$ as the analog of a velocity gradient in turbulence.

Thus, $`F_q`$ in the BTW sandpile has a scaling isomorphic to that of $`G_q`$ in turbulence. So far, theoretical approaches always tried to link avalanche to wave scaling . Since this last scaling is simple and exactly known at global level, attempts concentrated on the connection between avalanche scaling and this regime, just overlooking the very possibility of a different local regime and its consequences. We find that, in the range $`1<t<L^z`$, the wave series displays complex multiscaling structure, of the same type as that found for velocity fields in turbulence. This is not the case for the Manna sandpile , which does not show any form of multiscaling. In that model, due to the absence of correlation, the global regime of the structure function is reached as soon as $`t=1`$ and there is no window for local and inertial multiscaling regimes.

The intermittent local wave regime is certainly at the basis of the multiscaling already detected for global avalanche distributions of the 2D BTW model . An open problem remains that of connecting the multifractal exponents of the wave time series with those of avalanche distributions. Wave and avalanche time series are separated by an infinite coarse graining in time, for $`L\mathrm{}`$. At numerical level, one could in principle perform this coarse graining by partitioning waves into successive sets of $`n`$ elements . The corresponding block variables $`S_k=s_{n(k1)+1}+s_{n(k1)+2}+\mathrm{}+s_{nk}`$, $`k=1,2,\mathrm{}`$ are expected to have scaling properties approaching those of avalanche sizes in the limit $`n\mathrm{}`$. So, for coarse grained structure functions $`F_{q,n}=|S_{k+t}S_k|^q^{1/q}`$, the scaling of Eq. (2) should cross over, for increasing $`n`$, to the $`t`$-independent form $`F_{q,n}L^{\sigma _q/q}`$, where $`\sigma _q`$ is the multifractal moment exponent of the global avalanche size distribution $`P_{av}(s,L)`$, i.e. $`s^q_{P_{av}}L^{\sigma _q}`$ . Unfortunately, sampling limitations make such coarse graining extremely difficult to realize for large enough $`n`$.

The analogy with turbulence goes beyond the scaling properties alone. The pdf’s of velocity differences in turbulence show an increasing asymmetry between positive and negative increments as $`r`$ is reduced within the inertial range . This asymmetry, implying nonzero skewness, is an important feature of energy transfer across different scales. Fig. 4 reports a plot of the pdf of $`s_{k+1}s_k`$ for $`L=1024`$, which shows the expected slight asymmetry. One can also try to apply to such plots fitting schemes used in turbulence. We adopted a scheme proposed most recently, based on non–extensive statistics assumptions, and extremely successfull in reproducing fluid turbulence experimental data . Omitting details, to be published elsewhere, we just include in Fig. 4 the obtained best fit. This is very satisfactory within a substantial deviation range and the two parameters involved obey within less than a percent a relation established in Ref. .

Summarizing, BTW wave time series produce multiscaling schemes and pdf’s isomorphic to those of fully developed turbulence. This identifies the origin of the peculiar scaling of the 2D BTW sandpile, compared to similar models, and gives it a novel theoretical valence. Further elucidating why BTW waves behave like velocities of a turbulent flow, or why $`(s_{k+1}s_k)^2`$ can be viewed as a fluid energy dissipation rate, should provide deeper insight into this problematic model. Multiscalings consistent with Eqs. (2-3) have been previously identified also for a model of epitaxial growth in 1D  and, most recently, for the directed polymer in tilted columnar disorder in 2D . In the second case careful determinations of $`\alpha _q`$ and $`\zeta _q`$ revealed that the multiscaling is in fact just a bi-scaling, of the kind expected for Burgers turbulence , rather than for turbulence . Our results seem to definitely exclude bi-scaling in favor of genuine turbulent multiscaling. Thus, the analogies with turbulence discovered here include this crucial aspect. A possible connection between turbulence and self organized criticality was first conjectured when this paradigm was introduced . Our results show that indeed there exists a precise analogy in the case of the BTW sandpile, and that this analogy can be established through the concept of waves, which are in turn peculiar to a restricted class of models. However, this does not exclude that similar links with turbulence could be established on a different basis  or also for other models, in which for example waves are not defined .

Recent studies of the statistics of solar flares have correctly emphasized substantial differences existing between self-organized criticality and turbulence mechanisms, for the explanation of power law pdf’s . Such differences, which favor turbulence models, can ultimately be ascribed to the lack of long time correlations for the avalanches in the self organized critical models used so far to describe these phenomena. Quite remarkably, the distinction between the two mechanisms is however lost in the 2D BTW model, which has long time correlated bursts and genuine turbulence, if studied at the wave time scale. Self-organized critical models in which activity bursts (no matter whether waves or avalanches ) have long time correlations, would certainly be better candidates to describe flares. Of course, these correlations alone would not be necessarily enough to produce also turbulent scaling patterns.

Partial support from the European Network No. ERBFMRXCT980183 is acknowledged.