text:
# Untitled Document

A Calculator for Confidence Intervals

Roger Barlow

Department of Physics

Manchester University

England

Abstract

A calculator program has been written to give confidence intervals on branching ratios for rare decay modes (or similar quantities) calculated from the number of events observed, the acceptance factor, the background estimate and the associated errors. Results from different experiments (or different channels from the same experiment) can be combined. The calculator is available in http://www.slac.stanford.edu/~barlow/limits.html PACS code 13.85.Rm Keywords: Limits. Systematic errors. Combination of results. Rare processes.

1. Introduction

Many particle physics experiments measure a number $`n`$ of events of a particular type, and use it to set limits on a physical quantity $`R`$, typically a cross section or branching ratio. $`R`$ and $`n`$ are related by various factors such as the efficiency, dead-time, backgrounds, beam luminosity, and duration of the experiment.

For a Poisson process, the production of confidence limits on the mean $`\mu `$ from the observed $`n`$ has been well studied . Tables exist, or standard numerical techniques can be used, to solve exactly the desired probability equalities

$$\underset{r=0}{\overset{n}{}}P(r;\mu _+)=\alpha \underset{r=0}{\overset{n1}{}}P(r;\mu _{})=1\alpha $$
$`(1)`$
for the upper limit $`\mu _+`$ and the lower limit $`\mu _{}`$ at the $`(1\alpha )`$ confidence level. If the factor for efficiency etc., collectively known as the sensitivity $`S`$, is known exactly, then in the absence of background these limits translate directly into limits on $`R`$ using

$$\mu =RS$$
$`(2).`$
If the expected background $`b`$ is known exactly then it can be subtracted and the confidence limits on $`R`$ are given by

$$R_\pm =\frac{1}{S}(\mu _\pm b)$$
$`(3).`$
This note considers cases where there are uncertainties $`\sigma _S`$ and $`\sigma _b`$ in $`S`$ and $`b`$. It follows the approach of Cousins and Highland . Exact calculation is not possible, but a Monte Carlo technique can readily be used. A trial value for the signal is taken, then repeatedly it has a background mean drawn from a Gaussian of mean $`b`$ and standard deviation $`\sigma _b`$ added to it, is multiplied by a Sensitivity drawn from a Gaussian of mean $`S`$ and standard deviation $`\sigma _S`$, and is then used as the mean of a Poisson distribution from which a number of events $`r`$ is generated. By counting the fraction of cases in which $`r`$ is less than (or equal to) the observed number $`n`$, the probabilities $`1\alpha `$ (or $`\alpha `$) of Equation 1 are given.

2. Philosophy

This is basically a conventional frequentist limit. There are no prior assumptions about $`R`$. However uncertainties in $`S`$ and $`b`$ introduce a Bayesian viewpoint if (as they often do) they include factors like the accuracy of a theoretical model. Such errors are often admittedly not objectively defined: the experimenter makes statements like ‘we believe this is good to 10%’. Then the statement “$`S=0.50\pm 0.05`$”, if it is taken to mean that the sensitivity has a 68% probability of lying between 0.45 and 0.55, is ‘unscientific’ from a strict frequentist viewpoint .

This frequentist/Bayesian admixture is general and inevitable. It is defensible on the grounds that the main source of uncertainty – the Poisson statistics of the number of events – is treated in an objective frequentist fashion. The errors on the factors are small and have small consequences. The ambiguity due to prior distributions are merely one more uncertainty among many.

3. Results from differing priors

In interpreting a result one naturally writes

$$R=An$$
$`(4)`$
where the ‘appropriate’ factor $`A`$ is the inverse of the sensitivity $`S`$. This was done by the BaBar Statistics Working Group in their report for the collaboration .

It was found that probabilities obtained using the formalism of equation (2) and equation (4), smearing the factor $`S`$ or $`A`$ by the same relative amount, give slightly different results. Although one is merely the inverse of the other, the relationship is nontrivial. A Gaussian distribution for $`A`$ is not equivalent to a Gaussian distribution for $`S`$. The experimenter may have an a priori reason for believing that one or the other has a Gaussian uncertainty, however this is not the case in general and the choice is usually merely one of convenience.

A statistician would say  that the appropriate technique is to use an ‘invariant’ prior according to the suggestion of Jeffreys, i.e. to work with a uniform prior in the quantity for which the Fisher information is constant. For a scale factor $`A`$ this means a prior uniform in $`\mathrm{ln}A`$ or, equivalently, $`\mathrm{ln}S=\mathrm{ln}A`$.

For example, with no error on the sensitivity the probability of a signal 5.0 giving a result of 3 events or less is 26.5%. With increasing uncertainty the probability rises, but by slightly different amounts.

| % error | Cousins & | Jeffreys | BaBar |
| --- | --- | --- | --- |
| on $`S`$ | Highland |  | SWG |
| 0 | 0.265 | 0.265 | 0.265 |
| 10% | 0.272 | 0.269 | 0.266 |
| 20% | 0.291 | 0.277 | 0.267 |
| 30% | 0.316 | 0.291 | 0.273 |

Table 1: Probability for $`n3`$ with increasing uncertainty in the acceptance

Conversely, with 30% uncertainty the limit of 5.00 gives a 29.1% probability with a Jeffreys prior; to get the same probability with the Cousins and Highland hypothesis the limit has to be raised to 5.22, whereas for the BaBar SWG hypothesis it falls to 4.86. These differences are not large (and 30% is somewhat extreme) but this ambiguity does affect values at the level of precision with which they are generally quoted, so it cannot be ignored. All three values are given by the program, so the user can choose which to use (indeed, they are forced to be aware of the choice).

4. Combining results

When different experiments, or different channels in the same experiment, give limits, it is desirable to combine them. There is a possible ambiguity in the procedure. For, say, an upper limit from a single experiment one determines a value $`R_+`$ such that the probability of observing $`n`$ events or less is equal to the desired (small) level, as in Equation 1. But for a double channel, with results $`(n_1,n_2)`$, the expression ‘less’ is undefined: is $`(4,1)`$ less than $`(3,3)`$ ?

It could be suggested that ‘less than’ for a pair signify that both values should be less, or either value. However these are clearly unworkable: if the two experiments are identical (in luminosity, efficiency, and background) then there is clearly no information obtainable from the distribution between the two: $`n_1+n_2`$ is a sufficient statistic for $`(n_1,n_2)`$, and ‘less than’ can be defined with reference to the sum of the two results. Any definition of ‘less than’ must satisfy this.

We propose that ‘less than’ be defined for a multiple result using the value of the Branching ratio that one would deduce from it. This is given by a maximum likelihood estimator as (, Equation 129)

$$\underset{i}{}\frac{n_iS_i}{S_iR+b_i}S_i=0$$
$`(5)`$
where $`S_i`$ is the sensitivity factor for experiment $`i`$, and $`b_i`$ is the expected background. Sets of results can be ranked according to the value of $`R`$ given by this equation. Thus to run the simulation for several experiments one can first solve Equation (5) iteratively to find $`R_{data}`$ for the data values, and then for a given limit $`R_+`$ repeatedly generate random sets of result values $`\{n_i\}`$, calculate the resulting $`R`$, and count the fraction of cases that $`RR_{data}`$. In practice this can be achieved without an iterative solution for each result, as the quantity in Equation 5 is a monotonically decreasing function of $`R`$, so having found $`R_{data}`$ then for each set of results one merely calculates the quantity

$$Q=\underset{i}{}\frac{n_iS_i}{S_iR_{data}+b_i}S_i$$
$`(6)`$
and the sign of $`Q`$ gives the sign of $`RR_{data}`$.

This goes over to the conventional case if all experiments have the same $`S_i`$ and $`b_i`$.

5. Usage

The user is initally presented with a screen for a single experiment (Figure 1.) The number of observed events is entered, as are the expected background and its error (which, like all errors, is given as an absolute value), the overall sensitivity $`S`$ with its error, the number of Monte Carlo events to be used in the evaluation, and a guess for the limit. A button is then pressed to run the Monte Carlo simulation, and evaluate the probability for either an upper or a lower limit. (The difference is merely that one includes the $`r=n`$ result in the probability and the other does not, according to the two parts of Equation 1.)

Figure 1: Example of use with a single result

The probability for this limit, under the three different assumptions about priors, is displayed. The user can then iterate towards whatever limit type they choose, e.g. for a 90% upper limit they adjust the limit guess until the probability is shown as 0.10; for a 90% lower limit they would aim for 0.90. Only the ‘Limit guess’ data entry field need be modified for each trial, although as the limit is approached greater accuracy may be desired and the number of Monte Carlo events increased accordingly. (The calculation takes a few seconds for 1,000,000 events). In the example shown, 3 events were seen with zero background and no sensitivity uncertainty. With 10,000 Monte Carlo events the limit is quickly established at around 6.7. Increasing the number to 1,000,000 establishes it as 6.68 (which could have been obtained from tables anyway).

It is helpful to remember that the Poisson mean used is the product of the limit guess and the (smeared) sensitivity factor. It may be convenient to work in ‘event-count’ units, taking the sensitivity as unity ( with an error which is the the proportional error). The resulting limit is equivalent to a number of events and the sensitivity factor is applied afterwards: if in the above example the data sample corresponds to $`20fb^1`$ and the efficiency is 100% this establishes the limit as $`6.68÷20=0.334fb`$. Alternatively the Sensitivity factor could have been entered directly as 20, and the limit guess found to be 0.334.

If the efficiency were, say, $`50\pm 5\%`$ this could be accommodated either by setting the sensitivity error to 0.1 and dividing the resulting limit (risen to 6.81, using the Cousins and Highland prior) by a further 0.5, to give $`0.681fb`$, or by changing the sensitivity to 10 with an error of 1 and extracting a limit of 0.681 directly.

If the background is estimated as $`0.5\pm 0`$ then the limit falls from 6.81 to to 6.31. For a background of $`0.5\pm 0.2`$ it falls to 6.29.

If further experiments are desired (up to a limit of 10) there is a button to add data fields for them. There is now also a specific factor and error for each experiment, and this is entered in a further two columns, as shown in Figure 2.

Figure 2: Example of use with two results

For multiple results ‘event count’ units are not workable. The sensitivity is split into an overall value and a specific factor for each channel: the Poisson mean for an experiment is the product of the three values entered for the limit, the smeared overall sensitivity, and the smeared specific sensitivity.

The ambiguity in the factors (if the overall sensitivity is scaled up and each channel’s fraction scaled down by the same amount this makes no difference) is there for convenience in the handling of errors. The common acceptance is varied together for all channels, while the specific factors vary independently. For example, if a single experiment is measuring the rate of $`K_s^0`$ production, using both the $`\pi ^+\pi ^{}`$ and $`\pi ^0\pi ^0`$ decay modes, then the overall luminosity figure and global trigger efficiency (with associated errors) form the common value. The two specific factors are formed from the two branching ratios and the individual reconstruction efficiencies for making $`K^0`$ particles from tracks and from clusters. On the other hand, if two separate experiments were measuring the rate of $`Z`$ production through a decay mode to $`\mu ^+\mu ^{}`$, the common acceptance would contain the branching ratio $`Z\mu ^+\mu ^{}`$ (and associated error) and the specific factors contain the luminosities etc of the two experiments.

