text:
# Percolation fractal exponents without fractal geometry

## Abstract

Classically, percolation critical exponents are linked to the power laws that characterize percolation cluster fractal properties. It is found here that the gradient percolation power laws are conserved even for extreme gradient values for which the frontier of the infinite cluster is no more fractal. In particular the exponent 7/4 which was recently demonstrated to be the exact value for the dimension of the so-called ”hull” or external perimeter of the incipient percolation cluster keeps its value in describing the width and length of gradient percolation frontiers whatever the gradient value. Its origin is then not to be found in the thermodynamic limit. The comparison between numerical results and the exact results that can be obtained analytically for extreme values of the gradient suggests that there exist a unique power law from size 1 to infinity which describes the gradient percolation frontier, this law becoming a scaling law in the large system limit.

Spreading of objects in space with a gradient of probability is most common. From chemical composition gradients to the distribution of plants which depends of their solar exposure, probability gradients exist in many inhomogeneous systems. In fact inhomogeneity is a rule in nature whereas most of the systems that physicists are studying are homogeneous as they are thought to be more simple to understand. In particular phase transitions or critical phenomena are studied in that framework, the simplest being percolation transition.

In this work, the opposite situation, a strongly inhomogeneous system is studied. We report the discovery that the validity of some percolation exponents can be extended to situations which are very far from the large homogeneous system limit. This is found in the frame of gradient percolation, a situation first introduced in the study of diffusion front. Surprisingly these exponents, up to now believed to pertain to large systems, are also verified in a limit that could be called the small system limit (S.S.L.) where some of the gradient percolation properties can be computed analytically. This letter presents a report concerning the 2D square lattice. The same results and discussion have been obtained for the triangular lattice and will be published elsewhere, together with the details of the analytical calculations for both lattices.

The gradient percolation (G.P.) situation is depicted in FIG. 1. The figure gives an example of a random distribution of points on a lattice with a linear gradient of concentration in the vertical direction. It is a 2D square lattice of size $`L_g\times L`$, where each point $`(x,y)`$ is occupied with probability $`p(x)=1x/L_g`$ ($`x`$ being the vertical direction in the figure). In gradient percolation there is always an infinite cluster of occupied sites as there is a region where $`p`$ is larger than the standard percolation (S.P.) threshold $`p_c`$. There is also an infinite cluster of empty sites as there is a region where $`p`$ is smaller than $`p_c`$. The object of interest is the G.P. front, the external limit (or frontier) of the infinite occupied cluster. It is constituted by the sites which belong to the occupied cluster and are first-nearest neighbours with empty sites belonging to the infinite empty cluster. It is shown in grey in FIG. 1. This front is a random object with an average position $`x_f`$, a statistical width $`\sigma _f`$ and a total length $`N_f`$. In so far that the G.P. front and the S.P. external perimeter (often called hull) have the same geometry their fractal dimension was first conjectured to be exactly equal to 7/4 in . This result was then demonstrated heuristically by Saleur and Duplantier and very recently it was proved mathematically by Smirnov.

The early G.P. studies were focussed at finding its relation with standard percolation. Let first recall the definitions. For $`0xL_g`$, $`n_f(x)`$ is the mean number of points of the front lying on the line $`x`$ per unit horizontal length. It measures the front density at distance $`x`$. The length $`N_f`$, the position $`x_f`$ and the width $`\sigma _f`$ of the front are then defined in terms of the $`n_f(x)`$ by

$$N_f=L\underset{x=0}{\overset{L_g}{}}n_f(x),x_f=\frac{_{x=0}^{L_g}xn_f(x)}{_{x=0}^{L_g}n_f(x)},$$

$$\mathrm{and}\sigma _f^2=\frac{_{x=0}^{L_g}(xx_f)^2n_f(x)}{_{x=0}^{L_g}n_f(x)}.$$

It was found that the mean front was situated at a distance where the density of occupation was very close to $`p_c`$ or $`p(x_f)p_c`$. This was verified numerically with such precision that the gradient percolation method is now often used to compute percolation thresholds. It was also found that:

1. The width $`\sigma _f`$ depends on $`L_g`$ through a power law $`\sigma _f(L_g)^{\nu /(1+\nu )}`$ where $`\nu =4/3`$ is the correlation length exponent  in dimension $`d=2`$ so that $`\sigma _f=(L_g)^{4/7}`$. The width $`\sigma _f`$ was also shown to be a percolation correlation length.
2. Secondly it was found that the front was fractal with a dimension $`D_f`$, numerically determined, close to $`1.75`$. The front length followed a power law $`N_f(L_g)^{\alpha _N}`$ with $`\alpha _N=(D_f1)\nu /(1+\nu )`$.
3. But also, it was numerically observed that the sum of these two exponents was very close to 1. If true this meaned that $`\nu /(1+\nu )+(D_f1)\nu /(1+\nu )=1`$ or $`D_f=1+1/\nu `$ = 7/4. This is how it was conjectured in  that $`D_f=7/4`$.

In that sense the G.P. power laws were thought to be linked to the S.P. exponent $`\nu `$ and to the fractality of the percolation cluster hull. Up to now, these facts were considered to be strictly valid only in the large system limit.

However, if true, and we know now that 7/4 is the exact value, there follows an intriguing relation, namely $`\sigma _f^{D_f}`$ is exactly proportional to $`L_g^1`$. This means that the number of surface particles within the correlation length is exactly proportional to $`L_g`$. This is particularly striking for diffusion fronts. Diffusion of particles from a source results in a concentration gradient and an associated G.P. situation. In that frame, the above result means that, if $`L_g`$ particles have diffused on a vertical row, there is on average the same (or a constant fraction of) number of particles on the correlated surface (surface content of a box with a lateral size equal to the statistical width). This fact seems a priori to have nothing to do with scaling, percolation and the thermodynamic limit. From this point of view it is possibly the consequence of a conservation law and if such a conservation exists, it should apply also for extreme gradients corresponding to $`L_g`$ of a few units. In particular it should apply to the very extreme $`L_g=1`$, $`2`$ and $`3`$ for which exact values of $`x_f`$, $`N_f`$ and $`\sigma _f`$ can be calculated analytically.

For $`L_g=2`$ or $`L_g=3`$, given a site on a line $`x`$, one can describe all the configurations such that the point belongs to the front, and compute their probability. For example, for $`L_g=2`$ all the occupied sites on the line $`x=1`$ belong to the front, and a site on the line $`x=0`$ will belong to the front if at least one of its three neighbours on the line $`x=1`$ is empty. Thus we get in this case, $`n_f(1)=1/2`$ and $`n_f(0)=1(1/2)^3=7/8`$. With the same kind of arguments, we can make the computations for $`L_g=3`$ but the geometry of connected sets is more complex as there are more configurations to consider. This however can be done exactly. One obtains: for $`L_g=1`$, $`\sigma _f=0`$;
$``$ for $`L_g=2`$:

$$N_f/L=\frac{11}{8},x_f=\frac{4}{11}\mathrm{and}\sigma _f=\frac{2\sqrt{7}}{11}.$$

$``$ for $`L_g=3`$:

$$N_f/L=\frac{9401}{5832},x_f=\frac{6966}{9401}\mathrm{and}\sigma _f=\frac{9\sqrt{576049}}{9401}.$$

The problem is then to compare the numerical G.P. laws to these exact values. As will be shown, the numerical results verify the above power laws with such precision that the question arises of the existence of a simple mathematical power law extending from $`L_g=1`$ to infinity. To try to answer this question we proceed in 2 steps. First we test these laws on the numerical results obtained for $`L_g`$ between 4 and 50 for the square and triangular lattices by searching the best numerical power laws followed by the width. Considering arbitrary exponents $`\alpha `$ between $`1.6`$ and $`1.9`$, we study $`\sigma _f^\alpha `$ as a function of $`L_g`$ between $`4`$ to $`50`$. For each $`\alpha `$ value, there is a best line $`\sigma _f^\alpha =a_\alpha (L_g+b_\alpha )`$ fitting the numerical $`\sigma _f^\alpha `$. The introduction of the term $`b_\alpha `$ is justified by the fact that when a power law is verified for large systems it includes always the possibility that a small (as compared to the system size) term could contribute but in a negligible manner. But here the size itself is small or very small. On the other hand, one should remark that for $`L_g=1`$ the width is strictly 0 so that some negative value of $`b_\alpha `$ should be present. In the next step, the mean error $`d(\alpha )`$, defined by $`d(\alpha )^2=(1/47)_{L_g=4}^{50}(\sigma _f(L_g)^\alpha a_\alpha (L_g+b_\alpha ))^2`$, is measured numerically as a function of $`\alpha `$. The results are shown in FIG. 2. There is a clear minimum for $`\alpha =1.75`$, showing that this exponent gives the best power law fit. Once the best fit with the empirical data is made one has the best values for the parameters $`a`$ and $`b`$: $`a=0.297`$ and $`b=1.09`$. Note that $`b`$ should be strictly equal to $`1`$ in order to obtain a null width for the trivial case $`L_g=1`$.