6. Application to a related problem

A different problem arises where an overall limit is required from two (or more) channels for which the efficiencies are known, but the branching ratio is not. For example, suppose the $`X`$ particle is hypothesised to decay to $`\mu ^+\mu ^{}`$ and $`e^+e^{}`$ in an unknown ratio, but with (different) efficiencies for detecting $`X\mu ^+\mu ^{}`$ and $`Xe^+e^{}`$ which are credibly estimated from Monte Carlo. From the observed numbers $`N_\mu `$ and $`N_e`$ what can be said about the upper limit for $`X`$ production?

This has been surprisingly difficult to answer. However we can now do so, at least in principle. There are two parameters: $`N_X`$, the number of $`X`$ particles, and $`p`$, the probability of decaying into muons; the probability of decaying into electrons is $`1p`$. We use the fact that in the full definition of confidence levels, the two parts of Equation 1 are inequalities: for the upper limit $`N_+`$ on $`N_X`$ the probability of obtaining this result, or worse, is $`\alpha `$, or less, whatever the value of $`p`$. For a given $`N_X`$ the probability must be maximised by varying $`p`$,

$$P(N_X)=Max(P(N_X,p))$$

and then $`N_X`$ can be varied until $`P(N_X)=\alpha `$. It may be that there is a simple way of finding the value of $`p`$ which maximises $`P(N_X,p)`$, though this is complicated as when $`p`$ changes not only does the probability for each $`(n_1,n_2)`$ value change, but the region which counts as ‘less than’ the data result changes discontinuously. However in default of a solution, the limit can be found using this calculator.

7. Possible future developments

It would be possible to relieve the user of the need to guess limits by doing the iteration within the program. This would require a more complicated front interface (specifying the confidence level required and the interval type: upper, lower, central, shortest-distance…) and an automatic scanning strategy to locate the limit to some desired accuracy.

The calculator does not consider the problems that arise when the signal $`n`$ is similar to, or smaller than, the background $`b`$, i.e. when there has manifestly been a fluctuation downwards in the background process. This cannot be taken into account in the standard frequentist framework. For such cases extension of the calculation to use either a flat-prior Bayesian formula formula  or the frequentist Feldman-Cousins technique  would be possible.

The expected background contributions are specified in terms of numbers of observed events. The model assumes that the error on this quantity is independent of the errors on acceptance. This corresponds to cases where backgrounds are estimated, for example, from observed sidebands. If one has a background from a known physics process affected by a factor in common with the signal (e.g. background from Drell-Yan muon pairs and a serious uncertainty in the muon detection efficiency) then the model could be adapted to do this.

8. References

 E. L. Crow and R. S. Gardner, Biometrika 46 441 (1959)

 A. Stuart and J. K. Ord, Kendall’s Advanced Theory of Statistics 5th Edition, Volume 2, Oxford University Press (1991)

 R. D. Cousins and V. L. Highland, Nucl. Instr & Meth. A320 331 (1992)

 R. von Mises, Probability, Statistics and Truth, Dover Publications (1981)

 The Babar Statistics Working Group, Recommended Statistical Procedures for Babar, BaBar Analysis Document # 318. (Version of October 18, 2001) http://www.slac.stanford.edu/BFROOT/www/Statistics/Report/report.ps

 D. R. Cox and D. V. Hinkley Theoretical Statistics, p 378, Chapman and Hall,

London (1974)

D. S. Sivia Data Analysis: a Bayesian Tutorial, p 113, Clarendon Press,

Oxford(1996)

 G. Cowan, Journal of Physics G: Nucl. Part. Phys. 27, p 1375 (2001)

 O. Helene, Nucl. Instr & Meth. A212 319 (1983)

 G. J. Feldman and R. D. Cousins, Phys. Rev. D57 3873 (1998)


format: text

text:
# Ultra-efficient Cooling in Ferromagnet-Superconductor Microrefrigerators

## Abstract

A promising scheme for electron microrefrigeration based on ferromagnet-superconductor contacts is presented. In this setup, cooling power densities up to 600 nW/$`\mu `$m2 can be achieved leading to electronic temperature reductions largely exceeding those obtained with existing superconductor-normal metal tunnel contacts. Half-metallic CrO2/Al bilayers are indicated as ideal candidates for the implementation of the device.

Heat transport through metal-superconductor interfaces can be applied to microcooling . The physical mechanism underlying this electronic cooling is quite simple. When a normal metal (N) is brought in contact with a superconductor (S), quasi-particle transport is effective only at energies larger than the S gap ($`>\mathrm{\Delta }`$). In fact, owing to the existence of the gap in the energy spectrum of the superconductor, at a given bias across the SN system only electrons possessing energies higher than the gap match available single-particle states and can transfer into the S portion of the junction (see Fig. 1(a)). This selective transfer of ”hot” carriers leads to the lowering of the effective electron temperature of the N electrode, even in the regime when electrons are thermally decoupled from the lattice. This situation can be experimentally realized with SN tunnel junctions where transport is dominated by quasi-particle dynamics. This unique property of Superconductor-Insulator-Normal metal (SIN) contacts was successfully employed for the realization of microcoolers . Intrinsically, however, SIN devices present large values of the contact resistance ($`R_n`$) which hinder carrier transfer and lead to a severe limitation in the achievable cooling powers .

Decreasing $`R_n`$ by using high-transparency interfaces is not a viable route to increase cooling power. In this case, in fact, another process dominates carrier transfer across the junction. At low $`R_n`$ values, carrier transfer is made possible by Andreev reflection (AR) even at energies $`<\mathrm{\Delta }`$ . In AR electrons may coherently evolve into holes that retrace the incoming electron time-reversed path and transfer a Cooper pair in S (see Fig. 1(b)). This does lead to increased conductivity, but does not contribute to thermal transport through the system. In practical devices a compromise must be reached and an optimal transparency value can be determined for a given microrefrigerator structure .

In this Letter we propose a cooling mechanism alternative to traditional SIN junctions which will be shown to be highly efficient. Our proposal is based on the use of a thin ferromagnetic layer (F) in good electric contact with S as shown in Fig.1(c). The physical basis of the FS microcooler operation resides in the spin-band splitting characteristic of F. In fact the electron (e) involved in AR and its phase-matched hole (h) must belong to opposite spin bands (Fig.1(d)); thus, suppression of the Andreev current  occurs in an FS junction and its intensity depends on the degree of the F-layer spin polarization $`𝒫`$. In the limit of large $`𝒫`$ and good metallic contact between the F and S electrodes, we observe a drastic suppression of the subgap current, while keeping efficient hot-carrier transfer leading to a considerable thermal flux $`𝒥`$. Cooling power densities are almost two orders of magnitude larger than those of optimized SIN junctions. This efficiency enhancement translates in a dramatic reduction of the final achievable electron temperature.

The device is sketched in Fig. 1(c). It consists of a NFS microcooler biased at a voltage $`V`$. Since the refrigeration stems from the FS portion, we shall focus on this junction only. Let us consider a ballistic FS contact with a perfectly transparent interface . In such a junction the overall cooling power $`𝒥`$ is equal to the sum of terms of the form

$$\frac{1}{h}𝑑\left\{\left[1BA\right]eV\left[1B+A\right]\right\}\left(f_Ff_S\right)$$
(1)
each relative to a given angle of incidence and spin species. In (1), $`f_{F,S}()`$ is the electron Fermi distribution function in F(S), whereas $`A()`$ and $`B()`$ are the Andreev and normal reflection probabilities obtained, in close analogy to Refs., as solutions of the Bogoliubov-de Gennes equations. The S pairing potential is approximated as a step-function , while F is described with the usual Stoner model . Differences in the effective-mass and Fermi-velocity values between F and S are neglected. Note that $`𝒥`$ is a function of spin polarization (via the scattering probabilities), bias voltage and temperature.

In order to evaluate the final electron temperature we need to consider also the mechanisms that can transfer energy into the ferromagnet. The main contribution is due to electron-phonon coupling and is given by  $`𝒥_l=\mathrm{\Sigma }𝒱(T_e^5T_l^5)`$, where $`𝒱`$ is the volume of the F electrode, $`T_{e(l)}`$ is the electron (lattice) temperature and $`\mathrm{\Sigma }`$ is a material-dependent parameter of the order $`10^9`$ WK-5 $`\mu `$m-3. For our purposes we can neglect the power dissipated in the resistance of the F layer . The final electron temperature $`T_e`$ is determined by the energy-balance equation $`𝒥(𝒫,V,T_e,T_l)+𝒥_l(T_e,T_l)=0`$. In writing this expression we set the superconductor temperature equal to the bath temperature $`T_l`$. Note that this is an idealized assumption since the presence of hot quasi-particles in S, that may strongly decrease the real cooling capability of the device, is expected. In real applications, however, one can overcome this problem through the exploitation of “quasi-particle traps” attached to the superconductor itself as described in Ref. .

Figure 2(a) shows $`T_e`$ as a function of bias for two starting lattice temperatures (i.e., the system temperatures at $`V=0`$) $`T_l=200`$ mK (solid lines) and $`T_l=300`$ mK (dotted lines). The calculation was performed for half-metallic ($`𝒫=1`$) FS (thick lines) and for SIN (thin lines) microcoolers, assuming Al as superconductor ($`\mathrm{\Delta }_{Al}=180\mu eV`$). Other junctions parameters were taken according to experimental values . In particular for the SIN-junction specific contact resistance we assumed a value $`0.4`$ k$`\mathrm{\Omega }`$$`\mu `$m2, corresponding to very high-quality Cu/Al2O3/Al junctions . Figure 2(a) shows the remarkable $`T_e`$ reduction provided by the FS cooler with respect to the NIS cooler. What shown in Fig. 2(a) is not meant to be a comparison with the best optimized SIN cooler (i.e., a SINIS refrigerator), however it helps in establishing the main differences between the single FS cooling element and the single SIN one. Even starting from 200 mK, for which the tunnel junction provides its largest cooling effect (a temperature reduction of about 10% at $`eV\mathrm{\Delta }`$), the FS cooler yields $`T_e`$ of the order of 10 mK (a temperature reduction of about 95%). This marked difference stems from the high contact resistance of the SIN junction that strongly affects its performance and must be compared to specific contact resistances as low as $`10^3\mathrm{\Omega }\mu `$m2 that are currently achieved in highly transmissive FS junctions . For $`eV<\mathrm{\Delta }`$, the FS-junction cooling power is large and, contrary to the SIN device, not localized only at $`eV\mathrm{\Delta }`$. For $`eV>\mathrm{\Delta }`$, instead, the increase in the electron temperature is quite rapid owing to the high current driven through the FS contact.

The large cooling-power surface density $`𝒥_𝒜`$ makes the FS microrefrigerator a high performance device. Figure 2(b) shows $`𝒥_𝒜(𝒫)`$ calculated at each optimal bias voltage, for some values of the lattice temperature. In the inset, the same quantity is plotted for a SIN junction versus interface transmission probability $`𝒟`$ . The comparison shows that for $`𝒫=1`$ at $`T0.4\mathrm{\Delta }/k_B`$ it is possible to achieve power densities up to 600 nW/$`\mu `$m2, i.e. 30 times larger than those for the SIN junction at the optimized transmissivity ($`𝒟3\times 10^2`$ at $`T0.3\mathrm{\Delta }/k_B`$). In real applications, however, SIN interface transmissivity is much smaller ($`𝒟10^6÷10^5`$) , limiting the achievable cooling power density to some pW/$`\mu `$m2.

The proposed cooling mechanism is robust against various possible sources that can lower quasi-particle transmission or enhance AR. We have evaluated the impact of incomplete polarization of the ferromagnet ($`𝒫<1`$), of low transmittance FS contacts and of spin-flip processes. Figure 3 shows the calculated dimensionless cooling power $`𝒥(V)`$ for some $`𝒫`$ values at $`T=0.4\mathrm{\Delta }/k_B`$. For each value of polarization there exists an optimal bias voltage which maximizes $`𝒥`$. For $`𝒫=1`$, $`𝒥(V)`$ is maximized around $`V\mathrm{\Delta }/e`$. Notably, even for $`𝒫=94\%`$ there still is a positive heat current. The inset of Fig. 3 shows $`𝒥(T)`$ calculated at the optimum bias voltage for some values of $`𝒫`$. The heat current is maximum for $`T0.4\mathrm{\Delta }/k_B`$. In addition, lowering $`𝒫`$ leads to a reduction of the useful temperature window for electron cooling. The same behavior was verified in the presence of spin-flip scattering events. We calculated this by allowing a random misalignment of magnetic moments within a given angle from the perfect alignment, which produces a decrease in the current polarization. (There is, however, no evidence that such processes occur in CrO2/S nano-junctions which are the system we suggest for the implementation of the FS microcooler.) An important figure of merit of the refrigerator is represented by its coefficient of performance (COP), which is the ratio between the cooling power and the total input power. Our calculations show that for $`𝒫=1`$ in the $`0.2÷0.4\mathrm{\Delta }/k_B`$ temperature range the COP, evaluated at the optimum bias voltage, exceeds 20$`\%`$ reaching its highest value (23$`\%`$) at about $`T0.3\mathrm{\Delta }/k_B`$. For $`𝒫=96\%`$ we obtain about $`10\%`$ efficiency at the same temperature.

In order to simulate a non-perfectly transparent FS contact caused by fabrication procedures or mismatch in material parameters , we have inserted an insulating barrier at the interface. The numerical simulation reveals that for $`𝒫=1`$ and low transmittance interface (e.g., $`𝒟=0.1`$), the variation in the final electron temperature is not appreciable.

The FS junction discussed so far is part of the NFS device of Fig. 1(c). We estimated the impact of the F-layer thickness on the NFS-device cooling power. We verified that cooling effects comparable to the FS case can be reached for thickness values of the order of a few nm (i.e. corresponding to the magnetic length $`\mathrm{}v_F/h_0`$, $`h_0`$ being the F exchange field and $`v_F`$ the Fermi velocity).

In conclusion we should like to indicate half-metallic CrO2 (chromium dioxide)  in combination with Al as the ideal candidate for the implementation of these devices. The well-developed CrO2 technology  can in fact be used to realize complex FS arrays coupled with normal electrodes for optimal device geometries.

We acknowledge a useful correspondence with J. P. Pekola. This work was supported by INFM under the PAIS projects EISS and TIN.


text:
# Stable vortex-antivortex molecules in mesoscopic superconducting triangles

\[

## Abstract

A thermodynamically stable vortex-antivortex pattern has been revealed in mesoscopic type I superconducting triangles, contrary to type II superconductors where similar patterns are unstable. The stable vortex-antivortex “molecule” appears due to the interplay between two factors: a repulsive vortex-antivortex interaction in type I superconductors and the vortex confinement in the triangle.

\]

Symmetrically-confined vortex matter in superconductors, superfluids and Bose-Einstein condensates offers unique possibilities to study the interplay between the $`C_{\mathrm{}}`$ symmetry of the magnetic field and the discrete symmetry of the boundary conditions. More specifically, superconductivity in mesoscopic equilateral triangles, squares etc. in the presence of a magnetic field nucleates by conserving the imposed symmetry ($`C_3`$, $`C_4`$) of the boundary conditions and the applied vorticity. As a result, in an equilateral triangle, for example, in an applied magnetic field $`H`$ generating two flux quanta, $`2\mathrm{\Phi }_0`$, superconductivity appears as the $`C_3`$-symmetric combination $`3\mathrm{\Phi }_0\mathrm{\Phi }_0`$ (further on denoted as “$`31`$”) of three vortices and one antivortex in the center. These symmetry-induced antivortices can be important not only for superconductors but also for symmetrically confined superfluids and Bose-Einstein condensates. Since the order parameter patterns reported in Refs.  have been obtained in the framework of the linearized Ginzburg-Landau (GL) theory, this approach is valid only close to the nucleation line $`T_c(H)`$. Can then these novel symmetry-induced vortex-antivortex patterns survive deep in the superconducting state? Several attempts have been already made to answer this crucial question. In the limit of an extreme type II superconductor ($`\kappa 1`$), it has been shown that for a thin-film square, a configuration of one antivortex in the center and four vortices on the diagonals of the square is unstable away from the phase boundary . According to the analysis based on the coupled nonlinear GL equations, the vortex-antivortex pairs are unstable and no antivortices appear spontaneously at the $`T,H`$ points far away from the $`T_c(H)`$ line. Possible scenarios of penetration of a vortex into a mesoscopic superconducting triangle with increasing magnetic field have been studied in Ref. . Two different states were considered: a single vortex state and a state in the form of a symmetric combination of three vortices and an antivortex with vorticity $`L_{av}=2`$ (“$`32`$” combination). The calculations have shown that while a single vortex enters the triangle through a midpoint of one side, the “$`32`$” combination turns out to be energetically favorable when the vortices are close to the center of the triangle. Equilibrium is achieved when a single vortex is in the center of the triangle. When approaching the phase boundary, the free energy of a single-vortex state tends to coincide with the free energy of the “$`32`$” combination, thus confirming conclusions  that formation of antivortices is possible in the close vicinity to the phase boundary.

The previous inferences  on vortex-antivortex states in mesoscopic structures seem to give us no hope to find stable vortex-antivortex configurations deeper in the superconducting state, considering them just as the features appearing in materials with $`\kappa 1`$ at the phase boundary together with superconductivity. Here we propose the new solution demonstrating the stability of the vortex-antivortex patterns. This solution is based on the simple conjecture made by one of the authors (VVM, ): the main source of the vortex-antivortex pattern instability, namely vortex-antivortex attraction, can be removed by taking – instead of type II – type I superconductors, where vortex-antivortex interaction becomes repulsive. Indeed, when passing through the dual point $`\kappa =1/\sqrt{2}`$, the vortex-vortex interaction changes the sign and becomes attractive at $`\kappa <1/\sqrt{2}`$. At the same time, the vortex-antivortex interaction becomes repulsive. Therefore, one can expect that presence of antivortices, together with confinement of vortices and antivortices due to a potential barrier at the boundaries, can stabilize novel vortex-antivortex patterns in a mesoscopic sample of type I superconductor. Optimizing the geometry and the sizes of mesoscopic samples, one can therefore fulfil the conditions necessary for the existence of stable vortex-antivortex configurations. For instance, presence of sharp corners is known to lead to a strongly inhomogeneous distribution of the superconducting order parameter in a mesoscopic sample. Enhanced superconducting condensate density at the corners prevents vortices from leaving the sample. Altogether, a triangular type I superconducting sample seems to be an appropriate candidate to search for a stable vortex-antivortex configuration.

Fig. 1 further illustrates this idea. In a type I superconductor with $`\xi >\lambda `$, vortex cores overlap first when vortices approach each other. This triggers vortex-vortex attraction (and vortex-antivortex repulsion). In contrast with that, in a type II superconductor with $`\xi <\lambda `$, local fields created by vortices overlap first when vortices approach each other, thus inducing vortex-vortex repulsion (and vortex-antivortex attraction). As a result, a vortex-antivortex combination, consisting of three vortices ($`L_{3v}=3`$) and one antivortex with vorticity $`L_{av}=1`$ (“$`31`$” combination, or $`3v+1av`$ molecule) is stable or unstable, respectively, in a mesoscopic type I or type II superconducting triangle.

To verify these intuitive considerations we investigate a mesoscopic equilateral triangle of a type I superconductor in an applied magnetic field. The side of the triangle is chosen ($`a=1\mu `$m) to be larger than both characteristic lengths of the GL theory, $`\xi `$ and $`\lambda `$. In order to provide a sample to be a type I superconductor, we should consider in fact a triangular prism with a height $`h\xi `$ (and also $`h\lambda `$). This allows us to avoid an increase of $`\kappa `$ resulting in a known effect when a mesoscopic sample made of material, which is type I superconductor in bulk, becomes effectively type II superconductor . For comparison, a mesoscopic triangle of type II superconductor will be also considered. In our calculations, we use the GL parameters of Pb (type I): $`\xi _{\text{Pb}}=82`$ nm, $`\lambda _{\text{Pb}}=39`$ nm, $`\kappa _{\text{Pb}}=0.48`$, and of Nb (type II): $`\xi _{\text{Nb}}=39`$ nm, $`\lambda _{\text{Nb}}=50`$ nm, $`\kappa _{\text{Nb}}=1.28`$ .

In the description of the superconducting properties of mesoscopic triangles we rely upon the GL equations for the order parameter $`\psi `$ and the vector potential A of the magnetic field $`\text{H}=\text{rot}\text{A}`$. In the dimensionless form, when keeping the temperature dependence explicitly, the GL equations are

$`\left(i\text{A}\right)^2\psi \psi \left[\left(1{\displaystyle \frac{T}{T_c}}\right)\psi ^2\right]=0,`$ (1)
$`\kappa ^2\mathrm{\Delta }\text{A}={\displaystyle \frac{i}{2}}\left(\psi ^{}\psi \psi \psi ^{}\right)+\text{A}\psi ^2.`$ (2)
Here the GL parameter $`\kappa =\lambda (T)/\xi (T)`$, and $`\xi (0)`$ serves as the unit length. The imposed boundary condition is

$$\text{n}\left(i\text{A}\right)\psi |_{boundary}=0.$$
(3)
Topological characteristics of solutions of the GL equations are determined by (anti)vortex core lines. One revolution along any closed path around such a line changes the phase of the order parameter by $`2\pi L`$, where $`L`$ is the winding number (vorticity) of a vortex or antivortex.

Close to the sample, the local magnetic field is distorted as compared to the applied magnetic field H. However, far away from the sample the distortion is negligible, and the symmetric gauge of the vector potential can be applied: $`\text{A}=\left[\text{H}\times \text{r}\right]/2.`$ The boundary condition, which equates the vector potential of the local field with that of the applied field, is justified at the boundaries of the simulation region, which is taken large enough to provide that all the changes of the magnetic field occur inside this region (see, for example, Refs. ). Further on, we choose the $`z`$-axis coinciding with the direction of the applied magnetic field, which is normal to the triangular base. Because $`h\lambda `$, a non-zero $`z`$-component of the magnetic field is appreciable in the vicinity of the bases, where the field skirts the prism partially penetrating it. Near the central cross-section in the $`xy`$-plane, the magnetic field and the order parameter are uniform in the $`z`$-direction, and the problem becomes effectively two-dimensional. The third dimension is taken into consideration by imposing that, inside the simluation region, the total magnetic flux through any cross-section, perpendicular to the $`z`$-axis, is the same. The GL equations (1), (2) with the above boundary conditions are solved numerically, using the finite-difference method, on a square mesh with the density of 200 nodes per side of the triangle.