An other verification of the extreme G.P. power laws can be obtained from the study of the front length or of the quantity $`(N_f/L)^{7/3}`$ as a function of $`L_g`$. In FIG. 3, the diamonds represent the values of $`(N_f/L)^{7/3}`$ and the best linear fit has equation $`Y=c(L_g+d)`$ with $`c=0.845`$ and $`d=0.88`$. This shows indeed that the exponents $`4/7`$ and $`3/7`$ can be used down to the steepest gradients for which the frontier is no more fractal.

One can then extrapolate the $`\sigma _f`$ values to the case $`L_g=1`$, $`2`$ and $`3`$. The results are given in Table I. One observes that the numerical extrapolations correspond to the exact values with an apparent good precision. However no firm conclusion can be drawn without discussion of the numerical uncertainties. The values shown on FIG. 3 are averaged over $`50`$ trials on a length $`L=10^5`$. As the fit occurs through a power law it is difficult to give a confidence interval for the coefficients $`a`$ and $`b`$ (obtained from a least square linear regression on the values of $`\sigma _f^{7/4}`$).

In order to obtain a better control on the numerical precision of $`a`$ and $`b`$ we made extensive computations of the two cases $`L_g=4`$ and $`L_g=5`$ with 100 trials on a length $`L=5.10^5`$. Doing so we obtain the mean values with their standard deviation: $`\sigma _f(4)^{7/4}=0.8658\pm 0.0009`$ and for $`\sigma _f(5)^{7/4}=1.1610\pm 0.0013`$. Thus if we compute the equation of the line $`\overline{a}(L_g+\overline{b})`$ which interpolates the two points $`(4,\sigma _f(4)^{7/4})`$ and $`(5,\sigma _f(5)^{7/4})`$, we obtain $`\overline{a}=0.2952\pm 0.0022`$ and $`\overline{b}=1.066\pm 0.041`$. This last result shows that the value $`1`$ is not incompatible with $`\overline{b}`$ and its statistical error. Given the numerical values for $`L_g=4`$ and $`5`$, we can also get extrapolated values for $`\sigma _f`$ for $`L_g`$ smaller, together with their confidence interval. The result (see Table I) is that the predicted values are very close to the exact ones. For $`(N_f/L)^{7/3}`$, in the same way we obtain a linear interpolate of the values for $`L_g=4`$ and $`5`$, with coefficients $`\overline{c}=0,893\pm 0.014`$ and $`\overline{d}=0.427\pm 0.074`$.

At this point one could conclude that the numerical results are compatible with the existence of a single mathematical law for the width dependence, this law working from $`L_g=1`$ to infinity. Note however that the quality of the random number generator should intervene. As long as a mathematical proof has not been given, it is thus impossible for us to conclude on the exact values of the coefficients $`a`$ and $`b`$.

The above discussion bears on the validity of the numerical results. But the question of a unique mathematical power law can also be studied from the point of view of the exact results only. One can note that, remarkably, the values of the term $`b`$ in the fit of the width are close to $`1`$. If the above power law exists from $`L_g=1`$ to infinity it suggests that the real value of $`b`$ is exactly -1 as the width is null in the trivial case $`L_g=1`$. As we have exact values, one can compute the equation of the line $`y=a(L_g+b)`$ defined by the two points $`(2,\sigma _f(2)^{7/4})`$ and $`(3,\sigma _f(3)^{7/4})`$. One obtains $`a=0.294`$ and $`b=1.055`$. These values are close to the values obtained from the above numerical fit ($`a=0.297`$ and $`b=1.09`$), and here again $`b`$ is close but not equal to $`1`$.

Why is there a small mismatch with the simplest law? The answer to this question is two-fold. First, it is possible that the exact power law is not valid for $`L_g=1`$ or both $`L_g=1`$ and $`2`$. These cases could be ”anormal” as for these values there is no Grossman-Aharony effect . It is then possible that these two cases do not enter the unique mathematical law. Secondly the small discrepancy could be related to the fact that the frontier definition considers only the occupied sites. It gives to these sites a privilege role whereas one should also consider the frontier of the empty cluster. In fact this is not new in percolation studies  where it was shown that the barycenter between the frontier of the occupied cluster and the frontier of the empty cluster was a more natural object. It notably permitted better computations of the percolation threshold. We have studied the statistical width of the local barycenter which can also be computed exactly for $`L_g=2`$ or $`3`$. The results show the same behavior as described above i.e. a $`b`$ value close but not equal to $`1`$. The question remains then open to define the nature of the geometrical object which really would show a corresponding $`b`$ value equal to $`1`$. The answer has certainly interesting consequences for the S.P. problem itself.

The above results have a simple practical consequence. Given an irregular non-fractal interface, for instance the grey line shown in Fig. 1, one can determine if it belongs to gradient percolation by measuring its statistical width $`\sigma `$ and the average value of $`N_f/L`$. Through the results described above one can find $`L_g`$ from the measured $`\sigma _f`$. Then one can check if $`N_f/L`$ satisfies approximately the G.P. power law as a function of $`L_g`$. This gives an intrinsic way to check if a given interface is of the gradient percolation type without previous knowledge of the gradient. This is important as there exist cases where irregular fronts like corrosion fronts belong to G. P. even though the gradient which built the interface is no more present and only the interface remains .

In summary, it has been shown that the classical power laws of gradient percolation can be extended to extreme gradients with the same fractal exponents although the systems present no fractal geometry. Four comments can be drawn on these results. First, extreme gradient situations can be found in diffused contacts between materials. The fact that the contact geometry can be described by the same set of exponent whatever the width of the diffused layer is certainly of help to understand the properties of these contacts. There is an intrinsic method to find if a given rough interface belongs to gradient percolation without knowledge of the gradient. Secondly our results imply that there exists a conservation law which stipulates that the length of the correlated frontier is strictly proportional to the gradient length. Thirdly, the fact that the same exponent has been found for the square and the triangular lattice even for extreme gradients, in other words for small systems, suggest that universality here is not related to the neglect of the microscopic details of the interactions. Here there is no coarse scale and still universality is verified. Finally, the fact that the exponents 4/7 and 3/7 are valid down the smallest $`L_g`$ values (or the steepest gradients) suggests that these exponents play the same type of role here that the exponent 1/2 intervening in the fluctuations of the sum of independent identical random variables. In that last case the exponent apply to any number of random variables starting from 1, 2 or 3 up to infinity. The exponents 4/7 and 3/7 may then play a stronger role that critical exponents which only exist only in the thermodynamic limit.

The Centre de Mathématiques et de leurs Applications and the Laboratoire de Physique de la Matière Condensée are “Unité Mixte de Recherches du Centre National de la Recherche Scientifique” no. 8536 and 7643.


format: text

text:
# A comparison of extremal optimization with flat-histogram dynamics for finding spin-glass ground states

## 1 Introduction

Optimization with methods motivated from real physical processes is an active field of research. Simulated annealing  and genetic algorithm  are two well-known examples. In particular, there have been a large variety of methods proposed to find spin-glass ground states . Recently, Boettcher and Percus  introduced ‘extremal optimization’ (EO) inspired by models of self-organized criticality , which gave impressive performance.

Most of the heuristic optimization methods (including simulated annealing and genetic algorithm) are designed to find ground states only, thus it is not possible to give correct thermodynamics from a simulation. On the other hand, multi-canonical ensemble simulation , $`1/k`$-sampling , parallel tempering , and recent flat-histogram dynamics  are constructed for equilibrium thermodynamics, but can also be used as methods for optimization. A study of optimization by flat-histogram algorithm on the two-dimensional spin-glass model is carried out in ref. . Unlike simulated annealing and other heuristic methods, we note that these methods do not have any adjustable parameters. It is useful to know the efficiencies of this second class of methods when used as an optimization tool.

In this paper, we make a comparative study of the extremal optimization and flat-histogram/equal-hit dynamics. We compare four algorithms: EO at $`\tau =1`$ with a continuous approximation in the probability of choosing a site, original EO at optimal value of $`\tau =1.15`$, single-spin-flip flat-histogram dynamics, and equal-hit algorithm with $`N`$-fold way. It is found that EO at the value $`\tau =1.15`$ is very good for both two- and three-dimensional Ising spin glasses. The equal-hit algorithm with $`N`$-fold way is also competitive. For large systems, equal-hit appears even slightly better than EO. It is useful to have the efficiency of EO but still give equilibrium results. To this end, we introduce a rejection step in EO, thus turning EO into an equilibrium simulation method.