The calculations are performed for temperature $`T/T_c=0.96`$. First, we find the magnetic field regions, in which states with antivortices are expected to take place in a type I superconducting triangle. The free energy calculations show that at low fields, from $`H_0=0`$ to $`0.14H_c(0)`$, where $`H_c(0)`$ is the thermodynamical critical field at zero temperature, no vortices penetrate the triangle. The order parameter distribution changes continuously from homogeneous at $`H_0=0`$ to a strongly inhomogeneous one, characterized by concentration of the superconducting phase in the corners. At $`H_00.14H_c(0)`$ one vortex enters the triangle. For vorticity $`L=1`$, the equilibrium is achieved when the vortex is in the center of the triangle. For magnetic fields within the region from $`H_0=0.25`$ to $`0.38H_c(0)`$, a giant vortex state with the total vorticity $`L=2`$ is energetically preferable. Further on, we focus on $`H_0=0.32H_c(0)`$, which corresponds to the states with the total vorticity $`L=2`$. This state can be represented by two possible configurations: (i) two vortices in the form of a multivortex or a giant vortex state; (ii) a symmetric combination of three vortices on the triangle bisectors and one antivortex in the center (“$`31`$” combination, or $`3v+1av`$ molecule). (Symmetric combinations with a larger number of vortices and antivortices such as “$`64`$”, “$`97`$” etc. possess a higher free energy than the “$`31`$” combination and are not considered.)

According to our calculations, it is the “$`31`$” combination that minimizes the free energy in case of a type I superconductor. In Fig. 2, the free energy for this combination is shown as a function of the distance $`d_v`$ counted from the center of the triangle along the bisectors to vortices. There are three minima of the free energy as a function of $`d_v`$. The first minimum, which is at $`1.29\xi (0)`$ from the center of the triangle, corresponds to a configuration when vortices are situated between the center of the triangle and the midpoints of the sides of the triangle. This is a saddle point for the free energy as a function of the coordinates $`(x,y)`$ in the plane of the triangle, and the state is unstable. The second minimum is reached when all the vortices are in the center of the triangle, and the vortex-antivortex combination degenerates to a giant vortex $`L_{2v}=2`$. This local minimum represents a metastable state.

The absolute minimum is reached when three vortices are situated between the center and the apexes of the triangle at $`1.83\xi (0)`$ from the antivortex in the center (Fig. 2). This vortex-antivortex molecule is thermodynamically stable. Its stability can be understood in the following way. The distribution of the squared modulus of the order parameter $`\psi (x,y)^2`$, which relates to the above stable vortex-antivortex molecule, is shown in Fig. 3a. Four zeros of $`\psi (x,y)^2`$ correspond to three vortices and one antivortex. The distribution of the magnetic field $`H(x,y)`$ consistent with the above pattern of $`\psi (x,y)^2`$ will be presented elswhere. The function $`\psi (x,y)^2`$ reaches its maximum value in the corners. These “islands” of the superconducting phase in the corners prohibit vortices, which are repelled by an antivortex in the center, from leaving the triangle through the corners. Thus, vortices, being confined in a type I superconducting triangle and interacting with an antivortex, form a stable vortex-antivortex molecule.

It is worth noting, that for our type I sample a strongly enhanced nucleation field appears due to the confinement of the superconducting condensate in the mesoscopic triangle. In fact, this provides the “soft” scenario for the nucleation of the order parameter, like in bulk type II superconductors.

The stable vortex-antivortex patterns are qualitatively different in case of a type II superconducting triangle. In Fig. 4, the free energy of the “$`31`$” combination is shown (curve “a”) as a function of the distance $`d_v`$ from the center of the triangle to vortices. The lowest minimum is reached when all the vortices are in the center of the triangle. This means that a giant vortex with vorticity $`L_{2v}=2`$ is energetically more favorable in a type II superconducting triangle. (The local minimum at $`d_v=11.57\xi (0)`$ represents an unstable state.) The corresponding distribution of the squared modulus of the order parameter $`\psi (x,y)^2`$ is plotted in Fig. 3b. Although a giant vortex state with $`L_{2v}=2`$ has a lower free energy than the “$`31`$” combination, the equilibrium is reached for a multivortex state when two vortices are at two different bisectors (cf. Ref. ) of the triangle (see Fig. 4, curve “b”). The order parameter pattern corresponding to this stable two-vortex state is shown in Fig. 3c.

In conclusion, we have found deep in the superconducting state a thermodynamically stable vortex-antivortex configuration for a mesoscopic type I superconducting triangle, although until now it has been thought that vortex-antivortex patterns are unstable and they can manifest themselves only in the close vicinity to the phase boundary. Vortex-antivortex arrays become unstable in a type II superconducting triangle, in accordance with previous reports. The stability of the vortex-antivortex molecules in type I superconducting triangles is due to the change of the sign in the vortex-vortex and vortex-antivortex interaction forces when passing through the dual point $`\kappa =1/\sqrt{2}`$, combined with the condensate confinement by the boundaries of the mesoscopic triangle.

This work has been supported by GOA BOF UA 2000, IUAP, the FWO-V projects Nos. G.0306.00, G.0274.01, WOG WO.025.99N (Belgium), and the ESF Programme VORTEX. Useful discussions with L. Van Look, M. Morelle and G. Teniers are acknowledged.


text:
# Scattering of Stars by Transient Spiral Waves

## 1. Introduction

Barbanis & Woltjer (1967) first showed that a spiral wave disturbance of finite time duration will, on average, scatter stars away from circular orbits. Many subsequent $`N`$-body simulations (e.g., Sellwood & Carlberg 1984) have confirmed that transient spiral waves do heat the disk. Earlier, Spitzer & Schwarzschild (1953) had shown that stars are also scattered by dense objects, such as giant molecular clouds, co-orbiting within the disk.

The theory of these scattering processes has been developed in many subsequent papers, and it is now generally agreed that: (1) transient spirals heat mostly the in-plane motions, while (2) giant molecular clouds do little heating themselves but act mostly as scatterers which try to isotropize peculiar velocities by converting some of the in-plane components to vertical.

Scattering by spiral waves is still not completely understood. The cumbersome epicyclic approach by Barbanis & Woltjer was supplanted by the more powerful action-angle machinery introduced by Lynden-Bell & Kalnajs (1972). Subsequent work by Dekker (1976), Carlberg & Sellwood (1985), and Fuchs (2001) has developed a formalism which is somewhat related to quasi-linear theory in plasma physics. Binney & Lacey (1988) derived diffusion coefficients in action space while Jenkins & Binney (1990) used solar neighbourhood data to determine the relative efficiency of scattering by clouds and by spirals. Little & Carlberg (1991) report a numerical study of the angular momentum changes of test particles under the influence of a transient, weak, imposed bar and report results in general agreement with the predictions of Lynden-Bell & Kalnajs.

The conclusions from all this work correctly predict orbit averaged changes caused by spirals: small net changes are expected at corotation, while larger non-zero net changes occur at Lindblad resonances. But it is not generally appreciated that typical changes at corotation due to a moderate amplitude transient spiral are large, and spirals can cause stars to migrate over substantial radial distances – in fact, several papers have concluded the opposite! Wielen (e.g. Wielen et al. 1996) is a notable exception to this generalization; he has long suspected that the Sun, because of its atypically high metallicity, has migrated some 2kpc outwards from its place of birth.

Here we report our on-going studies of scattering by spiral waves in $`N`$-body simulations, and describe the range of possible changes for individual particles. In particular, we find that transient spirals of moderate amplitude cause surprisingly large angular momentum changes, but no significant heating, in the vicinity of corotation.

## 2. Simulation

Figure 1 shows the evolution of a half-mass, $`V_0=`$const., or Mestel, disk which is gently tapered at its inner and outer edges. We set the parameter of the distribution function (e.g. Binney & Tremaine 1987; eq 4-162) to be $`q=11.44`$ in order that $`Q=1.5`$ in the untapered disk. Toomre (1981) predicts this disk is stable if smooth, but particle noise in simulations (which is also present in real galaxy disks!) always induces some spiral activity. The $`N=1`$M particles of the simulation are confined to a plane and their mutual gravitational forces are determined with the aid of a 2-D polar grid. Forces from azimuthal wavenumbers $`0m4`$ only are included, and are derived from the usual Plummer softening law, with scale length $`0.1R_0`$. We use $`V_0`$, the constant circular speed of the Mestel disk, as our velocity unit, and the mean radius of the inner taper $`R_0`$ as our length unit; thus our units of time and angular momentum are $`R_0/V_0`$ and $`R_0V_0`$ respectively.

Figure 2 shows the change in angular momenta from $`t=0`$ to $`t=800`$ of a representative sample of the particles. The changes are generally quite large (rms $`2.68`$), even though the spiral activity was very mild; non-axisymmetric variations in surface density exceed 10% for only two or three patterns. We will show later that the angled features are signatures of scattering at corotation. Closer analysis indicates that most particles are affected by more than one wave over any short period. Evidently, the changes shown in Figure 2 are far too complicated to understand in detail.

## 3. An “isolated” spiral

In order to have a chance to understand this behaviour, we turn to a much more restricted simulation which was deliberately crafted to excite an isolated spiral mode that stands out clearly from the noise. It uses the same $`Q=1.5`$, half-mass Mestel disk, but particle noise was suppressed by using a quiet start (Sellwood 1983), disturbance forces were restricted to $`m=2`$ only, and the single spiral pattern was seeded by a groove (Sellwood & Kahn 1991). The groove in angular momentum density is centered on $`L_z=6.5`$ and has a half-width in $`L_z`$ of 0.2.

The later part of the evolution of this simulation is shown in Figure 3. The single transient spiral pattern stands out well, although the evolution is, unfortunately, tainted at the end by a weak bar pattern that appears as the spiral decays. The bar feature seems unavoidable. Note that since disturbance forces include only the $`m=2`$ component, this simulation does not capture the full non-linear evolution.

### 3.1. Standard results for a steady wave

Jacobi’s invariant, $`E_J=E\mathrm{\Omega }_pL_z`$, is a conserved quantity for a test particle in a steadily rotating non-axisymmetric potential. Therefore, changes in angular momentum and energy are related as

$$\mathrm{\Delta }E=\mathrm{\Omega }_p\mathrm{\Delta }L_z.$$
(1)
It is helpful to draw the classical Lindblad diagram to understand scattering of particles by a steadily-rotating, non-axisymmetric perturbation. The solid curve in Figure 4 is the curve for circular orbits which separates the inaccessible region of the diagram (shaded) from the region where orbits are eccentric. The energy of random motion is the vertical distance from the circular orbit curve.

Because of equation (1), particles are scattered along lines of constant slope, $`\mathrm{\Omega }_p`$. In general, scattering along a line of slope $`\mathrm{\Omega }_p`$ leads to changes in the relative amounts of orbital to epicyclic energy for a particle, but because the slope is also the slope of tangent to the circular orbit curve at corotation, changes at corotation do not change the random energy of a particle to first order. Finally, there is a bias in the possible shifts, which must be away from circular orbits when not at corotation. This bias away from circular orbits is the root cause of heating by spiral waves away from corotation.

### 3.2. Effect of finite growth

The spiral wave in the simulation grows exponentially, saturates and decays, implying that the perturbing potential does not have a single pure frequency as assumed in equation (1). To see the effect of the finite duration of a growing wave, we approximate the perturbing potential as a purely growing mode

$$\mathrm{\Phi }_1(r,\varphi )e^{im\mathrm{\Omega }_pt}e^{\gamma t}\times H(t)[1H(tT)],$$
(2)
where the Heaviside functions limit the time interval, $`0<t<T`$, over which the perturbation is non-zero. Here, $`\mathrm{\Omega }_p`$ is the pattern speed of the $`m`$-armed mode, $`\gamma `$ is its growth rate, and $`m=2`$ in our case. The power spectrum of this perturbation is an expression with three terms, but once $`\gamma T1`$, the dominant term has a simple Lorentzian form; in this limit, the power spectrum, normalized to integrate to unity, is

$$\left|\mathrm{\Phi }_\mathrm{\Omega }(r,\varphi )\right|^2\frac{m}{\pi }\frac{\gamma }{m^2(\mathrm{\Omega }\mathrm{\Omega }_p)^2+\gamma ^2}.$$
(3)
Particles in the the simulation experience a spread of frequencies given by this expression, and thus the ratio $`\mathrm{\Delta }E/\mathrm{\Delta }L_z`$ measured from the particles should also have this spread. Accordingly, we determined $`\mathrm{\Delta }E`$ and $`\mathrm{\Delta }L_z`$ separately for every particle in the simulation between the start and $`t=300`$ and histogram their ratio in Figure 5; we have discarded particles with $`\mathrm{\Delta }L_z/L_z<0.01`$ for which $`\mathrm{\Delta }E/\mathrm{\Delta }L_z`$ is less precisely determined. The shaded histogram is for particles with initial $`L_z>2`$ while the unshaded is for all particles. Note that particles with $`L_z<2`$ are all inside the inner Lindblad resonance of the spiral.

The smooth curve in Figure 5 is equation (3) with the parameters $`\mathrm{\Omega }_p=0.141`$ and $`\gamma =0.047`$ measured from the simulation – i.e. the curve is not fitted to the histogram, but the curve and the shaded histogram are normalized to enclose the same integrated area. The correspondence between the prediction and the empirical values clearly shows that the large majority of particles with $`L_z>2`$ is affected only by the one pattern. The curve is not quite a perfect fit to the shaded histogram; some discrepancies are due to the fact that the wave saturated and decayed before $`t=300`$.

The scattering trajectories shown in Figure 4 assume a steadily rotating, constant amplitude pattern, but because the frequency is broadened by the finite growth time of the mode, particles will be scattered along lines having the range of slopes spread according to equation (3) about the mean, which has slope $`\mathrm{\Omega }_p`$.

### 3.3. Scattering by one wave

Figure 6 (left) shows the mean (line) and $`\pm 1\sigma `$ distribution of $`\mathrm{\Delta }L_z`$ as a function of initial $`L_z`$. Note that these changes are azimuthally averaged. Corotation and the two Lindblad resonances are marked. Most changes occur near corotation, but with huge scatter; most particles inside corotation gain $`L_z`$ while most outside lose $`L_z`$. There are roughly as many gainers as losers, while the average change for each particle is large enough that most cross this resonance; in effect, particles generally exchange places across corotation. Figure 6 (right) confirms that the large changes at corotation cause very little heating; significant heating occurs only in the inner disk, where it is very difficult to disentangle the contributions from the inner Lindblad resonance (ILR) and from the bar.

Figure 7 shows the same plots but now for the $`20\%`$ of particles that start with the smallest epicycle energies. For these, the changes around corotation are even larger, but now one can see a clear excess of losers between ILR and corotation and a small excess of gainers between corotation and outer Lindblad resonance (OLR). The panel on the right shows the increase in rms radial velocity of the same particles; some Lindblad resonance heating is noticeable, although this subset of the particles may not reveal the full story. Changes in both panels of Figure 7 are larger inside corotation than outside because the wave amplitude is appreciably greater inside than outside.

Note that the greatest changes occur on opposite sides of corotation for particles several times farther from the groove center than the groove half-width ($`0.2`$ in $`L_z`$). Thus, these changes are produced by the spiral response, and are not a feature caused by the narrow groove we used to seed the instability.

### 3.4. Azimuthal variations

Figure 8 shows $`\mathrm{\Delta }L_z`$ for individual particles as a function of their azimuthal angle at $`t=300`$. The angle, $`w_2`$, is the variable conjugate to the angular momentum (Lynden-Bell & Kalnajs 1972), and is determined from the instantaneous position of the particle assuming an undisturbed, axisymmetric potential. The four panels show different groups of particles whose initial angular momenta placed them just inside or outside corotation (top), near the ILR and the OLR (bottom). Only the $`44`$% of particles with smallest epicyclic energies are plotted; the patterns in these panels are less marked when all particles are plotted. The obvious double period is to be expected for a two-arm pattern.

The coherence of the changes with phase relative to the wave is striking. Particles are being moved in a systematic manner by the large-scale potential perturbation, to produce these striking ribbons in phase space. It is clear that a large part of the spread at corotation seen in Figure 6 is not random scatter, but is strongly correlated with azimuth.

## 4. Radial mixing of stars and gas

From these Figures, we conclude that (1) “Lindblad resonance” heating is minor and predominantly associated with small peculiar velocities, (2) wave damping begins well inside the resonances, and (3) that angular momentum changes at corotation dominate.

Evidently, the magnitude of the changes induced by a single spiral will depend strongly on its amplitude. That shown in Figure 3 reaches a peak overdensity of $`20\%`$ – 25%, making it about twice the intensity of the strongest features seen in Figure 1. Neither value seems unreasonably large for real galaxies, perhaps even rather conservative; if the typically larger intensity variations in K-band (e.g. Rix & Zaritsky 1995; Block & Puerari 1999) really scale as mass variations, then the angular momentum changes in real galaxies must exceed those reported here.

Assuming the time-dependence of this spiral to be typical of those in real galaxies, we find that a single spiral wave generally changes the angular momentum of a star near corotation by $`\pm 20\%`$. Many stars are simply moved across corotation to another orbit with quite different angular momentum, but very similar radial action. The changes near the Lindblad resonances, which cause most of the heating (Lynden-Bell & Kalnajs 1972), are quite minor.

The churning of stars by multiple spiral events (Figure 2) has important consequences for the chemical abundance of stars within disks: gradients, which were established at the time stars formed must become shallower over time as a result of this activity.

It is important to realize that gas clouds are also scattered by the spiral perturbation. To obtain some idea how the interstellar medium might respond to the spiral, we show in Figure 9 the evolution of rings of test particles, which began on initially circular orbits in the simulation shown in Figure 3. The particles have no mass, and therefore do not contribute to the potential, but are integrated forward in the evolving potential from the disk of massive particles. The rings illustrate how streams of cold gas would behave up until the time that rings either self-intersect, or cross another ring; from that moment, collisions (or pressure if the gas has a finite temperature) will intervene and alter the subsequent motion of gas. Even so, the radial extent acquired by some (initially circular) rings indicates that substantial radial flows of gas are expected with in-going and out-going streams largely avoiding each other until about the time that the spiral saturates.

## 5. Conclusions

The main effect of a transient spiral wave is to churn the angular momentum distribution, with stars mostly changing places across corotation. A succession of spiral waves of different pattern speeds causes a random walk in $`L_z`$ with a variable step size of rms $`20\%`$ of a star’s initial $`L_z`$. It is therefore perfectly possible that the Sun has migrated by 2 kpc from its birthplace, as argued on other grounds by Wielen et al. (1996).

Remarkably, such substantial changes are not associated with any increase in random motion near corotation, which is perhaps why they have not previously attracted much attention. They are also not really expected from previous work, such as the analyses presented by Lynden-Bell & Kalnajs (1972) and Carlberg & Sellwood (1985). Those analyses assumed small departures from the unperturbed orbits, were also averaged over all phases, and focused on the limit of a steady wave, for which lasting changes are localized to resonances. These assumptions lead to quite different expectations but, with hindsight, seem rather questionable.

The churning of the disk by spiral waves in this manner has extensive implications for chemical abundance gradients. If stars can migrate considerably from their places of birth, abundance gradients established when they were formed will be diluted. Furthermore, the gas in the disk will also be churned, and the coherent azimuthal dependence shown in Figure 9 shows that crossing of streams does not occur until considerable radial displacements have been achieved.

#### Acknowledgments.

We thank James Binney and Scott Tremaine for helpful discussions. This work was supported by NSF grant AST-0098282 to JAS and by grant SFRH/BD/3444/2000 to MP from Fundação para a Ciência e a Tecnologia, Portugal.

## References

Barbanis, B. & Woltjer, L. 1967, ApJ, 150, 461

Binney, J. & Lacey, C. 1988, MNRAS, 230, 597

Binney, J. & Tremaine, S. 1987, Galactic Dynamics (Princeton: Princeton University Press)

Block, D. L. & Puerari, I. 1999, A&A, 342, 627

Carlberg, R. G. & Sellwood, J. A. 1985, ApJ, 292, 79

Dekker, E. 1976, Phys. Rep., 24, 315

Fuchs, B., 2001, MNRAS, 325, 1637

Jenkins, A. & Binney, J. 1990, MNRAS, 245, 305

Little, B. & Carlberg, R. G. 1991, MNRAS, 251, 227

Lynden-Bell, D. & Kalnajs, A. J. 1972, MNRAS, 157, 1

Rix, H-W. & Zaritsky, D. 1995, ApJ, 447, 82

Sellwood, J. A. 1983, J. Comp. Phys., 50, 337

Sellwood, J. A. & Carlberg, R. G. 1984, ApJ, 282, 61

Sellwood, J. A. & Kahn, F. D. 1991, MNRAS, 250, 278

Spitzer, L. & Schwarzschild, M. 1953, ApJ, 118, 106

Toomre, A. 1981, in The Structure and Evolution of Normal Galaxies, ed. S. M. Fall & D. Lynden-Bell (Cambridge: Cambridge University Press), p. 111

Wielen, R., Fuchs, B. & Dettbarn, C. 1996, A&A, 314, 438


text:
# Anelastic spectroscopy study of the spin-glass and cluster spin-glass phases of La2-xSrxCuO4 (0.015<𝑥<0.03)

## I INTRODUCTION