## 2 Single-spin-flip algorithms

In the following, we specialize our discussion in the context of spin models, and particularly the spin-glass model . The energy function is defined by

$$E(\sigma )=\underset{i,j}{}J_{ij}\sigma _i\sigma _j,$$
(1)
where the spin $`\sigma _i`$ takes on value $`+1`$ or $`1`$ with $`i`$ varying over a hypercubic lattice in $`d`$ dimensions. The coupling constant $`J_{ij}`$ for each nearest neighbor pair $`i,j`$ takes on a random value of $`+J`$ and $`J`$ with equal probability. We impose a constraint $`_{i,j}J_{ij}=0`$. The spin glass is known to be one of the hardest problems  to find the state $`\sigma `$ that minimizes $`E(\sigma )`$.

A single-spin-flip with rejection is described by a Markov chain transition matrix of the form

$$W(\sigma \sigma ^{})=\delta _N(\sigma ,\sigma ^{})\frac{1}{N}a(\sigma \sigma ^{}),\sigma \sigma ^{},$$
(2)
where $`\delta _N(\sigma ,\sigma ^{})=1`$ if $`\sigma ^{}`$ is obtained from $`\sigma `$ by a single spin flip, and 0 otherwise. The factor $`1/N`$ represents the random selection of a spin, where $`N`$ ($`=L^d`$) is the number of spins. $`a(\sigma \sigma ^{})`$ is the flip rate. If we choose $`a(\sigma \sigma ^{})`$ according to Metropolis rate,

$$a(\sigma \sigma ^{})=\mathrm{min}(1,\frac{f\left(E(\sigma ^{})\right)}{f\left(E(\sigma )\right)}),$$
(3)
we can realize equilibrium distribution with the probability of states distributed according to $`f\left(E(\sigma )\right)`$. Some choices are:

$$f(E)=\{\begin{array}{cc}\mathrm{exp}\left(E/(k_BT)\right),\hfill & \text{canonical ensemble;}\hfill \\ 1/n(E),\hfill & \text{multicanonical ensemble;}\hfill \\ 1/_{\mathrm{}}^En(E^{})𝑑E^{},\hfill & 1/k\text{ sampling,}\hfill \end{array}$$
(4)
where $`T`$ is temperature, $`k_B`$ is Boltzmann constant, and $`n(E)`$ is density of states at energy $`E`$.

Arbitrary choice of the flip rate $`a(\sigma \sigma ^{})`$ in general would not give one important property of the equilibrium systems, i.e., the microcanonical property that the probability distribution of the configuration $`\sigma `$ is a function of energy $`E`$ only. For example, the original broad histogram rate

$$a(\sigma \sigma ^{})=\mathrm{min}(1,\frac{N_{Zk}(\sigma )}{N_k(\sigma )}),$$
(5)
and the random walk rate of Berg  do not have microcanonical property, where $`N_k(\sigma )`$ is the number of possible moves of class $`k`$ in the state $`\sigma `$; we associate a class for each site $`i`$ with a number from 0 to $`Z`$ ($`=2d`$) by a scaled energy change $`k=\left((E(\sigma ^{})E(\sigma ))/J+2Z\right)/4`$, i.e.,

$$k=\frac{1}{2}\underset{j}{}(J_{ij}\sigma _i\sigma _j+1).$$
(6)
$`N_k(\sigma )`$ is the number of such sites having a class $`k`$.

The single-spin-flip version with rejection can be turned into a rejection-free $`N`$-fold way  simulation where a class is chosen with probability

$$P_k=\frac{a_kN_k(\sigma )}{A(\sigma )},A(\sigma )=\underset{k=0}{\overset{Z}{}}a_kN_k(\sigma ),$$
(7)
where $`a_k`$ is $`a(\sigma \sigma ^{})`$ associated with an energy change indicated by class $`k`$. A spin in that class is picked up at random, and the flip is always accepted. Thermodynamic quantities need to be weighted by a factor $`1/A(\sigma )`$.

EO  is somewhat related to $`N`$-fold way in the sense that a class is chosen with some probability $`P_k`$, and a spin in that class is picked up and always flipped. The EO algorithm can be stated as follows: we classify the site by its ‘fitness’ $`k`$. There are $`Z+1`$ possible values for $`k`$. In the general EO, the sites are sorted according to the fitness $`k`$. Since there are only a small number of possible values in the $`\pm J`$ spin-glass model, the sorting is not necessary. We simply make a list of sites in each category. We pick a class according to the probability $`P_k`$, and then choose a spin in that class and flip with probability one. The corresponding transition matrix is then

$$W(\sigma \sigma ^{})=\delta _N(\sigma ,\sigma ^{})P_k\frac{1}{N_k},\sigma \sigma ^{}.$$
(8)
The original choice of EO is to take

$$P_k\underset{N_0+N_1+\mathrm{}+N_{k1}<iN_0+N_1+\mathrm{}+N_k}{}i^\tau ,$$
(9)
with $`\tau `$ being a parameter of the algorithm. We define the standard EO to be a continuous approximation to the above sum at $`\tau =1`$ with the analytical expression by

$$P_k=\frac{1}{\mathrm{ln}(N+1)}\mathrm{ln}\frac{1+\underset{j=0}{\overset{k}{}}N_j}{1+_{j=0}^{k1}N_j}.$$
(10)
The number 1 in the numerator and denominator are introduced somewhat arbitrarily to avoid divergence when the sum of $`N_k`$ is zero. An optimized EO will be the discrete version, Eq. (9), with $`\tau `$ that gives best performance; we use $`\tau =1.15`$ as recommended in ref . To realize the power-law distribution, we generate an integer $`i=\xi ^{1/(1\tau )}`$, where $`\xi `$ is a uniformly distributed random number between 0 and 1, and pick a corresponding site $`i`$ ordered by the class. We have used $`\mathrm{}`$ for the floor function.

## 3 Comparison of EO with flat-histogram and equal-hit algorithms

The flat-histogram algorithm  is a special choice of the flip rate

$$a^{\mathrm{FH}}(\sigma \sigma ^{})=\mathrm{min}(1,\frac{N_{Zk}(\sigma ^{})_E^{}}{N_k(\sigma )_E}),$$
(11)
where the angular brackets with subscript $`E`$ denote a microcanonical average of the quantity $`N_k(\sigma )`$ at energy $`E`$; the starting state $`\sigma `$ has energy $`E`$ and the final state $`\sigma ^{}`$ has energy $`E^{}`$. This particular choice of the rate gives a flat distribution for the energy histogram, $`H(E)=n(E)f(E)=\mathrm{const}`$, $`n(E)`$ is density of states. This is one way to realize multicanonical ensemble.

In the $`N`$-fold way equal-hit algorithm , we perform the usual $`N`$-fold-way move (thus rejection-free) which is constructed from the following single-spin-flip rate:

$$a^{\mathrm{EQ}}(\sigma \sigma ^{})=\mathrm{min}(1,\frac{A_EN_{Zk}(\sigma ^{})_E^{}}{A_E^{}N_k(\sigma )_E}).$$
(12)
We note that $`A_E^1=1/A_N`$, where $`\mathrm{}_N`$ is average over the $`N`$-fold way samples. In equal-hit algorithm, it is guaranteed by construction that we change state in every move, and the distribution of the visits to different energies is flat.

Since the microcanonical averages used in the flip rates are not known before the simulation, we use running average to replace the exact microcanonical average. It appears that this is a valid approximation and should converge to the correct values for sufficiently long runs. For a truly exact algorithm (in the sense of realizing microcanonical property), it is sufficient with a two-pass simulation. The first pass uses a running average; in the second pass, we use a multicanonical rate determined from the first pass.

Many different criteria are used to measure the effectiveness of an optimization algorithm, such as the fraction of cases for which ground states are found in a set of runs. The first-passage-time, the time in units of Monte Carlo sweeps that a ground state is found, starting with similar random configurations, should be a good measure of the algorithms’ efficiency. We consider sample average of the first-passage-time, although the distribution of it is also very useful. The computer CPU time is another useful criterion when comparing algorithms of very different types.

In the flat-histogram and equal-hit algorithms, we can sample positive as well as negative energies uniformly. In this study, we have restricted to the negative energy part, where moves to $`E>0`$ region are rejected. We compute the average time (first-passage time) for each lattice size and given algorithm in units of sweeps ($`N=L^d`$ basic moves) to find a ground state, starting from a random configuration of equal probability of spin up and spin down. For two-dimensional $`\pm J`$ spin glass, we determine the first-passage time $`t_g`$ by comparing the current value of energy with the exact value of ground state energy, obtained from the Spin Glass Server . Thus the results are unbiased. The average first-passage time $`t_g`$ for the two-dimensional Ising spin-glass model is shown in Fig. 1. Over $`10^3`$ realization of random coupling samples are used for averaging for each algorithm and size.

Since the ground state energies are usually not known in three dimensions, we consider instead the time for finding the lowest energy for a fixed amount of sweeps $`t`$, averaged over the coupling constants with the constraint $`J_{ij}=0`$. For any fixed running length $`t`$, results obtained are only a lower bound for $`t_g`$. We consider run lengths of $`10^4`$, $`10^5`$, $`10^6`$, etc, until the first-passage time converges for large $`t`$. This limiting time $`t_g`$ is reported for the three-dimensional Ising spin-glass model in Fig. 2.

To compare the efficiency of the four algorithms, the actual CPU times are also an important factor. For our implementation, it turns out that the optimized EO, standard EO, and $`N`$-fold way equal-hit all have about the same speed at 6 microsecond per spin flip on a 700 MHz Pentium, while the single-spin-flip flat-histogram algorithm takes 3 microsecond. There are several important features in this comparison, see Fig. 1. All of them have a first-passage time that grows rapidly with sizes. With the exception of the standard EO, they nearly have the same slope of about 6 on a double logarithmic scale. It is also interesting to compare the first-passage time with that of equilibrium tunneling time reported in ref. . EO gives excellent performance for small to moderate size systems. However, for large sizes, equal-hit is as good as EO, or even better. Flat-histogram is worse by some constant factor. On the other hand, the performance of the standard EO at $`\tau =1`$ is rather poor. This shows that the results of EO is rather sensitive to the value of $`\tau `$.

Another very interesting aspect of Fig. 1 is that the curves all look linear in the semi-logarithmic scale. This implies that $`t_g\mathrm{exp}(cL)`$ for some constant $`c1`$, not a power law in $`L`$. Thus, all of the algorithms are asymptotically inefficient. It would be very interesting if this numerical observation can be supported by some argument. Similar results for the three-dimensional Ising spin glass is presented in Fig. 2.

In Table 1, we report some typical data for the average first-passage time $`t_g`$, energy per site, length of the run, and number of random samples for four algorithms for the three-dimensional spin glasses. Since we use the same set of samples with the four algorithms, a lower energy indicates a better performance. The data show that equal-hit is comparable to EO at optimal $`\tau `$.

## 4 Turning EO into an equilibrium algorithm

The flat-histogram and equal-hit algorithms can be used for equilibrium simulation. With the help of counting the number of potential moves, $`N_k`$, a basic requirement for obtaining equilibrium property of the simulated model is the microcanonical property. Using the broad histogram equation ,

$$n(E)N_k(\sigma )_E=n(E^{})N_{Zk}(\sigma ^{})_E^{},k=\left((E^{}E)/J+2Z\right)/4,$$
(13)
we can obtain the density of states $`n(E)`$ of energy $`E`$, thus the equilibrium thermodynamic quantities, including free energy.

Unfortunately, the microcanonical property that the probability distribution $`P(\sigma )`$ of the configurations is a function of energy $`E`$ only is strongly violated in EO. The probabilities of the ground states cluster into groups, rather than uniformly distributed. Numerical tests show that $`P(\sigma )`$ is a function of $`E`$, $`N_k`$, as well as additional unknown parameters. To correct this problem, we introduce a rejection step in the EO algorithm, as follows:

$$W(\sigma \sigma ^{})=\delta _N(\sigma ,\sigma ^{})P_k\frac{1}{N_k}a(\sigma \sigma ^{}).$$
(14)
The acceptance rate $`a`$ is determined by imposing a detailed balance with an unknown probability distribution $`f\left(E(\sigma )\right)`$,

$$f\left(E(\sigma )\right)W(\sigma \sigma ^{})=f\left(E(\sigma ^{})\right)W(\sigma ^{}\sigma ).$$
(15)
This gives an equation for the rate $`a`$:

$$f(E)P_k\frac{1}{N_k(\sigma )}a(\sigma \sigma ^{})=f(E^{})P_{Zk}^{}\frac{1}{N_{Zk}(\sigma ^{})}a(\sigma ^{}\sigma ).$$
(16)
The prime on $`P^{}`$ indicates that it is a set of $`P`$ values calculated from the state $`\sigma ^{}`$. A solution to this equation is a Metropolis-type choice:

$$a(\sigma \sigma ^{})=\mathrm{min}(1,\frac{f(E^{})P_{Zk}^{}/N_{Zk}(\sigma ^{})}{f(E)P_k/N_k(\sigma )}).$$
(17)
To implement this, we need a two-pass simulation. The first pass determines the function $`f(E)`$. The procedure is by no means unique. Here, we collect histogram of energy $`H(E)`$ as well as statistics for $`N_k_E`$ from an incorrect simulate of the original EO. Then we determine an approximate density of states with the help of the broad histogram equation, Eq. (13). The function $`f`$ is computed from $`f(E)=H(E)/n(E)`$. The EO with rejection is implemented in the second pass. The above procedure should be applicable for any model and any optimization algorithms that has a ‘steady state’. If the microcanonical property is only slightly violated, it will give a correct equilibrium algorithm with $`a`$ nearly equal to 1. Thus, we hope to have a method that is efficient for optimization, and yet at the same time, give correct equilibrium statistics.

Indeed, with the above method the microcanonical property is restored. Due to the rejection step, the dynamics is slightly changed. A consequence is that the histogram in the second pass shifted towards high energy side, thus the efficiency of the original EO is lost.

## 5 Conclusion

From this study, we show that equal-hit algorithm is an excellent candidate for ground state search. At the same time, it also offers the possibility for equilibrium calculations, such as the computation of the ground state entropy. We also show how optimization algorithms like EO can be turned into equilibrium algorithms by introducing a rejection step. All the algorithms studied here give rather rapid increase of $`t_g`$ with sizes, thus it is important and challenging to find algorithms that reduce this growth. Perhaps, algorithms based on single-spin-flip have their fundamental limitations.

## Acknowledgements

J.-S. W. thanks the hospitality of Tokyo Metropolitan University during part of his sabbatical leave stay. We also thank N. Kawashima and K. Chen for discussions. We thank M. Iwamatsu for drawing our attention to EO algorithm.


text:
# Scaling Cosmology

## 1 Introduction

As is widely known, current observational evidence heavily favors an accelerating and spatially flat Friedmann-Lemaître-Robertson-Walker universe (for a pedagogical short update see ). Since normal matter fulfils the strong energy condition and cannot drive cosmic acceleration, recourse is often made either to a small cosmological constant ($`\mathrm{\Lambda }`$CDM model) or to an almost evenly distributed source of energy called “dark energy” or “quintessence” with equation of state $`p_X=w_X\rho _X`$ where $`1w_X<0`$, such that it makes the pressure negative enough to render the deceleration parameter negative (see e.g. ). (Obviously, the quantity $`w_X`$ depends on the particular form assumed by the potential of the self-interacting quintessence scalar field). Since cold dark matter (i.e., dust) and quintessence decay with the expansion at different rates the question arises: “why the ratio between CDM and quintessence energies should be of the same order today?”. In other words, “where the relationship $`(\rho _M/\rho _X)_0=𝒪(1)`$ comes from?” This is in essence the “coincidence problem” . As usual, the zero subindex means present time.

Certain kind of models invoke that both components (dark matter and dark energy) may not be separately conserved due to some (unknown) coupling between each other. This proposal has been explored ,  and looks promising as a suitable mutual interaction can make both components redshift coherently. However, because of problems of their own (as the inability to recover the dust era when going back in time) neither of these proposals, as they stand, can be regarded as the final answer.

As suggested by Dalal et al. it seems rather more advisable to use the scanty observational information we possess to constrain the quintessence field from a minimum of theoretical input than trying to get a detailed fit to these data from any specific potential . These authors introduced a generalized class of dark energy models characterized by a non-canonical scaling of the ratio of the densities of dark matter and dark energy with the scale factor of the Robertson-Walker metric. They suggest a phenomenological form

$$\frac{\rho _M}{\rho _X}a^\xi $$
(1)
for the ratio of the dark matter density $`\rho _M`$ to the density $`\rho _X`$ of the dark energy, where the scaling parameter $`\xi `$ is regarded as a new variable. For an equation of state $`p_X=\rho _X`$ of the dark energy component a value $`\xi =3`$ amounts to the $`\mathrm{\Lambda }`$CDM model. A value $`\xi =0`$ represents a stationary ratio $`\rho _M/\rho _X=\mathrm{const}`$.

If the cosmological dynamics admits a stable, stationary solution $`\rho _M/\rho _X=\mathrm{const}`$, corresponding to $`\xi =0`$ and the present universe is already close to this state, there will be no coincidence problem. Consequently, according to , the deviation of the parameter $`\xi `$ from $`\xi =0`$ quantifies the severity of the problem. But it is not only the stationary solution which deserves interest. Any solution which deviates from $`\xi =3w_X`$ represents a testable, non-standard cosmological model and any solution with a scaling parameter $`\xi <3`$ will make the coincidence problem less severe. It is therefore desirable to have a physical mechanism that could give rise to such kind of deviations from the standard dynamics.

The purpose of this paper is to show that a departure from the standard $`\xi =3w_X`$ case can be obtained if cold dark matter and quintessential dark energy are no longer assumed to be separately conserved. More precisely, we shall demonstrate that a suitable interaction between dark matter and dark energy is able to produce any desired scaling. The specific parameter choice $`w_X=3`$ and $`\xi =1`$ is used to establish an exactly solvable toy model for a non-standard cosmological dynamics. Upcoming observations which will constrain cosmological models in a $`\xi w_X`$ plane, as discusseed in , will also put limits on such type of interactions.

## 2 Scaling solutions

We investigate a two–component system of cold dark matter (subindex $`M`$) and dark energy (subindex $`X`$) where

$$\rho =\rho _M+\rho _X\mathrm{and}p=p_M+p_X$$
(2)
are the total energy density and the total pressure, respectively. The components are assumed to possess the equations of state

$$p_M\rho _M\mathrm{and}p_X=w_X\rho _X.$$
(3)
We admit interactions between both components according to

$$\dot{\rho }_M+3H\rho _M=Q$$
(4)
and

$$\dot{\rho }_X+3H\left(1+w_X\right)\rho _X=Q,$$
(5)
where the coupling term $`Q`$ is to be determined below. It is convenient to introduce the quantities $`\mathrm{\Pi }_M`$ and $`\mathrm{\Pi }_X`$ by

$$Q3H\mathrm{\Pi }_M3H\mathrm{\Pi }_X,$$
(6)
with the help of which we can write ($`A`$ = M, X)

$$\dot{\rho }_A+3H\left(\rho _A+P_A\right)=0,P_A=p_A+\mathrm{\Pi }_A.$$
(7)
The coupling is then included via $`\mathrm{\Pi }_M=\mathrm{\Pi }_X`$.

To derive a specific expression for the interaction term let us consider the time evolution of the ratio $`\rho _M/\rho _X`$,

$$\left(\frac{\rho _M}{\rho _X}\right)^{}=\frac{\rho _M}{\rho _X}\left[\frac{\dot{\rho }_M}{\rho _M}\frac{\dot{\rho }_X}{\rho _X}\right].$$
(8)
¿From Eqs. (4)-(6) we obtain

$$\left(\frac{\rho _M}{\rho _X}\right)^{}=3H\frac{\rho _M}{\rho _X}\left[w_X\frac{\rho }{\rho _M\rho _X}\mathrm{\Pi }_M\right].$$
(9)
We look for solutions with the scaling behavior

$$\frac{\rho _M}{\rho _X}=r\left(\frac{a_0}{a}\right)^\xi .$$
(10)
Here, $`r`$ denotes the ratio of both components at the present time , i.e., at $`a=a_0`$, and the parameter $`\xi `$ is a constant. Inserting (10) into (9) and solving for $`\mathrm{\Pi }_M`$ we find

$$\mathrm{\Pi }_M=\mathrm{\Pi }_X=\left[\frac{\xi }{3}+w_X\right]\frac{a^\xi }{a^\xi +ra_0^\xi }\rho _M=\frac{\frac{\xi }{3}+w_X}{1+r\left(1+z\right)^\xi }\rho _M,$$
(11)
where $`1+za_0/a`$. This generalizes previous investigations for the case $`\xi =0`$ . There is a transfer of energy from the scalar field to the matter, i.e., $`Q>0`$, for $`w_X+(\xi /3)<0`$.

$`\mathrm{\Pi }_M`$ and $`\mathrm{\Pi }_X`$ are the effective pressures, equivalent to those interaction between both components given by the quantity $`Q`$ in (6), which guarantee a scaling solution (10). We arrive at the conclusion that by a suitable choice of the interaction between both components we may produce any desired scaling behavior of the energy densities.
The uncoupled case corresponding to $`\mathrm{\Pi }_M=0`$ is given by

$$\frac{\xi }{3}+w_X=0.$$
(12)
The $`\mathrm{\Lambda }`$CDM model is recovered as the special case with $`w_X=1`$ and $`\xi =3`$. The interacting models are parametrized by deviations from $`\xi =3w_X`$, equivalent to deviations from $`\xi =3`$ for $`w_X=1`$. ¿From (11) we find

$$\frac{\mathrm{\Pi }_M}{\rho _M}=\frac{\frac{\xi }{3}+w_X}{1+r\left(1+z\right)^\xi },$$
(13)
and

$$\frac{\mathrm{\Pi }_X}{\rho _X}=r\left(\frac{\xi }{3}+w_X\right)\frac{\left(1+z\right)^\xi }{1+r\left(1+z\right)^\xi }.$$
(14)
In the following we shall assume $`\xi 3w_X`$ and $`\xi >0`$, i.e., we consider departures from the standard case of separately conserved quantities. For $`1z`$, i.e., when $`\rho _X\rho _M`$ (according to (10)), we have

$$\frac{|\mathrm{\Pi }_M|}{\rho _M}1(1z),$$
(15)
and

$$\frac{\mathrm{\Pi }_X}{\rho _X}=\frac{\xi }{3}+w_X(1z).$$
(16)
While the ratio $`\mathrm{\Pi }_X/\rho _X`$ is constant in this limit and may be of the order of unity, the amount of $`\mathrm{\Pi }_M`$ is much smaller than $`\rho _M`$. For $`\frac{\xi }{3}+w_X<0`$ the $`X`$ component (that is dynamically unimportant for large $`z`$) looses energy which is transferred to the matter. While $`\mathrm{\Pi }_X`$ may be of the order of $`\rho _X`$, the quantity $`|\mathrm{\Pi }_M|`$ is negligible compared with $`\rho _M`$. Since the fractional quantities on the left-hand sides of Eqs. (15) and (16) quantify the amount of the coupling, this means, the dark matter does not feel the interaction, it is (almost) uncoupled. As the evolution proceeds, $`\mathrm{\Pi }_X/\rho _X`$ changes only slightly to a present value

$$\frac{\mathrm{\Pi }_X}{\rho _X}=\frac{r}{1+r}\left(\frac{\xi }{3}+w_X\right)\left(z=0\right),$$
(17)
whereas the corresponding ratio for the dark matter component becomes

$$\frac{\mathrm{\Pi }_M}{\rho _M}=\frac{\frac{\xi }{3}+w_X}{1+r}\left(z=0\right).$$
(18)
The point is that now the ratio $`|\mathrm{\Pi }_M|/\rho _M`$ may also be of the order of unity, i.e., the dark matter fluid feels the coupling as well. As far as the dark matter is concerned, the interaction has been switched on during the cosmic evolution. For $`\xi =0`$ we recover the relations of the previously discussed stationary solution . In the latter case the interaction does not depend on $`z`$.

Using the source terms corresponding to (6) in the balances (7), the latter can be integrated. For the matter energy density we find

$$\rho _M=\rho _M\left(a_0\right)\left[1+z\right]^{3\left(1+w_X\right)+\xi }\left[\frac{1+r\left(1+z\right)^\xi }{1+r}\right]^{1\frac{3w_X}{\xi }}.$$
(19)
The total energy density becomes

$$\rho =\rho _0\left[1+z\right]^{3\left(1+w_X\right)}\left[\frac{1+r\left(1+z\right)^\xi }{1+r}\right]^{\frac{3w_X}{\xi }},$$
(20)
where

$$\rho _0=\frac{r+1}{r}\rho _M\left(a_0\right).$$
(21)
Restricting ourselves to a universe with spatially flat sections, we obtain for the Hubble rate

$$H=\sqrt{\frac{8\pi G}{3}\rho _0}\left(1+z\right)^{\frac{3}{2}\left(1+w_X\right)}\left[\frac{1+r\left(1+z\right)^\xi }{1+r}\right]^{\frac{3w_X}{2\xi }}.$$
(22)
Likewise, the deceleration parameter $`q=\ddot{a}/\left(aH^2\right)`$ can be expressed as

$$q=\frac{1}{2}\frac{1+3w_X+r\left(1+z\right)^\xi }{1+r\left(1+z\right)^\xi }.$$
(23)
We find accelerated expansion for

$$1+3w_X+r\left(z+1\right)^\xi <0.$$
(24)
The redshift $`z_{acc}`$ at which the acceleration starts, is determined by

$$1+3w_X+r\left(z_{acc}+1\right)^\xi =0.$$
(25)
For $`w_X=1`$ (cosmological constant) we have accelerated expansion for

$$\left(z+1\right)^\xi <\frac{2}{r}z_{acc}=\left(\frac{2}{r}\right)^{\frac{1}{\xi }}1.$$
(26)
Both $`H`$ and $`q`$ depend on the set of parameters $`r`$, $`w_x`$, and $`\xi `$. The luminosity distance

$$d_\mathrm{L}=\left(1+z\right)_0^z\frac{\text{d}z}{H\left(z\right)},$$
(27)
as well as the angular distance $`d_\mathrm{A}(z)=(1+z)^2d_\mathrm{L}(z)`$, can be expressed in terms of them. The corresponding effective magnitude is

$$m_B^{eff}=M_{Bf}+5\mathrm{log}(H_0d_L),$$

where we have chosen $`M_{Bf}=3.4`$. Figure 1 depicts the effective magnitude $`m_B^{eff}`$ vs $`z`$ for $`\xi =1`$ and different values of $`w_X`$ with $`r=3/7`$. For $`\xi =3`$ almost the same figure would appear as for $`\xi =1`$. Figure 2 shows the magnitude differences for various combinations of $`\xi `$ and $`w_X`$. It becomes apparent that a much richer set of data (hopefully to be provided by the SNAP satellite) will be needed to discriminate between these models.

## 3 Special case of a scaling cosmology

The relations of the previous section considerably simplify for the special case $`\xi =1`$, $`w_X=1`$. The total energy density (20) reduces to

$$\rho =\frac{\rho _0}{\left(1+r\right)^3}\left[1+r\frac{a_0}{a}\right]^3.$$
(28)
For early times, $`aa_0`$, the energy density $`\rho `$ redshifts as dust, while for late times, $`a_0a`$, it tends to a constant value.
Likewise, the expressions for the components are

$$\rho _M=\frac{\rho _{M,0}}{\left(1+r\right)^2}\frac{a_0}{a}\left[1+r\frac{a_0}{a}\right]^2,\rho _X=\frac{\rho _{M,0}}{r\left(1+r\right)^2}\left[1+r\frac{a_0}{a}\right]^2.$$
(29)
Furthermore, the ratios which quantify the interactions among dark matter and dark energy become

$$\frac{\mathrm{\Pi }_M}{\rho _M}=\frac{2}{3}\frac{1}{1+r\frac{a_0}{a}}\text{and}\frac{\mathrm{\Pi }_X}{\rho _X}=\frac{2}{3}\frac{r\frac{a_0}{a}}{1+r\frac{a_0}{a}}.$$
(30)
The effective equation of state for the $`X`$ component is

$$P_X=\frac{1+\frac{1}{3}r\frac{a_0}{a}}{1+r\frac{a_0}{a}}\rho _X.$$
(31)
Assuming again $`r=3/7`$, this corresponds to a change from

$$P_X\frac{1}{3}\rho _X\left(z1\right)$$
(32)
at early times to

$$P_X\frac{4}{5}\rho _X\left(z=0\right)$$
(33)
at the present epoch. The effective pressure $`P_M`$ of the matter component changes from a negligible value at $`z1`$ to the present value

$$P_M=\mathrm{\Pi }_M\frac{7}{15}\rho _M\left(z=0\right).$$
(34)
The interaction has the effect that both components have a negative effective pressure.

The expression (28) for $`\rho `$ has to be contrasted with the energy density

$$\rho _{(\mathrm{\Lambda }CDM)}=\frac{\rho _0}{\left(1+r\right)}\left[1+r\left(\frac{a_0}{a}\right)^3\right]$$
(35)
for the $`\mathrm{\Lambda }`$CDM model. The sum of different powers in the latter is replaced by the power of a sum in our present model. The interaction in our model makes $`\rho _M`$ decay at a lower rate than in the uncoupled case. The dark energy density $`\rho _X`$, on the other hand, which would remain constant without interaction, decays as well as a consequence of the transfer of energy to the matter component. This feature is familiar from decaying cosmological constant models (see ).

The solution of the Friedmann equation with $`\rho `$ from (28) is

$`{\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{1}{3}}{\displaystyle \frac{\kappa \rho _0}{\left(1+r\right)^3}}\right]^{1/2}\left(tt_0\right)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{1+r}}}\left[1x\sqrt{{\displaystyle \frac{1+r}{1+x}}}\right]`$ (36)
$`+\mathrm{ln}\left\{{\displaystyle \frac{\sqrt{r}}{1+\sqrt{1+r}}}\left[x+\sqrt{x^2+1}\right]\right\},`$
where $`x=\sqrt{\frac{a}{ra_0}}`$. In the limit $`aa_0`$ we consistently recover the dust behavior $`at^{2/3}`$, while the scale factor approaches an exponential growth for $`aa_0`$. For the present special case $`w_X=1`$, $`\xi =1`$ the integral in (27) may be performed explicitly and yields

$$d_\mathrm{L}=\frac{2}{H_0}\left(1+z\right)\frac{1+r}{r}\left[1\frac{1}{\sqrt{1+\frac{rz}{1+r}}}\right].$$
(37)
For small redshifts we obtain, up to third order in $`z`$,

$$H_0^1d_\mathrm{L}z\left[1+\frac{1+\frac{r}{4}}{1+r}z\frac{1}{8}\frac{r}{\left(1+r\right)^2}\left(6+r\right)z^2\right].$$
(38)
Up to second order this expression coincides with the corresponding result for the $`\mathrm{\Lambda }`$CDM model. Differences occur only in the $`z^3`$ term. Here, the factor $`\left(6+r\right)`$ in the last expression replaces the factor $`\left(10+r\right)`$ of the $`\mathrm{\Lambda }`$CDM universe. The presently available SNIa data cannot discriminate between both models. Our model shares the merits of the $`\mathrm{\Lambda }`$CDM model but at the same time alleviates the coincidence problem.

## 4 Conclusions

Scaling solutions of the type $`\rho _M/\rho _X=r(a_0/a)^\xi `$ seem to be promising tools to deeper analyze the relationship between the two forms of energy dominating the current evolution of the universe, namely, dark matter and dark energy. We showed that a suitably chosen interaction between them can lead to any scaling behavior of the mentioned form. In the specific case $`w_X=1,\xi =1`$ the dynamics can be analytically integrated. Since the luminosity distance in this model differs from that of the $`\mathrm{\Lambda }`$CDM model only in third order in the redshift parameter $`z`$, it fits the present observations as well as the $`\mathrm{\Lambda }`$CDM model does. Further, the Einstein-de Sitter expansion law for a dust universe is recovered for large redshifts. On the other hand, the coincidence problem, although not solved, is less severe than for the $`\mathrm{\Lambda }`$CDM universe, which can be traced back to a continuous transfer of energy from the $`X`$ component to the CDM fluid.

Wile the available observational data are insufficient to discriminate between the models, it is to be expected that the SNAP satellite will provide us with a wealth of high redshift supernovae data able to do the job. Likewise, complementary observations regarding the angular distance between quasar pairs  and the evolution of cluster abundances  will further constraint the set of parameters entering the scaling models.

## Acknowledgements

The authors are indebted to David Rapetti for helping us in preparing the figures. This work was partially support by the NATO grant PST. CLG.977973, and the Spanish Ministry of Science and Technology under grant BFM 2000-C-03-01 and 2000-1322.


text:
# AGN fueling: the observational point of view

## 1. Fueling mechanisms

The essential issue to fuel gas into nuclei is to find efficient processes to transfer angular momentum. Since viscous torques are not efficient except very close to the accretion disk itself, the principal actors here are gravity torques due to non-axisymmetric structures:
1. Bars, spirals, m=1 waves
2. Galaxy interactions and mergers
These mechanisms are indeed observed when the AGN activity is very strong and accretion rate exceptional ($`>`$ 10 M/yr). For milder activities, it is difficult to see clear correlation with bars or companions. Nuclear starbursts and AGN have similar fueling requirements: they compete for gas supply and occur often simultaneously. The nuclear activity is however much more efficient to radiate energy and the fraction of energy contributed by the AGN (with respect to starbursts) increases with total luminosity.

## 2. Correlation with bars

### 2.1. Primary bars

Many studies have addressed the role of bars in nuclear activity (Simkin et al 1980, Dahari 1984, Moles et al 1995). Weak correlations are found in general, and it is often difficult to discriminate the role of biases due to morphological types, environment, etc.. Bars are more easily detected in the old stellar component, and in near infrared surveys, free of dust extinction: many more barred galaxies are discovered then, but Seyferts are not privileged (McLeod & Rieke 1995, Mulchaey & Regan 1997). The bottom line of most surveys is that there is no difference of bar frequency between AGN and non-active galaxies (see fig. 1).

A recent study comparing the CfA sample of Seyferts to a control sample finds slightly more bars in Seyferts (Peletier et al 1999, Knapen et al 2000): 79% + 7.5% of Seyferts are barred, while 59% + 9% control galaxies are barred. Paradoxically, there is also a lower fraction of strong bars in Seyferts with respect to the control sample (Shlosman et al 2000). This could be the consequence of bar destruction by central masses. As expected from gravity torques and consequent radial gas inflow, barred galaxies do possess larger molecular gas concentrations, revealed by CO observations (Sakamoto et al 1999). This gas inflow produces more nuclear starbursts, which correspond to conspicuous circumnuclear HII hot spots (Sersic 1963).

### 2.2. Decoupled stellar nuclear disks

However, the gas infall may stop at the inner Lindblad resonance of the large bar, and thus be prevented to fuel the nucleus. Are then nuclear bars the final processes to fuel the central regions? With high spatial resolution, it is now possible to investigate morphological features inside the central kiloparsec in galaxies. The first studies do not find more nuclear bars in Seyferts (e.g. Regan & Mulchaey 1999); the frequency of nuclear bars is 20-30%, the same in the control sample and for the AGN samples. Laine et al (2002) find a larger fraction of nested bars in Seyfert2 with respect to Seyfert1.

A striking feature is the high frequency of nuclear spirals discovered (Martini & Pogge 1999, Pogge & Martini 2002). The frequent presence of nuclear disks (dust, gas and young stars) point towards the action of bars and interactions, to produce the gravity torques and matter inflow. Emsellem et al (2001) discovered cold stellar nuclear disks in double-barred Seyfert galaxies, related to this inflow process. The observed stellar velocity dispersion reveals a characteristic drop in the center, interpreted in terms of recent nuclear star formation, as supported by numerical simulations (Wozniak et al 2002).

### 2.3. Why is there no correlation?

Many arguments could be advanced to explain the absence of strong correlation between bars and AGNs (e.g. Combes 2001). First, even if the bar produces gas inflow, there must be a sufficiently massive black hole at the nucleus for an AGN to appear, and this is not everywhere, in particular in view of the mass correlation between black holes and bulges: the condition is not fulfilled in late-type galaxies. Second, there could be a stellar bar, but with insufficient gas content to produce a strong enough fueling. Third, the time-scales for bar formation (and destruction), the gas fueling, and the subsequent nuclear activity could be different, and the phenomena could occur at slightly different epochs.

An interesting feature relating bars and nuclear activity is that Seyferts have more outer rings (by a factor 3-4) than non-active galaxies; these outer rings are thought to be the vestige of previous bars (Hunt & Malkan 1999). In their extended 12 microns sample of 891 galaxies, Hunt & Malkan (1999) find that 30% are AGN, 25% are interacting (perturbed). While HII/Starburst galaxies have more bars, AGN (LINERs and Seyferts) have no more bars but rings. LINERs have more inner rings, while Seyferts have 3-4 times more outer rings (cf fig. 2).

## 3. Statistics about AGN

AGN tend to lie in early-type galaxies (Terlevich et al 1987, Moles et al 1995). In an optical spectroscopic survey of 486 nearby galaxies, Ho et al (1997) detected 420 emission-lines nuclei (86% detection rate). Half of these objects can be classified as HII or star-forming nuclei, and half as some kind of AGN: Seyfert, LINERs and transition objects LINER/HII. A signature of Broad Line Region is found in 20% of the AGN, while Seyfert nuclei reside in about 10% of all galaxies. AGNs are found predominantly in luminous, early-type galaxies, while HII nuclei are in less luminous late-type objects (see fig. 3). The relation between bulge mass and black hole mass might be the reason why AGN are favored in early-types. Also the concentration of mass in the bulge favors the presence of an ILR, then nuclear bar and leads to easiest fueling.

In summary, the following points may explain that there is no significant correlation between bars and nuclear activity:

* to have an AGN a bar is not sufficient, there must already exist a massive black hole (which favors early-type objects)
* the fueling requires a large central mass concentration
* different time-scales between the bar growth and AGN activity cycles
* there exist other fueling mechanisms like dense nuclear star clusters (accounting for a time delay between the bar and fueling)

## 4. Galaxy interactions and mergers

Again, a strong correlation appears between AGN and the presence of tidal interactions only for very luminous cases (and accretion rates a few M/yr). QSOs appear to interact with companions more than usual (Hutchings & Morris 1995, Bahcall et al 1997). Seyferts can be fueled by 0.01 M/yr (like NGC 1068 for instance), and do not require external interactions.

The role of the interactions is to trigger gravitational instabilities in the galaxy disks, and the final fueling mechanism is a bar perturbation. Again the fueling is favored for early-type objects. The stability of the galaxies against perturbations depends essentially on the bulge/disk ratio (e.g. Mihos et al 1996), and the late-type objects are more susceptible to experience violent flows leading to starbursts. The fueling of AGN awaits the building of large mass concentration, the presence of ILRs, and an efficient fueling of the very center.

### 4.1. Frequency of companions

Early studies (Dahari 1984, 1985) found a correlation between AGN and the presence of companions: 15% of Seyferts have companions, compared to 3% of normal galaxies. While galaxies with companions present no more H$`\alpha `$ or radio power, they are more infra-red bright (McKenty et al 1989). Keel et al (1985) remark the presence of 25% of Seyfert in the close pairs of Arp Atlas, while the frequency of Seyfert is only 5% in a control sample. But the selection of the brightest objects in Arp Atlas might bias the conclusions. Other studies (Bushouse 1986, DeRobertis et al 1996) find a deficiency of Seyferts in interacting galaxies. A paradoxical result is the large fraction of Seyfert in the low surface brightness (LSB) galaxies (Sprayberry et al 1995), since they are unevolved objects generally in isolated environment. But the involved samples are small.

In a recent study, Schmitt (2001) does not find any statistical difference in companion frequency for systems of the various activity types. It has been claimed that Seyfert 2 have a larger number of companions than Seyfert 1 (Laurikanen & Salo 1995, Dultzin-Hacyan et al 1999), but this is controversial (Schmitt et al 2001). There could be an indirect link, some Seyfert 2 have an UV excess due to a starburst (Cid Fernandes et al 1998), and companions enhance the star formation. Narrow-line Seyf1 also do not reveal more companions (Krongold et al 2001). The percentage of galaxies with companions may sometimes appear higher for LINERs, or absorption-line galaxies with respect to Seyfert and H II galaxies. But when only galaxies of similar morphological types are considered, this difference vanishes. This suggests that the difference is due to a morphology-density effect (Schmitt 2001).

### 4.2. Compact groups

Shimada et al (2000) find no statistical difference between compact group galaxies and field galaxies, as far as AGN is concerned, and conclude that interactions do not trigger AGN activity nor starbursts. Wu et al (1998) draw the opposite conclusion, but for very luminous infrared galaxies. From their recent survey of 193 galaxies in 49 compact groups, Coziol et al (2000) confirm that AGN are located in the most luminous galaxies; AGN prefer the early-type galaxies and starbursts the late-types. AGN is the most frequent (41%) type of activity; when counting also the Starburst/AGN, the frequency rises to about 70%.

### 4.3. Radio Galaxies

Here also, interactions are obvious in the more powerful objects. FR-I have low-power, radio jets declining with radius; generally hosted by elliptical galaxies, they have rare interactions or tidal tails (less than 10%, Smith & Heckman 1989). FR-II are high power, double-lobe sources like Cygnus-A, strong in their lobe extremities: they show a high percentage of interactions from 32 to 100% (Hutchings 1987, Yates et al 1989). Most of them have two companions, and many star-forming regions.

Radio-loud QSOs have 4-5 times more neighbours, while radio-quiet 2 times more than a control sample. QSO are morphologically perturbed in 35-55% of the cases, radio-loud QSO in 70-80% (Disney et al 1995, Kirhakos et al 1999).

### 4.4. Starburst - AGN Connection

The role of gravitational interactions in luminous and ultra-luminous infrared galaxies have been established (Sanders & Mirabel 1996). In lower luminosity galaxies, only weak correlations are found (Krongold et al 2001). It appears that AGN and starbursts are two different phases of the same evolutionary process. In general Seyf2 have the most frequent starbursts (Levenson et al 2001).

Active galaxies have circumnuclear gas, it is correlated with star formation, and interactions (Storchi-Bergmann et al 2001) and with inner Hubble types from Malkan et al (1998). High spatial resolution in the center of galaxies have emphasized the differences between Seyf1 and Seyf2, that do not have the same hosts statistically. Seyferts 2 possess more circumnuclear dust, and their “inner” morphological types are shifted towards later type than the Seyferts 1. Some of the dust obscuration characteristic of Seyf2 could come from this galactic circumnuclear dust more than from an accreting torus (Malkan et al. 1998).

## 5. Orientation of accretion disks

To better understand the fueling mechanisms, it might be interesting to investigate the orientation of accretion disks observationally, with respect to the orientation of the galaxies themselves. Since the accretion disks are too small to be observed directly, a common tracer are radio jets, that are perpendicular to the accretion disks, and give their orientation. Their orientation has been compared with the rotation axis of the spiral galaxy host, but random orientations have been found (Clarke et al, 1998, Nagar & Wilson 1999, Kinney et al 2000).

Radio galaxies (in general of elliptical type) have been found to have jets perpendicular to their dust lanes (Kotanyi & Ekers 1979, van Dokkum & Franx 1995, de Koff et al 2000, for FRI only). But this is put into question in a recent study of 20 radio galaxies (Schmitt et al 2002): jets appear to have rather random orientation between 55 to 77 from the galaxy disks, while they appear to avoid an orientation too close to the disks. They are then “only roughly” perpendicular to the disks, but do not align strictly with the disk axis. This might not indicate a random orientation of the fueling through mergers and gas accretion, but this misalignement could be due to warping instabilities of the accretion disks, or warping of galaxy disks.

## 6. Conclusions

There is no correlation between AGN activity and presence of bars or interactions, except for the high luminosity objects. Many of the correlations claimed come from morphological type mismatch between AGN and control sample. AGN are predominantly found in early-type galaxies, while starbursts are found in late-types. But precisely this has to be taken into account! Since the morphology of galaxies change during their evolution.

A large set of observations compared to numerical simulations have shown that bars are only transient in galaxy disks, when gas is present. Bars produce radial inflows of gas able to destroy them, but they can re-form through gas accretion (Bournaud & Combes 2002). When there exists a central mass concentration (the galaxy looks more early ), it is more easy to fuel the nucleus. When gas is accreted all over the disk, starbursts are triggered, and the galaxy looks more late-type, but is not yet an AGN. A galaxy is in continuous evolution, and accretes mass all along its life. The evolution is governed by self-regulated processes, through bars and interactions.

The bulge-to-disk ratio and the gas fraction can evolve and therefore the morphological type somewhat oscillates from early to late-type, while secularly evolving to early-types. Any galaxy will be barred, active, and spend some time as an early-type or late-type; the frequency of bars in galaxies can be interpreted within this evolution (cf fig 4, Block et al 2002). It is when the gas has accumulated towards the center, when the type is early (although the “inner” type is more late) that nuclear activity is triggered. The morphological-type/AGN correlation is a consequence of this scenario.

## References

Bahcall, J. N., Kirhakos, S., Saxe, D. H., Schneider, D. P.: 1997, ApJ 479, 642

Block D., Bournaud F., Combes F. et al 2002, A&A (astro-ph/0209406)

Bournaud F., Combes F.: 2002, A&A 392, 83

Bushouse, H. A.: 1986, AJ 91, 255

Clarke, C. J., Kinney, A. L., Pringle, J. E.: 1998 ApJ 495, 189

Cid Fernandes R., Storchi-Bergmann T., Schmitt H.: 1998, MNRAS 297, 579

Combes F.: 2001, in “Starburst-AGN Connection”, eds I. Aretxaga et al World Sci.,223

Coziol R., Iovino A., de Carvalho R.R.: 2000, AJ 120, 47

Dahari, O.: 1984, AJ 89, 966

Dahari, O.: 1985, AJ 90, 1772

de Koff, S., Best, P., Baum, S. A. et al.: 2000, ApJS 129, 33

Derobertis, M., Hayhoe, K., Yee, H. K. C.: 1996 AAS 189 109.04

Disney, M. J., Boyce, P. J., Blades, J. C. et al.: 1995, Nature 376, 150

Dultzin-Hacyan D., Krongold, Y., Fuentes-Guridi, I., Marziani P.: 1999, ApJ 513, L111

Emsellem E., Greusard, D., Combes, F. et al.: 2001, A&A 368, 52

Ho L.C., Filippenko A.V., Sargent W.L.W.: 1997, ApJ 487, 568

Hunt L.K., Malkan M.A.: 1999 ApJ 516, 660

Hutchings, J. B.: 1987 ApJ 320, 522

Hutchings, J. B., Morris, S. C.: 1995, AJ 109, 928

Keel W.C., Kennicutt R.C., Hummel, E., van der Hulst J-M.: 1985, AJ 90, 708

Kinney A.L., Schmitt H.R., Clarke C.J. et al: 2000, ApJ 537, 152

Kirhakos, S., Bahcall, J. N., Schneider, D. P., Kristian, J.: 1999, ApJ 520, 67

Knapen, J. H., Shlosman, I., Peletier R.F.: 2000, ApJ 529, 93

Kotanyi C. G., Ekers, R. D.: 1979, A&A 73, L1

Krongold, Y., Dultzin-Hacyan, D., Marziani, P.: 2001, AJ 121, 702

Laine S., Shlosman I., Knapen J., Peletier R.F.: 2002, ApJ 567, 97

Laurikainen, E., Salo, H.: 1995, A&A 293, 683

Levenson N.A., Weaver K.A., Heckman T.M.: 2001, ApJ 550, 230

Malkan M.A., Varoujan G., Raymond T.: 1998, ApJS 117, 25

Martini P., Pogge R.W.: 1999, AJ 118, 2646

McKenty J.: 1989, ApJ 343, 125

Mcleod, K. K., Rieke, G. H.: 1995, ApJ 441, 96

Mihos, J.C., Hernquist, L.:1996, ApJ, 464, 641.

Moles, M., Marquez, I., Perez, E.: 1995, ApJ 438, 604

Mulchaey J.S., Regan M.W.: 1997, ApJ 482, L135

Nagar N.M., Wilson A.S.: 1999 ApJ 516, 97

Peletier, R. F., Knapen, J. H., Shlosman, I., et al.: 1999, ApJS 125, 363

Pogge R.W., Martini P.: 2002, ApJ 569, 624

Regan, M. W., Mulchaey, J. S.: 1999 AJ 117, 2676

Sakamoto, K., Okumura, S. K., Ishizuki, S., Scoville, N. Z.: 1999, ApJ 525, 691

Sanders D.B., Mirabel I.F.: 1996, ARAA 34, 749

Schmitt, H. R.: 2001, AJ 122, 2243

Schmitt, H. R., Antonucci R.R.J., Ulvestad J.S. et al.: 2001, ApJ 555, 663

Schmitt, H. R., Pringle, J. E., Clarke, C. J., Kinney, A. L.: 2002, ApJ 575, 150

Sersic J.J.: 1963, BAAA 6, 75

Shimada M., Ohyama Y., Nishiura S. et al.: 2000, AJ 119, 2664

Shlosman, I., Peletier, R. F., Knapen, J. H. : 2000, ApJ 535, L83

Simkin S.M., Su H.J., Schwarz M.P.: 1980, ApJ, 237, 404

Smith E., Heckman T.: 1989, ApJ 641, 358

Sprayberry, D., Impey, C. D., Bothun, G. D., Irwin, M. J.: 1995, AJ 109, 558

Storchi-Bergmann T., Gonzalez-Delgado R.M., Schmitt H.R. et al.: 2001, ApJ 559, 147

Terlevich, R., Melnick, J., Moles, M.: 1987, IAU 121, Kluwer, p.499

van Dokkum P.G. & Franx M.: 1995 AJ 110, 2027

Wozniak H., Combes F., Emsellem E. et al 2002, A&A in prep (and this meeting)

Wu H., Zou Z.L., Xia X.Y., Deng Z.G.: 1998, A&AS 132, 181

Yates, M. G., Miller, L., Peacock, J. A.: 1989, MNRAS 240, 129