The magnetic properties of the High $`T_c`$ Superconductors (HTCS) have attracted much interest also because these compounds are experimental realizations of a 2D Heisenberg antiferromagnet. In undoped La2CuO4 the Cu2+ spins order into a 3D antiferromagnetic (AF) state below $`T_N315`$K with the staggered magnetization in the $`ab`$ plane. When doping by substituting Sr for La, the long range AF order is rapidly suppressed and around $`x_c0.02`$ the Néel temperature drops to 0 K. In the long range AF region $`T_\text{N}\left(x\right)`$ follows a power law relationship with $`x`$, interpreted as an indication that the holes introduced by doping form walls separating domains of undoped material. Later work indicated that the spin degrees of freedom associated with the doped holes are distinct from the in-plane Cu2+ spin degrees of freedom that order themselves below $`T_\text{N}`$, and the localization of the doped holes allows the associated spins to progressively slow down and freeze. The state in which the doped spins freeze is usually referred to as a spin glass (SG) state, and occurs below $`T_f\left(x\right)\left(815\text{K}\right)x`$. For $`x>x_c`$ there is no long range AF order, but approaching $`T_g0.2`$ K$`/x`$ AF correlations develop within domains separated by charge walls, with the easy axes of the staggered magnetization uncorrelated between different clusters. This picture corresponds to a cluster spin-glass (CSG) state. The formation of the SG and CSG states are inferred from sharp maxima in the 139La NQR and $`\mu `$SR relaxation rates, which indicate the slowing of the AF fluctuations below the measuring frequency ($`10^710^8`$ Hz in those experiments), and from the observation of irreversibility, remnant magnetization, and scaling behavior in magnetic susceptibility experiments. Also anelastic spectroscopy can provide useful information about the magnetic properties of high HTCS. Above $`x_c`$, the elastic energy loss coefficient shows a rise below a temperature $`T_{on}`$ close to the $`T_g\left(x\right)`$ for freezing into the CSG state. The absorption is not peaked near $`T_g`$, but is step-like or at least displays a plateau, and therefore does not directly correspond to the peak in the dynamic spin susceptibility due to the spin freezing; rather, it has been attributed to the stress-induced changes of the sizes of the spin clusters, or equivalently to the motion of the domain walls.

Recently, Matsuda et al. reported a neutron scattering study of the magnetic correlations in La2-xSrxCuO4 for $`x<x_c`$, which suggests a different picture of the spin glass phase. In fact, they found that also at $`x<x_c`$ the 3D AF ordered phase coexists below $`30`$ K with domains of ”diagonal” stripe phase (with the hole stripes at $`45^\text{o}`$ with respect to the Cu-O bonds), as observed for $`x>x_c`$. According to these authors, the hole localization starting around 150 K involves an electronic phase separation into regions with $`x_10`$ and $`x_20.02`$. The volume fraction of the $`x_2=0.02`$ phase changes as a function of the Sr doping, in order to achieve the average $`x`$.

In the following we will report on anelastic spectroscopy measurements of lightly doped La2-xSrxCuO4, where the anelastic spectra for $`x<x_c`$ present the same features attributed to the domain wall motion into the CSG phase, although attenuated. The step-like rise of dissipation occurs near 10 K, as for $`x0.02`$, supporting the view of an electronic phase separation, at least for $`0.015<x<0.02`$.

## II EXPERIMENTAL AND RESULTS

The samples where prepared by standard solid state reaction as described in Ref.  and cut in bars approximately $`40\times 4\times 0.6`$ mm3. In the as-sintered state all the samples contained small amounts of interstitial O, which was outgassed by heating in vacuum up to $`790`$K. The complex Young’s modulus $`E\left(\omega \right)=E^{}+iE^{\prime \prime }`$, whose reciprocal is the elastic compliance $`S=E^1`$, was measured as a function of temperature by electrostatically exciting the flexural modes. The vibration amplitude was detected by a frequency modulation technique. The vibration frequency, $`\omega /2\pi `$, is proportional to $`\sqrt{E^{}}`$, while the elastic energy loss coefficient (or reciprocal of the mechanical $`Q`$) is given by $`Q^1(\omega ,T)=`$ $`E^{\prime \prime }/E^{}=`$ $`S^{\prime \prime }/S^{}`$, and was measured by the decay of the free oscillations or the width of the resonance peak. The imaginary part of the dynamic susceptibility $`S^{\prime \prime }`$ is related to the spectral density $`J_\epsilon (\omega ,T)=`$ $`𝑑te^{i\omega t}\epsilon \left(t\right)\epsilon \left(0\right)`$ of the macroscopic strain $`\epsilon `$ through the fluctuation-dissipation theorem, $`S^{\prime \prime }`$ $`\left(\omega /2k_\text{B}T\right)J_\epsilon `$. Magnetoelastic coupling can couple strain to the spin degrees of freedom, so that pseudodiffusive spin excitations (with null characteristic frequency) or the motion of magnetic boundaries can contribute to $`J_\epsilon `$. An elementary relaxation process with a thermally activated relaxation time $`\tau \left(T\right)`$ contributes with $`S^{\prime \prime }`$ $`\omega \tau /\left[1+\left(\omega \tau \right)^2\right]`$, peaked at the temperature at which $`\omega \tau =1`$ (which therefore increases with increasing $`\omega `$), while the relaxation of extended and interacting structures like domain walls generally gives rise to a broader dissipation curve.

The nominal compositions of the samples were: $`x_{\text{nom}}=`$ $`0.015`$, 0.016, 0.018, 0.024 and 0.030. The final Sr contents were checked from the temperature position of the step in the Young’s modulus $`E`$ due to the tetragonal (HTT) / orthorhombic (LTO) transition, which occurs at a temperature $`T_t\left(x\right)`$ decreasing with doping approximately as $`T_t\left(x\right)=\left(535x/0.217\right)`$ K (Ref. ). The shape of the anomaly in $`E`$ cannot be completely fitted by a simple model, since it includes both contributions from the coupling of the spontaneous strain with the soft mode and from the domain wall motion, but its position in temperature provides information about the level of doping. Figure 1 presents the logarithmic derivative of the Young’s modulus with respect to temperature, $`d\mathrm{ln}\left[E/E\left(0\right)\right]/dT`$, so that the steps become peaks whose widths provide upper limits to the spread of the local Sr concentration; the peaks are fitted with lorentzians. The Sr concentrations estimated from the peak positions and from $`T_t\left(x\right)`$ are $`x=0.0165`$, 0.0155, 0.0183, 0.0227 and 0.0273 respectively; the sample with $`x_{\text{nom}}=0.016`$ results to be less doped than the one with $`x_{\text{nom}}=0.015`$ and with a broader transition. There is considerable overlapping of the peaks of the samples with $`x_{\text{nom}}<0.02`$, but this should not be simply interpreted as a distribution of local $`x`$ spanning 0.015-0.018 for all these samples. In fact, the peaks in the derivative of $`E\left(T\right)`$ cannot be interpreted as the distribution function of the local $`x`$ after the $`T`$ scale is converted into an $`x`$ scale inverting $`T_t\left(x\right)`$, since the step of $`E\left(T\right)`$ has an intrinsic width due to the progressive phonon softening above $`T_t`$ and to the relaxational character of the domain wall motion below $`T_t`$. The undoped samples, for example, certainly do not have a spread of Sr concentration, but the widths of their lorentzians (not shown here) are $`1622`$ K, compared to 10 K of the sample with $`x_{\text{nom}}=0.016`$. In the following we will refer to the doping level obtained from the position of the step in the Young’s modulus.

Figure 2 presents the anelastic spectra of the five samples measured exciting the first flexural mode, around 1 kHz. The sharp rise of dissipation at the lowest temperatures is the tail of an intense peak attributed to the tunneling-driven tilt motion of a fraction of the O octahedra. The shift of the peak to lower temperature with increasing doping would be due to a direct coupling between the hole excitations and the tunneling-driven local tilts of the O octahedra: the more the hole excitations and the faster the local tilting. We are concerned with the steplike feature at $`T_g\left(x\right)0.2/x`$ for $`x>0.02`$, and below $`10`$ K for $`x<0.02`$. That such a dissipation rise occurs at $`T_g`$ and is steplike or at least consists of a plateau rather than a peak is particularly evident at higher doping, when the influence of the tail of the peak at lower temperature is less important; here only two curves with $`x0.02`$ are presented. The dissipation rise has been attributed to the stress-induced movement of the boundaries between the clusters of quasi-frozen antiferromagnetically correlated spins which form the CSG phase. The main result is that also the samples with $`x<0.02`$ present a similar feature around 10 K. With lowering $`x`$ below 0.02, the position of the dissipation step remains unchanged, while its intensity decreases; this is more difficult to see for the sample with $`x=0.0155`$, since the tail of the low temperature peak is shifted to higher temperature (consistent with the lowest doping in spite of the nominal $`x`$). Considering also previous data, it appears that, starting from high doping and lowering it, the dissipation step first rises in temperature according to $`T_g\left(x\right)`$ and increases in intensity down to $`x0.02`$; below that doping its temperature remains unchanged and the intensity decreases.

## III DISCUSSION

One would expect that the dissipation rise near 10 K for $`x<0.02`$ corresponds to the peak in the spin susceptibility occurring when the magnetic fluctuations slow down below the measuring frequency during the freezing process into the SG phase, as in the NQR measurements; also the temperature region would be as expected, since, according to the generally accepted phase diagram, at $`x=x_c`$ both $`T_f\left(x\right)`$ and $`T_g\left(x\right)`$ merge at $`1015`$ K. However, there are two features signalling the different nature of the absorption rise below $`10`$ K: one is the shape, certainly different from the sharp peak of the NQR relaxation rate, and the other is the total absence of the expected linear shift toward lower temperature when the doping is reduced. This is better demonstrated by the derivatives of the $`Q^1\left(T\right)`$ curves in Fig. 3, where the negative peaks at high temperature, labelled with the estimated doping values, correspond to the steps in $`Q^1\left(T\right)`$. The three peaks with $`x<0.02`$ even shift to higher temperature with decreasing $`x`$, although the effect is very small and most likely attributable to the increasing influence of the tail of the peak at lower temperature. Instead, according to the generally adopted phase diagram with $`T_f\left(x\right)x\times \left(815\text{K}\right)`$ the three peaks should display an overall shift by 2.3 K in the opposite direction; it is therefore clear from Fig. 3 that there is no relation between the acoustic dissipation step and $`T_fx`$.

The spectra of the samples with $`x<0.02`$ strictly resemble those for $`x>0.02`$, and the latter have been successfully interpreted in terms of onset below $`T_g`$ of the motion of the domain walls between different spin clusters in the CSG phase. The decrease of the amplitude $`\mathrm{\Delta }`$ of the absorption step above $`x=0.02`$ has been semiquantitatively explained considering a simple model of the CSG phase, with the Sr atoms acting as pinning points for the domain walls, which coincide with the hole stripes. The relaxation strength is expected to be of the form $`\mathrm{\Delta }nl^\alpha `$, where $`n`$ is the volume concentration of the domain walls, $`l`$ is the distance between pinning points, and $`\alpha `$ turns out $`2.5`$, intermediate between the case of the motion of dislocations and of walls between ferromagnetic domains. It was also possible to observe the pinning of the walls by the low-temperature tetragonal lattice modulation in samples doped with Ba instead of Sr. The mechanism of dissipative motion of the domain walls is well known for ferromagnetic materials, and is possible also for an AF state, if an anisotropic strain is coupled with the easy magnetization axis. In this case, the elastic energy of domains with different orientations of the staggered magnetization are differently affected by a shear stress, and the lower energy domains grow at the expenses of the higher energy ones providing the coupling between strain and domain wall motion. The same mechanism could not be applied to the SG state, even admitting the existence of hole stripes below 10 K, since they would move into a uniform long-range ordered AF matrix, instead of separating inequivalent domains, and other more subtle mechanisms should be invoked.

It has to be concluded that the anelastic spectra for $`x<0.02`$ (Fig. 2a) cannot be justified in terms of the spin-glass phase with $`T_f\left(x\right)`$ increasing linearly with $`x`$; rather, they can be explained in a straightforward way within the physical picture derived from the inelastic neutron scattering measurements. According to these observations, there is an electronic phase separation below $`x=0.02`$, where domains with fixed $`x_2=0.02`$ coexist with regions of undoped material with $`x_1=0`$. Such domains have sizes estimated in several hundreds of Å within the CuO2 planes, and display the same ”diagonal stripe” correlations which are observed in the CSG state for $`x0.02`$. The volume fraction of these domains has also been verified to increase linearly with $`x`$, as expected. In this phase separation picture, the elastic energy dissipation curves for $`x<0.02`$ simply contain the step-like increase due to the motion of the hole stripes within the $`x_20.02`$ domains, whose sizes are sufficiently large that they appear as the homogeneous phase $`x=0.02`$ with the same $`T_g10`$ K. The reduced volume fraction at lower doping simply results in a reduced amplitude of the elastic energy absorption.

The fact that the temperature at which our measurements reveal the presence of spin domains is much lower than the temperature of 30 K reported by Matsuda et al. also agrees with what is observed in the CSG state. In fact, the onset temperature for freezing toward the glass state, $`T_{on}`$, depends on the time scale at which the experimental technique probes the system, and decreases as the angular frequency $`\omega `$ of the probe decreases. In neutron scattering experiments $`\omega `$ is of the order of $`410^{11}`$ s-1, whilst for the anelastic spectroscopy $`\omega 10^4÷10^5`$ s-1, resulting in a factor 2.5 between the two $`T_{on}`$ for at $`x0.03`$; this is consistent with the $`T_{on}`$ $`=10÷15`$ K deduced from the curves of Fig. 2a.

The dissipation curves of Fig. 2a clearly indicate that the same dissipation mechanism of the CSG phase is operative also at $`x<0.02`$, and therefore are in good agreement with the phase separation picture from the neutron scattering experiment; however, these curves alone are not sufficient for clearly discriminating between a neat phase separation into $`x_10`$ and $`x_20.02`$ or a situation with a smooth transition from the CSG to the SG state over a wide concentration range around 0.02. For example, according to the theoretical model proposed by Gooding et al., at low temperature the holes localize near the Sr dopants and circulate over the four Cu atoms neighbors to Sr, inducing a distortion of the surrounding Cu spins, otherwise aligned according to the prevalent AF order parameter. The spin texture arising from the frustrated combination of the spin distortions from the various localized holes produces domains with differently oriented AF order parameters, which can be identified with the frozen AF spin clusters. In this model there is no clear boundary between the SG and CSG state, and the present data could also fit into such a description.

The possibility has also to be considered that the discrepancy between the present results below $`x0.02`$ and the canonical phase diagram is due to inhomogeneous doping of the samples at a microscopic level. Errors in evaluating the actual doping and inhomogeneous doping would produce particularly large shifts and uncertainty in the determination of $`T_f`$ and $`T_g`$, since these temperatures appear to merge around the value $`1015`$ K at $`x=x_c=0.02`$ , where they strongly depend on doping. Indeed, there is a wide scattering of experimental data in that region of the phase diagram, which can be due to uncertainties in $`x`$, but also to the fact that the transition between SG and CSG is not as sharp as supposed. The present data support the second explanation. In fact, we cannot exclude the possibility of inhomogeneous doping of our samples with certainty, only based on the width of the HTT/LTO phase transition, because we are not able to carry out a rigorous analysis of the curves of Fig. 1. Still, we note that if the high temperature of the step of samples with average $`x<0.02`$ is due to regions with $`x>0.02`$, then we should expect a similar spread of the local doping also for the samples with $`x>0.02`$, since all the samples have been prepared in the same way. Then, the sample with $`x=0.0273`$ should exhibit a partial onset of dissipation already above 10 K, due to regions with $`x0.02`$. The sample with $`x=0.0273`$, instead, displays a neat rise below 8 K, excluding the presence of regions with $`x0.02`$. It can be concluded that the present measurements demonstrate the presence of CSG domains also for $`x<0.02`$, and the insensitivity of the temperature of the acoustic anomaly on doping supports the proposal of an intrinsic phase separation at $`x<0.02`$.

## IV CONCLUSION

We measured the anelastic spectra of La2-xSrxCuO4 for both $`x<0.02`$ and $`x>0.02`$, spanning the region of the phase diagram where the transition from the SG to the CSG phase is expected. The same step-like feature in the elastic energy loss function that has been attributed to the domain wall motion in the CSG state is also found for $`0.015<x<0.02`$. The step cannot be related to the SG phase whose onset is usually assumed to be $`T_fx`$, since its temperature remains locked at the value found for $`x=0.02`$. The present data clearly show that CSG regions exist also for $`x<0.02`$ and therefore there is no neat separation between the SG and CSG regions. They can be naturally interpreted in the framework of the phase separation model for $`x<0.02`$, recently proposed by Matsuda et al. to explain their inelastic neutron scattering measurements.

## Acknowledgments

The authors thank Prof. A. Rigamonti for useful discussions. This work has been done in the framework of the Advanced Research Project SPIS of INFM.


text:
# Exact enumeration of the Critical States in the Oslo Model

## Abstract

We determine analytically the number $`N_{}(L)`$ of recurrent states in the $`1d`$ Oslo model as a function of system size $`L`$. The solution $`N_{}(L)=\frac{1+\sqrt{5}}{2\sqrt{5}}\left(\frac{3+\sqrt{5}}{2}\right)^L+\frac{\sqrt{5}1}{2\sqrt{5}}\left(\frac{3\sqrt{5}}{2}\right)^L`$ is in exact agreement with the number enumerated in computer simulations for $`L=110`$. For $`L1`$, the number of allowed metastable states in the attractor increases exponentially as $`N_{}(L)c_+\lambda _+^L`$, where $`\lambda _+=\frac{3+\sqrt{5}}{2}`$ is the golden mean. The system is non-ergodic in the sense that the states in the attractor are not equally probable.

Some non-equilibrium systems can be described by simple physical laws. Despite this being the case, we find complex and non-linear behavior in many such systems. A slowly driven granular system illustrates this phenomenon spectacularly. If grains are deposited onto a finite surface, they will eventually build into a pile. The pile does not collapse and flatten due to the friction between the grains. We call the different profiles of piles produced by this process, metastable states. The addition of grains and the dissipation events (avalanches) form a delicate feedback loop. Avalanches initiated by dropping grains onto the pile will ensure that the slope does not become too steep, while the addition of grains prevents the slope from becoming too shallow. Eventually, the system settles into an attractor where, the average avalanche size of the system diverges in response to small perturbations induced by adding grains. The distribution of avalanche sizes becomes scale invariant and a small perturbation may initiate small as well as large avalanches. As such, the system is intrinsically non-linear and displays characteristics similar to equilibrium systems poised at a critical point. This has been observed experimentally, for real granular systems where friction dominates over inertial effects Nature , as well as numerically BTW87 ; Oslo1 . It is an example of systems displaying self-organized criticality BTW87 ; PBak ; HJensen .

We study the attractor of a simple $`1d`$ granular model known as the Oslo model. It successfully models slowly driven granular systems displaying self-organized criticality Nature ; Oslo1 . The model has been shown to be a member of a large universality class of $`1d`$ systems, which includes the de-pinning transition of an interface dragged through a random medium and the deterministic Burridge-Knopoff train model for earthquakes Maya ; BK ; Train . Thus, the Oslo model plays the role of the “Ising model” for $`1d`$ self-organized critical systems and has been widely studied Maya ; Train ; Priezzhev ; ref1 ; ref2 ; ref3 ; ref4 ; ref5 ; ref6 ; ref7 ; ref8 ; ref9 ; ref10 ; ref11 . However, in contrast to various deterministic models of granular systems where exact results do exist Dhar1 ; Dhar2 , the Oslo model has resisted analytic treatment Priezzhev . Only a handful of exact results are known for non-equilibrium systems in general Derrida  and self-organized critical systems in particular Dhar1 ; Dhar2 ; Tang . Any analytic result for the Oslo model would be very valuable.

We are able to determine the number $`N_{}(L)`$ of states in the attractor of the Oslo Model. We show that $`N_{}(L)=\frac{1+\sqrt{5}}{2\sqrt{5}}\left(\frac{3+\sqrt{5}}{2}\right)^L+\frac{\sqrt{5}1}{2\sqrt{5}}\left(\frac{3\sqrt{5}}{2}\right)^L`$. For system sizes $`L1`$, the number increases exponentially fast as $`N_{}(L)c_+\lambda _+^L`$, where $`\lambda _+`$ is the golden mean. This is a clear signature of the complexity displayed by this simple model. Moreover, the system is non-ergodic. The recurrent states are not visited with the same frequency.

The Oslo model is defined on a one-dimensional discrete lattice consisting of $`L`$ sites, $`i=1,2,\mathrm{},L`$. There is a vertical wall at the left boundary next to site $`i=1`$. The system is open at the right boundary next to site $`i=L`$, where grains are allowed to leave the system. We will denote the respective heights of the pile starting from the vertical wall by an $`n`$-tuple $`𝐡=(h_1,h_2,\mathrm{},h_i,\mathrm{},h_L)`$, see Fig. 1.

In accordance with the experiment in Ref. Nature , grains are only inserted at $`i=1`$, the site next to the vertical wall. In the following discussion, it will be convenient to refer to local slopes instead of heights. The slopes are denoted by $`𝐳=(z_1,z_2,\mathrm{},z_i,\mathrm{},z_L)`$ with the definition $`z_i=h_ih_{i+1}`$ where $`h_{L+1}=0`$. Each site $`i`$ is assigned a critical slope, $`𝐳^c(t)=(z_1^c(t),z_2^c(t),\mathrm{},z_i^c(t),\mathrm{},z_L^c(t))`$. If the local slope exceeds the critical slope, $`z_i>z_i^c(t)`$, a grain at site $`i`$ will topple to the site $`i+1`$, that is $`h_ih_i1,h_{i+1}h_{i+1}+1`$. As a consequence, the slopes at neighboring sites $`i\pm 1`$ may also exceed their critical slopes and topple, consequently causing an avalanche that propagates until the system settles in a metastable state with $`z_iz_i^c`$ for all $`i`$.

The critical slope $`z_i^c(t)`$ is a function of time as it is chosen randomly every time site $`i`$ topples. We denote the minimum allowed value of the random critical slopes $`z^{c(min)}`$ and similarly the maximum allowed value $`z^{c(max)}`$. In general, we can use any discrete probability distribution of critical slopes. However, we will consider the simplest case where $`z_i^c`$ is chosen at random between the $`1`$ or $`2`$, $`z_i^c(t)\{1,2\}`$ when site $`i`$ topples, such that $`z^{c(min)}=1`$ and $`z^{c(max)}=2`$.

The algorithm for the Oslo model is defined as follows.

1. Initialize the critical slopes $`z_i^c`$ and the system in a arbitrary metastable state with $`z_iz_i^c`$ for all $`i`$.
2. Add a grain at site $`i=1`$, i.e., $`z_1z_1+1`$.
3. If $`z_i>z_i^c`$, the site relaxes and

$`z_i`$ $``$ $`z_i2`$
$`z_{i\pm 1}`$ $``$ $`z_{i\pm 1}+1`$
except for the sites at the vertical wall $`i=1`$ and the open boundary $`i=L`$ where

$`z_1`$ $``$ $`z_12z_Lz_L1`$
$`z_2`$ $``$ $`z_2+1z_{L1}z_{L1}+1.`$
The critical slope at site $`i`$ is chosen randomly if it topples, that is $`z_i^c1\text{or}\mathrm{\hspace{0.33em}\hspace{0.33em}2}`$. A new metastable state is reached when $`z_iz_i^c`$ for all $`i`$.
4. Proceed to 2 and reiterate.

Note the separation of time scales which is built into the definition of the model. The addition of grains can only take place when the system has settled down into a metastable state. Thus the response of the system (the avalanches) is fast when compared with the interval between perturbations.

Metastable states $`_j`$ are rest states of the system. Different metastable states occur due to the different possible sets of critical slopes. Metastable states can also be divided into transient states $`𝒯_j`$ and attractor states $`_j`$. Transient states are not reachable once the systems has entered into the attractor, which are recurrent metastable states. The index $`j`$ denotes discrete time steps associated with the long time scale of the system.

When adding a grain to a metastable state $`_j`$, it evolves into a new metastable state $`_{j+1}`$ by the relaxation rules given above. Symbolically we write $`_j_{j+1}`$, where the arrow is a shorthand notation for the operation of adding sand and, if necessary, relaxing the system until it reaches a new metastable state. Starting from, say the empty lattice $`𝒯_1`$, which is a transient state as it will never be encountered again, we have

$$𝒯_1\mathrm{}𝒯_n_1\mathrm{}_{j1}_j_{j+1}\mathrm{}$$

After $`n`$ additions and associated relaxations (if any), the system reaches the attractor of the dynamics. For example, if we set $`z_i^c`$ independent of time, we arrive at the BTW model in one dimension BTW87 . Any grain introduced into the system after the attractor state has been reached will simply flow to the open boundary fall out. In this case, there is only one attractor state with $`z_i=z_i^ci`$. We want to enumerate the number of different states in the attractor $`\left\{_j\right\}`$ for the Oslo model, where $`z_i^c`$ is chosen randomly between $`1`$ or $`2`$.

We adopt the following nomenclature for sites

$`z_i<z^{c(min)}`$ $`\text{sinks,}z_i=z^{c(min)}\text{stable,}`$
$`z^{c(min)}<z_iz^{c(max)}`$ $`\text{critical,}z_i>z^{c(max)}\text{supercritical}.`$
Due to the way the problem is defined, $`z_i0`$ always. Also, supercritical sites are by definition not allowed in any metastable state. Fig. 2 shows a diagrammatic relation of these classifications.

From the relaxation rules, we find that states in the attractor are subject to the following two constraints.

Starting from the open boundary $`i=L`$, the first site ($`i<L`$) that is not a stable site ($`z_iz^{c(min)}`$) must be a critical site.

The first site that is not a stable state to the left of the sink must be a critical site. Equivalently, we cannot find a recurrent state in which two sinks are separated only by stable sites. The exception is at the closed boundary $`i=1`$, where the first site that is not stable must be the boundary itself.

We can show that these constraints must be fulfilled by assuming that each is not true and examining how that system would enter such as state. For example, Fig. 3(a) shows a case where the the first rule C1 does not hold. The possible preceding configurations are displayed in Fig. 3(b) and (c). Since these are not states in the attractor, the state displayed in Fig. 3(a) is not in the attractor. One could argue likewise for the second constraint, C2.

Note however, that by iterating these cases, we have only shown that these constraints are necessary but not sufficient. However, we make the conjecture that these constraints are also sufficient for enumerating all attractor states from exact agreement with simulations for system sizes $`L=1`$ to $`10`$, see Fig. 4.

For a system with $`L`$ lattice sites, we want to find the number of states $`N_{}(L)`$ that obey the constraints above. We can simplify the problem further by first considering the number states with a fixed height at site $`i=1`$, $`h_1`$. We introduce the parameter $`q=h_1L=_{i=1}^Lz_iL`$ as it simplifies the equations involved considerably. Observe that the range of $`q`$ is $`0qL`$. Let $`K(L,q)`$ describe the number of states in a system of $`L`$ sites and $`h_1=q+L`$. $`N_{}(L)`$ can then be determined by $`N_{}(L)=_{q=0}^LK(L,q)`$.

To generate all states $`\{𝐳\}`$ that obey the counting constraints for states in the attractor, consider the situation where we know all the allowed profiles for system sizes $`L^{}<L`$. We denote the family of $`n`$-tuples that represent allowed profiles for a system of size $`L`$ by $`(L)`$. $`N_{}(L)`$ is the number of $`n`$-tuples in $`(L)`$. We can generate $`(L)=\{𝐡\}`$ by extending $`(L1)`$ by one site at the closed boundary. This may be done using the following procedure - extend each $`𝐳^{}`$ to $`𝐳`$ by writing $`z_{i+1}=z_i^{}`$. Consider one of three possible values for the new site next to the closed boundary, $`z_1`$. $`z_1=2`$ for critical sites, $`z_1=1`$ for stable sites and $`z_1=0`$ for sinks. Lastly, check that the new states in each of the three cases obey the counting constraints (C1) and (C2).

Consider the three possible values for $`z_1`$, that give $`(L)`$ by extending $`(L1)`$ corresponding to three possible cases such that $`K(L,q)=_{i=1}^3K^{(i)}(L,q)`$.

1. When $`z_1=2`$, there is a critical site next to the closed boundary. The number of allowed profiles for a system of size $`L`$ and parameter $`q`$ must be the same as the number of allowed profiles with a system size of $`L^{}=L1`$. We relate $`q^{}`$ to $`q`$ by observing that $`q=h_1L=h_1^{}(L1)+1=q^{}+1`$. Therefore we assert that when $`z_1=2`$, $`K^{(1)}(L,q)=K(L1,q1)`$.
2. Similarly, when $`z_1=1`$ corresponding to a stable site next to the closed boundary all profiles for a system size $`L1`$ with $`q^{}=q`$ are allowed. In this case, $`K^{(2)}(L,q)=K(L1,q)`$.
3. Finally, the situation when $`z_1=0`$ means that a sink resides next to the closed boundary. The counting constraints state that the sink next to the closed boundary must be followed by a critical site with any number of stable sites in between. In this case, we consider $`𝐳^{}=(0,2,\mathrm{})`$ extending $`(L2)`$, $`𝐳^{}=(0,1,2,\mathrm{})`$ extending $`(L3)`$ and so on to give $`L`$. This involves summing all allowed metastable states $`(L^{})`$ with $`L^{}L2`$. As above, we find that $`q^{}=q`$ for all $`L^{}`$. Hence

$$K^{(3)}(L,q)=\underset{i=1}{\overset{L2}{}}K(i,q).$$

Note that this equation is a simple summation because of the appropriate choice of $`q`$.

Summing the three cases,

$`K(L,q)`$ $`=`$ $`K(L1,q1)+K(L1,q)+{\displaystyle \underset{i=1}{\overset{L2}{}}}K(i,q)`$
$`=`$ $`{\displaystyle \underset{i=1}{\overset{L1}{}}}K(i,q)+K(L1,q1)`$
$`=`$ $`2K(L1,q)+K(L1,q1)K(L2,q1).`$
Hence we calculate the number of states as a function of system size by summing through $`q`$

$`N_{}(L)`$ $`=`$ $`{\displaystyle \underset{q=0}{\overset{L}{}}}K(L,q)`$
$`=`$ $`3{\displaystyle \underset{q=0}{\overset{L1}{}}}K(L1,q){\displaystyle \underset{q=0}{\overset{L2}{}}}K(L2,q)`$
$`=`$ $`3N_{}(L1)N_{}(L2).`$
We have used the fact that $`0qL`$. $`N_{}(L)=3N_{}(L1)N_{}(L2)`$ is a Fibonacci like relation. Substituting a trial solution $`N_{}(L)=c_+\lambda _+^L+c_{}\lambda _{}^L`$, and demanding the initial conditions $`N_{}(1)=2`$ and $`N_{}(2)=5`$, we obtain the solution, see Fig 4.

$$N_{}(L)=\frac{1+\sqrt{5}}{2\sqrt{5}}\left(\frac{3+\sqrt{5}}{2}\right)^L+\frac{\sqrt{5}1}{2\sqrt{5}}\left(\frac{3\sqrt{5}}{2}\right)^L.$$
(1)
The golden mean $`\lambda _+=\frac{3+\sqrt{5}}{2}2.618\mathrm{}`$. Computer simulations were used to test these results for systems of size $`1`$ to $`10`$. They were found to be in full agreement.

This result should be compared with the size of the attractor for the Abelian Bak-Tang-Wiesenfeld sandpile determined by Dhar Dhar1 . He showed that $`N_{}=\text{Det}\mathrm{\Delta }`$, where $`\mathrm{\Delta }`$ is the toppling matrix. In $`2D`$ and for large system sizes, $`N_{}`$ increases exponentially fast with the number of sites in the system $`(3.21\mathrm{})^{L\times L}`$ Creutz . Moreover, it was proven that the BTW model is ergodic Dhar1 . We find, in contrast, that the Oslo model is non-ergodic in the sense that the states in the attractor are not visited with the same frequency. This is most easily demonstrated by showing that the probability of visiting stationary states in very small systems (e.g., $`L=2`$ which has a total of $`5`$ states in the attractor) is non-uniform. The number of states in the attractor $`N_{}`$ does not depend on the probability distribution with which $`z_c^i`$ is chosen to be $`1`$ or $`2`$, except for the trivial cases of $`P(z_c^i=1)=0`$ or $`1`$, but the distribution will affect the accessibility of a given microstate. Finally, note that the Oslo model can be mapped into an interface dragged through a random medium Maya . Thus, the results obtained for the Oslo model are directly applicable to the latter model as well.

Possible extensions. A pile can be divided into an inactive zone consisting of grains that never will take part in the dynamics and the active zone. The width $`\lambda _L`$ of the active zone plays an important role as it is related to various physical phenomena. E.g. the average transit times of individual grains has been measured to increase with system size as $`T_LL\lambda _LL^{1+\chi }`$. For the Oslo model $`\chi 0.25`$ while experimentally $`\chi 0.5`$ Oslo1 ; Nature . We believe that using a similar approach one might be able to calculate analytically the distribution of the height of the pile at the vertical wall $`P(h_1)`$ and the critical exponent $`\chi `$. This would be of great importance as this exponent will be the same in all the models belonging to the universality class of the Oslo Model.

We are grateful to S.C. Khoo, R. Jack and Á. Corral for their helpful comments and insights. K. C. gratefully acknowledges the financial support of U.K. EPSRC through Grants No. GR/R44683/01 and GR/L95267/01.

To whom correspondence should be addressed, E-mail: k.christensen@ic.ac.uk.