text:
# On singular probability densities generated by extremal dynamics

## I Introduction

The Bak-Sneppen (BS) model was proposed as a possible explanation of mass extinctions observed in the fossil record , and was recently adapted to model experimental data on bacterial populations . Independent of its biological interpretation, the model has atracted much attention as a prototype of self-organized criticality (SOC) under extremal dynamics . The model has been studied through various approaches, including simulation , theoretical analysis , probabilistic analysis (run time statistics) , renormalization group , field theory  and mean-field theory . Some variants have been proposed, for example the anisotropic BS model . In this paper, we study the consequences of extremal dynamics, using mean-field theory and simulation. With this aim, we propose variants of the model and analyse how varying the updating rule affects the stationary probability density and the critical behaviour.

In the evolutionary interpretation of the BS model, each site $`i`$ represents a “niche” occupied by a single species, and bears a real-valued variable $`x_i`$ representing the “fitness” of this species. (In the present context “fitness” denotes a propensity to resist extinction: if $`x_i>x_j`$ then species $`j`$ goes extinct before $`i`$, so that $`x_i`$ may be termed a “barrier to extinction”.) At each step the site with the smallest $`x_i`$, and its nearest neighbors, are replaced with randomly chosen values. The replacement of the neighboring variables with new random values may be interpreted as a sudden unpredictable change in fitness when a nearby niche (which might have borne a predator, or a food source for the species in question), is suddenly colonized by a new species. Selection, at each step, of the global minimum of the $`\{x_i\}`$ (“extremal dynamics”) represents a highly nonlocal process, and would appear to require an external agent with complete information regarding the state of the system at each moment. Applications of the model in evolution studies are reviewed in .

In the original updating rule, neighbors are randomly affected by the extinction of an interacting species. There is no a priori reason to expect that evolution should obey this specific rule on a specific lattice. Thus, we ask: what happens if the extinction of the least adapted species favors (or prevents) the extinction of other species? Moreover, what is the signature of extremal dynamics that can, in principle, allow us to recognize it in the real world?

Due to extremal dynamics, the BS model exhibits scale invariance in the stationary state, in which several quantites display power-law behaviour . Simulations show that the stationary distribution of barriers follows a step function, being zero (in the infinite-size limit) for $`x<x^{}0.66702(8)`$ . Relaxing the extremal condition leads to a smooth probability density and loss of scale invariance .

A striking feature of the BS model is that a simple updating rule leads to a singular stationary probability density $`p(x)`$. An intriguing question therefore arises, as to how changes in the updating rule affect this density, an issue that has not, to our knowledge, been investigated previously. In this work we examine the consequences of rules in which one or more sites are updated according to $`xx^{}=f(x)`$, instead of being replaced by a random number. (Here $`f`$ maps the interval  to itself.) We find that this can provoke dramatic changes in the stationary probability density. In the extremal dynamics limit, the variant models belong to the same universality class as the original. We find that the hallmarks of extremal dynamics are that i) the stationary probability density is singular, and ii) with probability one, the extremal $`x_i`$ (i.e., the next variable to be updated) belongs to the ‘prohibited’ region in which $`p(x)=0`$. Using a two-site mean-field approximation, we also find evidence that nontrivial correlations are restricted to the prohibited region.

The balance of this paper is organized as follows. Section II introduces the models, which are then analyzed using mean-field like approaches in Sec. III. In Sec. IV we present simulation results, and summarize our conclusions in Sec. V.

## II Models

The Bak-Sneppen  model is a discrete-time Markov process on a d-dimensional lattice of $`L^d`$ sites, with periodic boundaries. At each site we define a real-valued variable $`x_i(0)`$. Initially, these variables are independently assigned random values uniform on \[0,1). At time 1, the site $`m`$ bearing the minimum of all the numbers $`\{x_i(0)\}`$ is identified, and it, along with its $`2d`$ nearest neighbors, are given new random values, again drawn independently from the interval \[0,1). (In the one dimensional case considered here this amounts to: $`x_m(1)=\eta `$, $`x_{m+1}(1)=\eta ^{}`$, and $`x_{m1}(1)=\eta ^{\prime \prime }`$, where $`\eta `$, $`\eta ^{}`$, and $`\eta ^{\prime \prime }`$ are independent and uniformly distributed on \[0,1); for $`|jm|>1`$, $`x_j(1)=x_j(0)`$.) At step 2 this process is repeated, with $`m`$ representing the site with the global minimum of the variables $`\{x_i(1)\}`$, and so on. In the random neighbor version of the model, the process is realized on a complete graph (all sites are considered neighbors); two randomly selected sites are updated in addition to the minimum $`m`$.

We now define three modified Bak-Sneppen models, that differ from the original only in the way that the barriers $`x_i`$ evolve. In one, the site $`M`$ bearing the maximum of the $`\{x_i\}`$ is replaced with a random number $`\eta `$ and the two nearest neighbors are replaced with their own square roots: $`x_M(t+1)=\eta `$ and $`x_{M\pm 1}(t+1)=\sqrt{x_{M\pm 1}(t)}`$. We shall refer to this as the ‘radical’ variant. In the second variant, the site with the maximum value is replaced by its own value squared, while its two nearest neighbors receive random numbers $`\eta `$ and $`\eta ^{}`$: $`x_M(t+1)=[x_M(t)]^2`$, $`x_{M+1}(t+1)=\eta `$, and $`x_{M1}(t+1)=\eta ^{}`$. This will be called the ‘centered square’ version. Finally, we define a ‘peripheral square’ variant, in which one of the nearest neighbors of $`M`$ is squared, while $`M`$ and its other neighbor are replaced with random numbers: $`x_M^{}(t+1)=[x_M^{}(t)]^2`$, $`x_M(t+1)=\eta `$ and $`x_{M^{\prime \prime }}(t+1)=\eta ^{}`$. (Here $`M^{}=M+\sigma `$ and $`M^{\prime \prime }=M\sigma `$, where $`\sigma `$ is a random variable that assumes values of +1 and -1 with equal probability.)

The motivation for studying these variants is twofold. First, the Bak-Sneppen model is notable for exhibiting a singular stationary probability density, and it is of interest to examine the effect of changes in the dynamical rule on this density and on the critical behavior. If we introduce a deterministic function $`f(x)`$ as part of the updating rule, it is desirable that $`f`$ map the interval  to itself, making functions of the form $`f(x)=x^\alpha `$ a natural choice. In this context we note that the variants feature what may be called ‘migration’, that is, the systematic movement of certain variables $`x_i`$ within the interval $`[0,1)`$. In the radical variant the migration is from the populated region (smaller $`x`$) to the ‘excluded’ region (larger $`x`$), whereas in the square versions migration occurs in the opposite direction.

Secondly, the variants admit interpretation as evolutionary processes. (In the modified models, we have for convenience defined the site with the maximum variable as the most vulnerable, so that small $`x_i`$ now corresponds to high fitness.) In the radical variant, replacement of the least-fit species provokes a reduction in the fitness of its neighbors, without leading to their immediate extinction. Thus, some memory of the fitness of the neighboring species is retained. The radical variant therefore seems a plausible modification of the original model, in the biological context. In the peripheral square variant the extinction of the least-fit species provokes extinction of one neighbor, and increased fitness of the other. Finally, in the centered-square variant, $`x_Mx_M^2`$ represents an increase in the fitness of the least viable element of the system, while its neighbors go extinct.

## III Mean-field theory

### A Original Model

We develop mean-field approximations for the original and modified models, along the lines of Refs.  and . To begin, we relax the extremal condition introducing a flipping rate of $`\mathrm{\Gamma }e^{\beta x_i}`$ at site $`i`$, where $`\mathrm{\Gamma }^1`$ is a characteristic time, irrelevant to stationary properties, and which we set equal to one ($`\mathrm{\Gamma }=1`$). Call this regularized system the ‘finite-temperature’ model. The extremal dynamics of the original model is recovered in the zero-temperature limit, $`\beta \mathrm{}`$.

Consider the probability density $`p(x)`$. In the finite-temperature version of the original model, $`p(x)`$ satisfies

$$\frac{dp(x,t)}{dt}=e^{\beta x}p(x,t)2_0^1e^{\beta y}p(x,y,t)𝑑y+3_0^1e^{\beta y}p(y,t)𝑑y,$$
(1)
where $`p(x,y,t)`$ is the joint density for a pair of nearest-neighbor sites and $`p(y,t)`$ is the one-site marginal density. Invoking the mean-field factorization $`p(x,y)=p(x)p(y)`$ (we suppress the time argument from here on), we find:

$$\frac{dp(x)}{dt}=[e^{\beta x}+2I(\beta )]p(x)+3I(\beta ),$$
(2)
where

$$I(\beta )_0^1e^{\beta y}p(y)𝑑y,$$
(3)
represents the overall flipping rate. Eq. (2) is a nonlinear equation for $`p`$, in which the density at $`x`$ is coupled to $`p`$ at all other arguments via $`I(\beta )`$.

In the stationary state we have

$$p(x)=\frac{3I}{2I+e^{\beta x}}.$$
(4)
Multiplying by $`e^{\beta x}`$ and integrating over the range of $`x`$, we find

$$I(\beta )=(e^{2\beta /3}1)/[2e^\beta (1e^{\beta /3})],$$
(5)
and thus

$$p_{st}(x)=\frac{3}{2}\frac{1e^{2\beta /3}}{1e^{2\beta /3}+e^{\beta x}(e^{\beta /3}1)}.$$
(6)
This solution is plotted for various $`\beta `$ values in Fig. 1.

In the limit $`\beta \mathrm{}`$ the solution becomes a step function:

$$p_{st}(x)=\frac{3}{2}\mathrm{\Theta }(x1/3)\mathrm{\Theta }(1x).$$
(7)
Thus the mean-field approach correctly predicts a step-function singularity for the probability density, although it places the critical barrier at $`x^{}=1/3`$, whereas it actually falls at 0.66702(8) . On the other hand, the slightest relaxation of the extremal condition destroys the singularity , since $`p(x)`$ is a smooth function for $`\beta <\mathrm{}`$. The rate of convergence to the step-function is generally exponential with $`\beta `$, away from the discontinuity. The curves for various $`\beta `$ values exhibit an approximate crossing near $`x=1/3`$. The derivative at this point however diverges only linearly with $`\beta `$: $`(dp_{st}/dx)_{x=1/3}3\beta /8`$ for large $`\beta `$.

Using Eq. (2), we find, in the limit of large $`\beta `$, that the relaxation time for a small disturbance from the stationary solution grows $`e^{\beta /3}`$. (By ‘small’ we mean $`I(\beta )e^{\beta /3}/2`$.)

The following observation will prove useful in the discussion of the modified models. If we assume that, in the limit $`\beta \mathrm{}`$, $`p_{st}(x)`$ is identically zero for $`x<x^{}`$, and that $`p_{st}C>0`$ on some interval $`[x^{},a]`$ (in other words, the density suffers a jump discontinuity at $`x^{}`$), then $`I=_0^1e^{\beta y}p(y)𝑑ye^{\beta x^{}}`$ and so $`e^{\beta x}/Ie^{\beta (xx^{})}0`$ for $`x>x^{}`$. Then Eq. (4) reduces to the step-function expression, Eq. (7). Note however that $`lim_\beta \mathrm{}e^{\beta x^{}}I(\beta )=1/2`$ not $`3/2\beta `$, as would be found by naively inserting the limiting density, Eq.(7), in Eq.(3). This means that the dominant contribution to $`I`$ is due to the interval $`[0,x^{}]`$, even in the limit $`\beta \mathrm{}`$, which is readily seen if we write Eq.(6), for large $`\beta `$, as

$$p_{st}(x)\frac{3}{2}[\mathrm{\Theta }(xx^{})+\mathrm{\Theta }(x^{}x)e^{\beta (xx^{})}].$$
(8)
In the limit $`\beta \mathrm{}`$, sites with $`x<x^{}`$ constitute a set of probability zero, but the site $`m`$ selected for extinction belongs to this set with probability one. This is a singular property of the Bak-Sneppen model in the infinite-size limit, as discussed in the next subsection. (The infinite-size limit is implicit in mean-field theory.)

Although in the modified models we are unable to find an analytical solution for finite $`\beta `$, it is possible to integrate the mean-field equation numerically. Due to the factor $`e^{\beta x}`$, for large $`\beta `$, a very small time step would be needed to avoid instability in the usual integration methods (e.g., Euler or Runge-Kutta). We circumvent this difficulty using a partial integration method . To apply this method to the MF equation for the original model, we write Eq. (2) in the form

$$\frac{dp(x,t)}{dt}=f(t)p(x,t)+g(t),$$
(9)
where $`f(t)=e^{\beta x}+2I(t)`$ and $`g(t)=3I(t)`$. The formal solution is

$$p(x,t)=\mathrm{exp}[_0^t𝑑uf(u)]\left\{p(x,0)+_0^t𝑑t^{}\mathrm{exp}\left[_0^t^{}𝑑t^{\prime \prime }f(t^{\prime \prime })\right]g(t^{})\right\}.$$
(10)
For a small time interval $`\mathrm{\Delta }t`$, we find

$`p(x,\mathrm{\Delta }t)`$ $``$ $`e^{f(0)\mathrm{\Delta }t}\left\{p(x,0)+g(0){\displaystyle _0^{\mathrm{\Delta }t}}𝑑t^{}e^{f(0)t^{}}\right\}`$ (11)
$`=`$ $`e^{f(0)\mathrm{\Delta }t}p(x,0)+{\displaystyle \frac{g(0)}{f(0)}}(1e^{f(0)\mathrm{\Delta }t}).`$ (12)
This relation can be iterated to find the evolution of $`p(x,t)`$ from a given initial distribution, which converges to the stationary density.

### B BS model on a finite complete graph

Mean-field theory is exact for the “random neighbor” model, which may also be thought of as the BS model on a complete graph, i.e., one in which all pairs of sites are neighbors. (When $`x_m`$ is updated, two of these neighbors are chosen at random for updating as well.) In this subsection we analyze the BS model with extremal dynamics on a complete graph of $`N`$ sites.

Since sites are assigned independent random numbers, the $`x_i`$ are independent, identically distributed random variables drawn from the density $`p(x,t)`$. Define the distribution function $`P(x,t)=_0^xp(y,t)𝑑y`$. The probability that the next site to be updated, $`x_m`$, lies between zero and $`x`$ is:

$$\text{Prob}[x_mx]=1\left[1P(x)\right]^N,$$
(13)
i.e., one less the probability that the minimum is larger than $`x`$. The probability that a randomly chosen neighbor has $`x_ix`$ is simply $`P(x)`$, and the probability that one of the updated sites receives a number $`x`$ is $`x`$. At each step, therefore, the expected change in the number of sites with $`x_ix`$ is $`3x2P(x)\text{Prob}[x_mx]`$, which implies

$$\frac{dP(x,t)}{dt}=\left\{1\left[1P(x)\right]^N\right\}2P(x)+3x.$$
(14)
(Here we have taken the time unit to represent $`N`$ updates.)

Letting $`Q1P`$, we have in the stationary state,

$$Q^N+2Q3(1x)=0.$$
(15)
(Note that $`Q(0)=1`$, $`Q(1)=0`$, and $`dQ/dx0`$.) Numerical solution (Fig. 2) shows that for large $`N`$, $`P_N(x)`$ approaches a singular function that is zero for $`x<1/3`$, while for $`x>1/3`$, $`P(x)=3(x1/3)/2`$. It is straightforward to show that for fixed $`x`$ and $`N\mathrm{}`$,

$$P(x)\{\begin{array}{c}1(13x)^{1/N},x<1/3\hfill \\ \\ \frac{3x1}{2}+\frac{1}{2}[\frac{3}{2}(1x)]^N,x>1/3.\hfill \end{array}$$
(16)
For $`x=1/3`$ we have $`\mathrm{ln}P\mathrm{ln}N+\mathrm{ln}\mathrm{ln}(N/2)`$ (plus terms of lower order in $`N`$) as $`N\mathrm{}`$. Of interest is the exponential convergence of $`P`$ to its limiting form for $`x>1/3`$, compared with algebraic convergence for $`x<1/3`$. Note also that $`\text{Prob}[x_m1/3]12/N`$, so that the minimum indeed belongs to the excluded region with probability one, when $`N\mathrm{}`$.

### C Pair approximation

The analysis of the finite-temperature model is readily extended to the pair level, in which one studies the evolution of the two-site joint probability density $`p(x,y)`$. Our starting point is the following exact relation, obtained using the same reasoning that led to Eq. (1):

$`{\displaystyle \frac{dp(x,y)}{dt}}`$ $`=`$ $`\left(e^{\beta x}+e^{\beta y}\right)p(x,y){\displaystyle _0^1}e^{\beta u}\left[p(x,y,u)+p(u,x,y)\right]𝑑u`$ (17)
$`+`$ $`{\displaystyle _0^1}{\displaystyle _0^1}\left(e^{\beta u}+e^{\beta v}\right)p(u,v)𝑑u𝑑v`$ (18)
$`+`$ $`{\displaystyle _0^1}{\displaystyle _0^1}\left[p(x,u,v)+p(v,u,y)\right]e^{\beta v}𝑑u𝑑v.`$ (19)
Now invoking the pair factorization,

$$p(x,y,u)\frac{p(x,y)p(y,u)}{p(y)},$$
(20)
Eq. (19) becomes

$`{\displaystyle \frac{dp(x,y)}{dt}}`$ $`=`$ $`\left[e^{\beta x}+e^{\beta y}+K(x)+K(y)\right]p(x,y)+2I`$ (21)
$`+`$ $`{\displaystyle _0^1}\left[p(x,u)+p(y,u)\right]K(u)𝑑u,`$ (22)
with

$$K(x)=J(x)/p(x),$$
(23)
$$J(x)=_0^1p(x,u)e^{\beta u}𝑑u$$
(24)
and

$$I=_0^1p(x)e^{\beta x}𝑑x=_0^1J(x)𝑑x.$$
(25)
To find the stationary solution numerically, we note that

$$p(x,y)=\frac{2I+_0^1\left[p(x,u)+p(y,u)\right]K(u)𝑑u}{e^{\beta x}+e^{\beta y}+K(x)+K(y)}.$$
(26)
Starting from an arbitrary normalized density $`p_0`$ (for example, uniform on $`[0,1]\times [0,1]`$), we generate $`p_1`$ by evaluating the r.h.s. of Eq. (26) using $`p_0`$ and normalizing the resulting expression. This procedure is then iterated until it converges to the stationary density. We find that for large $`\beta `$, the stationary marginal density approaches the step-function solution

$$p(x)=\frac{1}{1x^{}}\mathrm{\Theta }(xx^{})\mathrm{\Theta }(1x)$$
(27)
with $`x^{}0.47186`$, a considerable improvement over the site approximation. In this limit, the joint distribution is the product of two identical one-site distributions. Once again, the portion of the unit square that has probability zero is in fact responsable for all transitions. In the region $`D\{(x,y)|0<x,y<x^{}\}`$, the two variables are correlated, as shown by the nonzero correlation coefficient $`\rho \text{cov}(x,y)/\sqrt{\text{var}(x)\text{var}(y)}`$. In the pair approximation, $`\rho 0.327`$ for $`(x,y)D`$, as $`\beta \mathrm{}`$.

### D Radical Variant

In the radical model, the probability density $`p(x)=p_X(x)`$ satisfies:

$`{\displaystyle \frac{dp_X(x)}{dt}}`$ $`=`$ $`e^{\beta x}p_X(x)2{\displaystyle _0^1}e^{\beta y}p_X(x,y)𝑑y+{\displaystyle _0^1}e^{\beta y}p_X(y)𝑑y`$ (28)
$`+`$ $`2p_{X^{1/2}}(x){\displaystyle _0^1}e^{\beta y}p_X(y)𝑑y,`$ (29)
where $`p_{X^{1/2}}(x)=2xp_X(x^2)`$. With the definition $`I(\beta )_0^1e^{\beta y}p(y)𝑑y`$ and the mean-field assumption $`p(x,y)=p(x)p(y)`$, Eq. (29) reduces to

$$\frac{dp(x)}{dt}=e^{\beta x}p(x)+I(\beta )[2p(x)+1+4xp(x^2)].$$
(30)
Given the step-function form of the stationary density for the original model, it is reasonable to expect that in this case as well, $`p(x)`$ will have a jump discontinuity for $`\beta \mathrm{}`$, and be zero for $`x>x^{}`$. We can get some insight into the nature of the density as $`\beta \mathrm{}`$ by first observing that if $`p(x)C<\mathrm{}`$ as $`x1`$, then $`ICe^\beta /\beta `$ for large $`\beta `$, and therefore $`e^\beta I0`$ as $`\beta \mathrm{}`$. This in fact holds unless $`p(x)`$ has a $`\delta `$-like contribution at $`x=1`$, which is not expected since it is precisely the largest values of $`x`$ that are removed in the dynamics. Write the stationary solution to Eq. (30) as

$$p(x)=\frac{\overline{I}[1+4xp(x^2)]}{1+2\overline{I}},$$
(31)
with $`\overline{I}=e^{\beta x}I`$. This gives $`p(1)=0`$ in the limit $`\beta \mathrm{}`$. Similarly, supposing that $`xp(x^2)0`$ as $`x0`$, we have $`p(0)=1/2`$ in this limit. A similar line of reasoning can be used to show that $`dp/dx|_{x=1}=0`$ in the $`\beta \mathrm{}`$ limit.

The preceding discussion suggests that in the limit $`\beta \mathrm{}`$, $`p_{st}(x)`$ is identically zero over some finite range $`[x^{},1]`$. Assuming this to be so, we have $`Ie^{\beta x^{}}`$ and in the limit $`\beta \mathrm{}`$ the stationary solution is given by the functional equation:

$$2p(x)4xp(x^2)=1,$$
(32)
for $`0x<x^{}`$. Writing this as $`p(x)=1/2+2xp(x^2)`$, we can iterate to find $`p(x^2)=1/2+2x^2p(x^4)`$, $`p(x^4)=1/2+2x^4p(x^8)`$, $`p(x^8)=1/2+2x^8p(x^{16})`$, etc., which suggests the solution

$$p_1(x)=\frac{1}{2}+x+2x^3+4x^7+8x^{15}+\mathrm{}=\underset{i=0}{\overset{\mathrm{}}{}}2^{i1}x^{2^i1}.$$
(33)
Substituting this ‘lacunary series’ in Eq.(32), one readily verifies that $`p_1(x)`$ is a solution. Similarly, rewriting Eq.(32) as $`p(x^2)=1/4x+p(x)/2x`$, one finds $`p(x)=1/4x^{1/2}+p(x^{1/2})/2x^{1/2}`$, $`p(x^{1/2})=1/4x^{1/4}+p(x^{1/4})/2x^{1/4}`$, etc., leading to a second solution:

$$p_2(x)=\frac{1}{4x^{1/2}}\frac{1}{8x^{3/4}}\frac{1}{16x^{7/8}}\frac{1}{32x^{15/16}}\mathrm{}=\underset{i=1}{\overset{\mathrm{}}{}}2^{(i+1)}x^{(12^i)}.$$
(34)
We now search a solution of the form $`p(x)=Ap_1(x)+Bp_2(x)`$; substituting in Eq. (32), yields the condition $`A+B=1`$. This linear combination must however be normalizable. The relevant integrals are:

$$_0^x^{}p_1(x)𝑑x=\frac{1}{2}\underset{n=0}{\overset{\mathrm{}}{}}x_{}^{}{}_{}{}^{2^n}=\frac{1}{2}(x^{}+x_{}^{}{}_{}{}^{2}+x_{}^{}{}_{}{}^{4}+x_{}^{}{}_{}{}^{8}+\mathrm{}),$$
(35)
$$_0^x^{}p_2(x)𝑑x=\frac{1}{2}\underset{n=1}{\overset{\mathrm{}}{}}x_{}^{}{}_{}{}^{2^n}=\frac{1}{2}(1+x_{}^{}{}_{}{}^{1/2}+x_{}^{}{}_{}{}^{1/4}+x_{}^{}{}_{}{}^{1/8}+\mathrm{}).$$
(36)
Since $`0<x^{}<1`$, the first integral converges while the second diverges, so that $`p_2(x)`$ is not normalizable, implying $`B=0`$. A normalized, positive solution is $`p(x)=p_1(x)`$ for $`x<x^{}=0.793189`$ and $`p(x)=0`$ for $`x>x^{}`$. ($`x^{}`$ is determined by normalization.) This solution, for infinite $`\beta `$, is plotted in Fig. 3. Numerical integration of Eq. (30) through the method outlined in Sec. IIIA yields results consistent with this solution, as may again be seen in the figure. Finally, simulation of the random-neighbor version of the model yields a stationary distribution consistent with this expression (see Fig. 4).

We will now analyze the radical model on a finite complete graph and show that its probability density (equation 33) can be derived via a different path. Define the distribution function $`Q(x,t)=_x^1p(y,t)𝑑y`$. By the same reasoning developed in section III B, the probability that $`x_m`$, the next site to be updated, lies between $`x`$ and 1 is

$$\text{Prob}[x_mx]=1\left[1Q(x)\right]^N.$$
(37)
The probability that a randomly chosen neighbor has $`x_ix`$ is $`Q(x)`$, while the probability that a neighbor receives a barrier $`x`$ is $`Q(x^2)`$. The probability that $`x_m`$ receives a new value between $`x`$ and 1 is simply $`1x`$. The expected change in the number of sites with $`x_ix`$ is $`(1x)+2Q(x^2)2Q(x)\{1\left[1Q(x)\right]^N\}`$, so that

$$\frac{dQ(x,t)}{dt}=(1x)+2Q(x^2)2Q(x)\left\{1\left[1Q(x)\right]^N\right\}.$$
(38)
In the sationary state, the definition $`P1Q`$ leads to

$$x+2P(x^2)2P(x)P(x)^N=0.$$
(39)
(Note that $`P(0)=0`$ and $`P(1)=1`$). Since $`0P(x)<1`$ for $`x<x^{}`$, $`lim_N\mathrm{}P(x)^N=0`$ and the probability density in an infinite system obeys the functional equation $`x+2P(x^2)2P(x)=0`$. Using the iterative approach, we find two solutions, one divergent (which is rejected), and the finite solution

$$P(x)=\frac{1}{2}\underset{n=0}{\overset{\mathrm{}}{}}x^{2^n},$$
(40)
i.e., the integral of $`p_1(x)`$, Eq. (33). Normalization then demands that $`P=1`$ for $`xx^{}`$.

Numerical solution of equation (38) (see Fig. 5) shows that, for large $`N`$, $`P_N(x)`$ approaches the singular function described above. For fixed $`x>x_{}^{}{}_{}{}^{1/2}`$, and $`N\mathrm{}`$, we find $`P(x)1x^{1/N}`$.

### E Centered Square Variant

The probability density in the finite-temperature version of the centered square variant obeys:

$$\frac{dp_X(x)}{dt}=e^{\beta x}p_X(x)2_0^1e^{\beta y}p_X(x,y)𝑑y+2_0^1e^{\beta y}p_X(y)𝑑y+\frac{e^{\beta \sqrt{x}}}{2\sqrt{x}}p_X(\sqrt{x}).$$
(41)
Mean-field factorization leads to

$$\frac{dp(x)}{dt}=e^{\beta x}p(x)2I(\beta )[p(x)1]+\frac{e^{\beta \sqrt{x}}}{2\sqrt{x}}p(\sqrt{x}),$$
(42)
with $`I(\beta )`$ as defined in Sec. IIID. In the stationary state, we find

$$p(x)=\frac{2+\frac{1}{2\sqrt{x}}p(\sqrt{x})e^{\beta \sqrt{x}}/I}{2+e^{\beta x}/I}.$$
(43)
The hypothesis that in the limit $`\beta \mathrm{}`$ , the stationary density $`p(x)=0`$ for $`x>x^{}`$ implies $`I(\beta )e^{\beta x^{}}`$. Therefore

$$\underset{\beta \mathrm{}}{lim}\frac{e^{\beta x}}{I}=\{\begin{array}{c}0ifx<x^{}\\ \mathrm{}ifx>x^{}\end{array}$$
(44)
and

$$\underset{\beta \mathrm{}}{lim}\frac{e^{\beta \sqrt{x}}}{I}=\{\begin{array}{c}0ifx<x_{}^{}{}_{}{}^{2}\\ \mathrm{}ifx>x_{}^{}{}_{}{}^{2}\end{array}$$
(45)
This, combined with Eq. (43), implies that $`p(x)=1`$ for $`x[0,x_{}^{}{}_{}{}^{2}]`$, and that

$$p(x)=\underset{\beta \mathrm{}}{lim}\left[1+\frac{1}{4\sqrt{x}}\frac{e^{\beta \sqrt{x}}}{I}p(\sqrt{x})\right]\mathrm{for}x[x_{}^{}{}_{}{}^{2},x^{}]$$
(46)
The iterative process used in the preceding section is not useful here as it leads to pathological solutions, due to the divergent ratio $`e^{\beta \sqrt{x}}/I`$ in this interval. Observe that if, for $`x[x^{},x^{1/2}]`$, $`p(x)=g(x^2)I/e^{\beta x}`$, as $`\beta \mathrm{}`$, with $`g(x^2)`$ finite, then $`p(x)=1+g(x)/4\sqrt{x}`$ for $`x[x_{}^{}{}_{}{}^{2},x^{}]`$. (Note that this means that $`p(x)0`$ as $`\beta \mathrm{}`$ for $`x`$ in the interval $`[x^{},x_{}^{}{}_{}{}^{1/2}]`$, consistent with the hypothesis that $`p(x)0`$ for $`x>x^{}`$.) The function $`g(x)`$ is however yet to be determined. Attempting the simplest solution, $`g(x)=`$ constant, we find a surprisingly reasonable result, as shown in Fig. 6. Next, allowing $`g(x)=ax+b`$, with $`a`$ and $`b`$ constant, yields excellent agreement with the numerical integration of Eq. (42) as also shown in Fig. 6. We do not have an argument why $`g(x)`$ should take this form.

The threshold $`x^{}`$ can be determined in a simple way. Since the fraction of barriers in the interval $`[x^{},1]`$ is constant in the stationary state, the mean number of barriers removed from this interval at each time step must equal the mean number inserted. The probability that the maximum $`x_M`$ lies in $`[x^{},1]`$ is 1, while its random neighbors are certainly below $`x^{}`$. Each updated neighbor has a probability $`(1x^{})`$ of receiving a barrier in the interval $`[x^{},1]`$, while the maximum remains in $`[x^{},1]`$ with probability $`p_1`$, the probability that $`x_M[x_{}^{}{}_{}{}^{1/2},1]`$. Thus, we have

$$1=2(1x^{})+p_1.$$
(47)
This reasoning can be repeated for the intervals $`[x_{}^{}{}_{}{}^{1/2},1]`$, $`[x_{}^{}{}_{}{}^{1/4},1]`$, …, leading to

$`p_1=2(1x_{}^{}{}_{}{}^{1/2})+p_2,`$ (48)
$`p_2=2(1x_{}^{}{}_{}{}^{1/4})+p_3\mathrm{}`$ (49)
where $`p_n`$ is the probability that $`x_M[x_{}^{}{}_{}{}^{1/2^n},1]`$. Substituting this result in equation (47), we find

$$1=2(1x^{})+2(1x_{}^{}{}_{}{}^{1/2})+2(1x_{}^{}{}_{}{}^{1/4})+2(1x_{}^{}{}_{}{}^{1/8})+\mathrm{},$$
(50)
which provides $`x^{}=0.761072`$. Finally, we note that normalization implies $`_{x_{}^{}{}_{}{}^{2}}^x^{}g(x)/4\sqrt{x}𝑑x=1x^{}=0.238928`$, providing a constraint on the function $`g(x)`$.

### F Peripheral Square Variant

We now apply the mean-field analysis to the peripheral square variant. In this case, the probability density satisfies:

$`{\displaystyle \frac{dp_X(x)}{dt}}`$ $`=`$ $`e^{\beta x}p_X(x)2{\displaystyle _0^1}e^{\beta y}p_X(x,y)𝑑y+2{\displaystyle _0^1}e^{\beta y}p_X(y)𝑑y`$ (51)
$`+`$ $`p_{X^2}(x){\displaystyle _0^1}e^{\beta y}p_X(y)𝑑y,`$ (52)
where $`p_{X^2}(x)=p_X(x^{1/2})/2x^{1/2}`$. Under the mean-field factorization this reduces to

$$\frac{dp(x)}{dt}=e^{\beta x}p(x)+I(\beta )[2p(x)+2+\frac{1}{2x^{1/2}}p(x^{1/2})],$$
(53)
with $`I(\beta )`$ as given in Sec. IIIC. The hypothesis that, for $`\beta \mathrm{}`$, the stationary density $`p(x)=0`$ for $`x>x^{}`$ implies $`I(\beta )e^{\beta x^{}}`$, leading to the functional equation:

$$p(x)\frac{1}{4x^{1/2}}p(x^{1/2})=1(\beta \mathrm{}),$$
(54)
for $`x<x^{}`$.

As in the centered-square variant, the iterative method does not yield a useful solution, and we pursue a different approach. Let $`p(x)=f(x)`$ on the interval $`x_{}^{}{}_{}{}^{2}x<x^{}`$, where the function $`f(x)`$ and the constant $`x^{}`$ are yet to be determined. Using Eq. (54), we find that

$`p(x)=1+{\displaystyle \frac{f(x^{1/2})}{4x^{1/2}}}`$ $`\mathrm{for}x_{}^{}{}_{}{}^{4}x<x_{}^{}{}_{}{}^{2},`$ (55)
$`p(x)=1+{\displaystyle \frac{1}{4x^{1/2}}}+{\displaystyle \frac{f(x^{1/4})}{4^2x^{3/4}}}`$ $`\mathrm{for}x_{}^{}{}_{}{}^{8}x<x_{}^{}{}_{}{}^{4},`$ (56)
$`p(x)=1+{\displaystyle \frac{1}{4x^{1/2}}}+{\displaystyle \frac{1}{4^2x^{3/4}}}+{\displaystyle \frac{f(x^{1/8})}{4^3x^{7/8}}}`$ $`\mathrm{for}x_{}^{}{}_{}{}^{16}x<x_{}^{}{}_{}{}^{8},`$ $`\mathrm{etc}.`$ (57)
Thus we have found a family of solutions $`p(x)`$ to Eq.(54), which is quite general since $`f(x)`$ is still undetermined.

We now show that $`f(x)=1`$ on the interval $`[x^2,x^{})`$. To begin, we note that in the stationary state, Eq. (53) implies

$$\left(\frac{e^{\beta x}}{I}+2\right)p(x)=2+\frac{p(\sqrt{x})}{2\sqrt{x}}.$$
(58)
In particular, for $`x=1`$ we have

$$\left(\frac{e^\beta }{I}+\frac{3}{2}\right)p(1)=2,$$
(59)
so that if $`IAe^{\beta x^{}}`$, then $`p(1)2Ae^{\beta (1x^{})}`$. It is readily seen that $`p(x)2Ae^{\beta (xx^{})}`$ satisfies Eq. (58) for $`x^{}<x1`$. The same equation then leads to $`p(x)=1`$ as $`\beta \mathrm{}`$, for $`x^2<x<x^{}`$. We may then develop the full solution using Eqs. (55) to (57); normalization implies $`x^{}=1/2`$. The result is the function $`p(x)`$ plotted in Fig. 7, which is in good agreement with the density found via simulation of the random-neighbor model. This solution is discontinous at $`x^{}`$, $`x^2`$, $`x^4`$,…, and exhibits an integrable divergence at $`x=0`$.

The value of $`x^{}`$ can be confirmed through the reasoning developed in previous subsection for the peripheral square model, which in this case implies $`1=2(1x^{})`$ so that $`x^{}=1/2`$.

## IV Numerical simulation and Critical exponents

We now compare the mean-field theory predictions with simulation results. We estimate the probability density $`p(x)`$ on the basis of a histogram of barrier frequencies, dividing  into 100 subintervals. Histograms are accumulated after $`N_{st}`$ time steps, as required for the system (a ring of $`N`$ sites) to relax to the stationary state. In order to improve statistics, we average over $`N_r`$ realizations.

The simulation results for the original BS model are shown in Fig. 8. We observe qualitative agreement between nearest-neighbor (NN) and random-neighbor (RN) versions. (The simulation parameters are: $`N=1000`$, $`N_{st}=10^6`$, $`N_r=10^3`$ (NN); $`N=2000`$, $`N_{st}=10^5`$, $`N_r=10^3`$ (RN).) Here and in all other cases, the simulation result for the random neighbor model appears to converge (as expected) to the mean-field prediction. In light of the discussion in Sec. IIIB, it is reasonable to regard the rounding of the step function as a finite-size effect.

Fig. 4 presents similar results for the radical variant. We notice that the NN and RN versions exhibit qualitatively similar probability densities, differing mainly in the value of the threshold $`x^{}`$. (In this case we use $`N_{st}=10^6`$ and $`N_r=10^3`$, with system sizes $`N=1000`$ (NN) and 2000 (RN).)

Simulation results for the centered and peripheral square versions are shown in Fig. 9 and 10, resp.. For the centered square, the NN and RN, probability densities are quite similar, differing mainly in the value of the threshold and in the inclination of the central portion (Fig. 9). (The simulation parameters are $`N=2000`$, $`N_{st}=10^7`$, $`N_r=500`$ (NN), $`N=10000`$, $`N_{st}=10^6`$, $`N_r=10^3`$ (RN).) On the other hand, for the peripheral square version, the nearest-neighbor and random-neighbor densities are somewhat different, since the latter exhibits various steps, while only one such step is evident in the former (see Fig. 10). Further study is needed to determine whether this is represents a qualitative difference between the two formulations, or is instead due to finite size and/or finite numerical resolution. (Note that for the square models we use larger lattices, which were necessary to observe clear singularities. The simulation parameters are: $`N=1000`$, $`N_{st}=10^6`$, $`N_r=10^3`$ (NN), $`N=50000`$, $`N_{st}=10^7`$, $`N_r=20`$ (RN).)

Several quantities are known to display power-law behaviour in the Bak-Sneppen model . In particular, we studied the distribution $`P_J(r)`$ that sucessive updated sites are separated by a distance $`r`$. In the original model, $`P_J(r)r^\pi `$, with $`\pi =3.23(2)`$ . (Figures in parentheses denote uncertainties.) We performed simulations of the modified models using $`10^9`$ time steps on lattices of 2000 or more sites, yielding (see Fig. 11), $`\pi =3.22(2)`$ for the original BS model; $`\pi =3.27(2)`$ for the radical variant; $`\pi =3.25(2)`$ and $`3.24(2)`$ for the centered and peripheral square variants, respectively. Thus we find strong evidence that all of the variants introduced here belong to the same universality class as the original model.

## V Conclusions

In order to understand the implications of extremal dynamics, we propose several modified Bak-Sneppen models. Although different updating rules lead to completely different probability densities, they are always singular at one or more points. The step-like singularities appear in the limiting densities, either as $`\beta \mathrm{}`$ in the finite temperature model (on an infinite lattice), or as the system size $`N\mathrm{}`$, under extremal dynamics. Thus the double limit of infinite size and zero temperature is required for the BS model or it variants to generate a singular probability density.

A remarkable feature common to the original model and all the variants considered, is that in the infinite-size limit, and under extremal dynamics, a certain interval $`D[0,1]`$ has probability zero, and yet contains the extremal site with probability one. If we regard the density of active sites, $`\text{Prob}[xD]`$, as the order parameter $`\rho `$, then the BS model and its variants are seen to realize the ‘SOC limit’  $`\rho 0^+`$. Correlations between site variables $`x_i,x_j,\mathrm{},x_n`$ are zero unless one or more of these values falls in $`D`$.

Our conclusions regarding the form of the stationary densities are based on mean-field analyses that are exact for the random-neighbor versions, as is verified numerically. We also present a pair approximation for the original model. An important point is that mean-field theory captures the form of the probability density correctly, as shown via simulations of the modified models on a one-dimensional lattice. The latter suggest that the critical exponents are independent of dynamics.

The Bak-Sneppen model appears to be a prototype for a large universality class, since countless variants, beyond those presented here, are possible. We expect any dynamics that respects the symmetry of the original model (that is, spatial isotropy), and that does not introduce new conserved quantities, to have the same exponents as the original model. This is interesting, since the same critical behavior may subsist, as it were, upon stationary distributions of very different forms. Aside from its intrinsic interest, the question of universality is important for applications, since the precise form of the dynamics in a specific setting (e.g., evolution) is generally unknown, and probably quite different from that of the original model. Power laws and a singular stationary density only appear in the extremal dynamics limit, which may be difficult to realize in spatially extended natural systems.

ACKNOWLEDGMENTS

We thank CNPq and CAPES, Brazil, for finacial suport.


format: text

text:
# Search for quantum transducers between electromagnetic and gravitational radiation: A measurement of an upper limit on the transducer conversion efficiency of yttrium barium copper oxide

## I Acknowledgments

We thank Dung-Hai Lee, Jon-Magne Leinaas, Robert Littlejohn, Joel Moore, Richard Packard, Paul Richards, Daniel Solli, Achilles Speliotopoulos, Neal Snyderman, and Sandy Weinreb for stimulating discussions. This work was supported in part by the ONR.

## II Figure Captions

Figure 1: Quantum transducer between electromagnetic (EM) and gravitational (GR) radiation, consisting of a quantum fluid, such as the quantum Hall fluid, which possesses charge and spin. The minimal-coupling rule for an electron coupled to curved spacetime via its charge and spin, results in two processes. In (a) an EM plane wave is converted upon reflection from the quantum fluid into a GR plane wave; in (b), which is the reciprocal or time-reversed process, a GR plane wave is converted upon reflection from the quantum fluid into an EM plane wave.

Figure 2: Simplified schematic of a Hertz-type experiment, in which gravitational radiation at 12 GHz could be emitted and received using two superconductors. The “Microwave Source” generated electromagnetic radiation at 12 GHz (“EM wave”), which impinged on Superconductor A, could be converted upon reflection into gravitational radiation (“GR wave”). The GR wave, but not the EM wave, could pass through the “Faraday Cages.” In the far field of Superconductor A, Superconductor B could reconvert upon reflection the GR wave back into an EM wave at 12 GHz, which could then be detected by the “Microwave Detector.”

Figure 3: The T-antenna (expanded view on the left) used as antennas inside the “Source Can” and the “Receiver Can.” The YBCO samples were oriented so that a GR microwave beam could be directed from one YBCO sample to the other along a straight line of sight.

Figure 4: Data from the Hertz-type gravity-wave experiment using YBCO superconductors as transducers between EM and GR radiation. In (a), the cover lids were off both the source and the receiver cans, so that a small leakage signal (the central spike) could serve to test the system. In (b), both cover lids were on the cans, but no detectable signal of coupling between the cans could be seen above the noise. Both YBCO samples were immersed in liquid nitrogen for these data.


text:
# Density perturbations in a brane-world universe with dark radiation

## I Introduction

In brane-world cosmology, which has emerged in the context of recent developments in M theory, the observable universe is a 1+3-dimensional “brane” surface embedded in a 1+3+$`d`$-dimensional “bulk” spacetime. Fields and particles in the non-gravitational sector are confined to the brane, while gravity propagates in the bulk. Simple phenomenological brane-world models are the cosmological generalizations  of the Randall-Sundrum 2 model , which has a self-gravitating brane with Minkowski geometry embedded in an infinite part of 5-dimensional Anti de Sitter space (AdS5). Gravity on the brane is prevented from “leaking” into the infinite extra dimension at low energies via “warping” of the metric by the negative bulk cosmological constant, $`\mathrm{\Lambda }_5=6/\mathrm{}^2,`$ where $`\mathrm{}`$ is the curvature radius of AdS5. Modes of the 5D graviton have an effective mass on the brane, since the projection onto the brane of a null 5D graviton momentum vector is in general a 4D timelike momentum vector. The massless mode arises when the projection is null, and corresponds to the 4D graviton, which is dominant at low energies. The massive Kaluza-Klein (KK) modes produce corrections to the gravitational potential in the weak-field static limit, which are $`O(\mathrm{}^2/r^2)`$ . Table-top experiments currently impose an upper bound $`\mathrm{}0.1`$ mm.

On the brane, the negative $`\mathrm{\Lambda }_5`$ is offset by the positive brane tension $`\lambda `$, and the effective cosmological constant on the brane is

$$\mathrm{\Lambda }=\frac{1}{2}(\mathrm{\Lambda }_5+\kappa ^2\lambda ),$$
(1)
where $`\kappa ^2=8\pi /M_4^2`$ and $`M_410^{19}`$GeV is the effective Planck scale on the brane. Fine-tuning can set $`\mathrm{\Lambda }=0`$. Because of the large extra dimension, the true fundamental gravity scale can be as low as $``$ TeV in some brane-world scenarios, but in generalized RS2 models it is higher:

$$M_5^3=\frac{M_4^2}{\mathrm{}},\mathrm{}<0.1\mathrm{mm}M_5>10^5\mathrm{TeV},\lambda >(100\mathrm{GeV})^4.$$
(2)
At high energies ($`\rho \lambda `$) in the early universe, gravity becomes 5-dimensional and there are significant corrections to standard cosmological dynamics.

The unperturbed cosmological brane-world is a Friedmann brane in a Schwarzschild-AdS5 bulk  (where the expansion of the universe can be interpreted as motion of the brane in the bulk), satisfying the usual energy conservation equation but modified Friedmann equations. For a spatially flat universe without cosmological constant, these are

$`H^2`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{3}}\rho \left(1+{\displaystyle \frac{\rho }{2\lambda }}\right)+{\displaystyle \frac{m}{a^4}},`$ (3)
$`\dot{H}`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{2}}\rho (1+w)\left(1+{\displaystyle \frac{\rho }{\lambda }}\right)2{\displaystyle \frac{m}{a^4}},`$ (4)
where $`w=p/\rho `$ and $`m`$ is a constant, proportional to the mass of the bulk black hole. To avoid a naked singularity in the bulk, we take $`m0`$. The tidal, Coulomb effect of the 5-dimensional black hole on the brane is an effective radiative ($`a^4`$) term, the so-called “dark radiation”. Nucleosynthesis is sensitive to the inclusion of additional relativistic energies not thermally coupled to the radiation plasma, and so is the cosmic microwave background. This places limits on the amount of dark radiation compatible with observations

$$\left(\frac{\rho ^{}}{\rho _{\mathrm{rad}}}\right)_{\mathrm{nuc}}0.03,$$
(5)
so that $`\rho ^{}\rho _{\mathrm{rad}}`$ always, and the dark radiation is negligible by the matter era. The quadratic correction terms $`\rho ^2/\lambda `$ are dominant at high energies, $`\rho \lambda `$, and negligible at low energies. Since $`\lambda 1\mathrm{MeV}`$, the quadratic effects are negligible by nucleosynthesis. At late times, the standard 4-dimensional evolution is therefore recovered.

In the absence of dark radiation, the exact solution (with $`\dot{w}=0=m`$) is

$$a=\text{const}[t(t+t_\lambda )]^{1/3(w+1)},t_\lambda =\frac{M_4}{\sqrt{\pi \lambda }}<10^9\mathrm{sec}.$$
(6)
In the radiation era, $`w=\frac{1}{3}`$, the exact solution with dark radiation ($`m0`$) is

$$a=\text{const}[t(t+t_\lambda )]^{1/4},t_\lambda =\frac{\sqrt{3}M_4}{4\sqrt{\pi \lambda }(1+3m/\kappa ^2\rho a^4)}.$$
(7)
At late times, $`tt_\lambda `$, we recover the standard solution $`at^{2/3(w+1)}`$, whereas at early times (high energies), $`tt_\lambda `$, the evolution is very different: $`at^{1/3(w+1)}`$.

The high-energy brane-world correction provides increased Hubble damping in inflation; instead of $`H\sqrt{V}`$ as in general relativity, the modified Friedmann equation (3) shows that $`HV`$ when $`V\lambda `$. This means in particular that slow-roll inflation is possible even for potentials $`V(\phi )`$ that would be too steep in standard cosmology . In addition, large-scale scalar perturbations generated by high-energy slow-roll inflation have an enhanced amplitude $`A_\mathrm{s}`$ compared with the standard general relativity amplitude $`(A_\mathrm{s})_{\mathrm{gr}}`$ :

$$A_\mathrm{s}^2\left(A_\mathrm{s}^2\right)_{\mathrm{gr}}\left(\frac{V}{\lambda }\right)^2.$$
(8)
The background dynamics of brane-world cosmology are simple because they are effectively 4-dimensional. Cosmological perturbations introduce truly 5-dimensional degrees of freedom, so that the 5-dimensional bulk perturbation equations must be solved in order to solve for perturbations on the brane. These 5-dimensional equations are partial differential equations for the 3-dimensional Fourier modes, subject to complicated initial and boundary conditions. However, on large scales, if one can neglect gradient terms, the density or curvature perturbations can be determined without knowing the 5-dimensional solutions . The qualitative 5-dimensional effects on large-scale density perturbations during inflation and the transition to radiation-domination were investigated via a toy model in . Here we generalise the results of  to include the dark radiation term $`m/a^4`$ in the background. The presence of this term in the background can have a significant effect on the perturbations.

## II Field equations

The effective field equations on the brane are derived from the 5D field equations in the bulk, $`{}_{}{}^{(5)}G_{AB}^{}=\mathrm{\Lambda }_5^{(5)}g_{AB}`$, by projecting the 5D curvature (using the Gauss-Codazzi equations), and then imposing the Darmois-Israel junction conditions at the brane (with $`Z_2`$-symmetry) :

$$G_{ab}=\mathrm{\Lambda }g_{ab}+\kappa ^2T_{ab}+6\frac{\kappa ^2}{\lambda }𝒮_{ab}_{ab}.$$
(9)
Here $`𝒮_{ab}(T_{ab})^2`$ is the high-energy correction term, which is negligible for $`\rho \lambda `$, while $`_{ab}`$ is the projection of the bulk Weyl tensor, encoding corrections from 5D graviton (KK) effects (and giving the dark radiation in the background). From the brane-observer viewpoint, the energy-momentum corrections in $`𝒮_{ab}`$ are local, whereas the KK corrections in $`_{ab}`$ are nonlocal , since they incorporate 5D gravity wave modes. These nonlocal corrections cannot be determined purely from data on the brane, and so the effective field equations are not a closed system. One needs to supplement them by 5D equations governing $`_{ab}`$, which are obtained from the 5D Einstein and Bianchi equations .

The trace free $`_{ab}`$ contributes an effective KK energy density $`\rho ^{}`$, pressure $`\rho ^{}/3`$, momentum density $`q_a^{}`$ and anisotropic stress $`\pi _{ab}^{}`$ on the brane,

$$\frac{1}{\kappa ^2}_{ab}=\rho ^{}\left(u_au_b+\frac{1}{3}h_{ab}\right)+q_a^{}u_b+q_b^{}u_a+\pi _{ab}^{},$$
(10)
where $`u^a`$ is a physically determined 4-velocity on the brane and $`h_{ab}=g_{ab}+u_au_b`$ projects into the comoving rest space at each event. The brane “feels” the bulk gravitational field through these terms. In the background, $`q_a^{}=0=\pi _{ab}^{}`$, since only the dark radiation term $`\rho ^{}`$ ($`ma^4`$) is compatible with Friedmann symmetry. The KK momentum density defines a velocity $`v_a`$ of the Weyl “fluid” relative to $`u^a`$, by $`q_a^{}=\rho ^{}v_a`$.

For a perfect fluid or scalar field, we choose $`u^a`$ as the frame in which there is no energy flux. Then the effective total energy density and pressure are

$`\rho ^{\text{eff}}`$ $`=`$ $`\rho \left(1+{\displaystyle \frac{\rho }{2\lambda }}+{\displaystyle \frac{\rho ^{}}{\rho }}\right),`$ (11)
$`p^{\text{eff}}`$ $`=`$ $`\rho \left[w+(1+2w){\displaystyle \frac{\rho }{2\lambda }}+{\displaystyle \frac{\rho ^{}}{3\rho }}\right].`$ (12)
Energy-momentum conservation, $`^bT_{ab}=0`$, holds on the brane. Together with the 4D Bianchi identity, this implies that $`^a_{ab}=6\kappa ^2^a𝒮_{ab}/\lambda `$, which shows qualitatively how 1+3 spacetime variations in the matter-radiation on the brane can source KK modes. These KK “conservation equations” in linearized form are

$`\dot{\rho }^{}+{\displaystyle \frac{4}{3}}\mathrm{\Theta }\rho ^{}+\text{D}^aq_a^{}=0,`$ (13)
$`\dot{q}_a^{}+4Hq_a^{}+{\displaystyle \frac{1}{3}}\text{D}_a\rho ^{}+{\displaystyle \frac{4}{3}}\rho ^{}A_a+\text{D}^b\pi _{ab}^{}={\displaystyle \frac{(\rho +p)}{\lambda }}\text{D}_a\rho ,`$ (14)
where $`\mathrm{\Theta }`$ is the volume expansion rate ($`=3H`$ in the background), $`A_a`$ is the 4-acceleration, and $`\text{D}_a`$ is the covariant spatial derivative in the comoving rest space. Spatial inhomogeneity ($`\text{D}_a\rho 0`$) is a source for KK modes. Qualitatively and geometrically this can be understood as follows : the non-uniform 5D gravitational field generated by inhomogeneous 4D matter-radiation contributes to the 5D Weyl tensor, which nonlocally “backreacts” on the brane via its projection $`_{ab}`$. Equation (14) shows that the source term is suppressed at low energies, and during quasi-de Sitter inflation on the brane.

Equations (13) and (14) are propagation equations for $`\rho ^{}`$ and $`q_a^{}`$. There is no propagation equation on the brane for $`\pi _{ab}^{}`$, so that one cannot determine the KK modes purely from data on the brane. In the 1+3-covariant description of braneworld perturbations , the KK anisotropic stress $`\pi _{ab}^{}`$ is isolated as the term that must be determined from 5D equations. Once $`\pi _{ab}^{}`$ is determined in this way, the 1+3 perturbation equations on the brane form a closed system. The KK terms act as source terms modifying the standard general relativity perturbation equations, together with the local high-energy corrections.

## III  Density perturbation equations

We define density and expansion (velocity) perturbation scalars, as in general relativity ,

$$\mathrm{\Delta }=\frac{a^2}{\rho }\text{D}^2\rho ,Z=a^2\text{D}^2\mathrm{\Theta },$$
(15)
and then define dimensionless KK perturbation scalars ,

$$U=\frac{a^2}{\rho }\text{D}^2\rho ^{},Q=\frac{a}{\rho }\text{D}^2q^{},\mathrm{\Pi }=\frac{1}{\rho }\text{D}^2\pi ^{},$$
(16)
where the scalar potentials $`q^{}`$ and $`\pi ^{}`$ are defined by $`q_a^{}=\text{D}_aq^{}`$, $`\pi _{ab}^{}=(\text{D}_a\text{D}_b\frac{1}{3}h_{ab}\text{D}^2)\pi ^{}`$. The dark radiation fluctuation $`U`$ is present even if there is no dark radiation in the background ($`\rho ^{}=0`$). It leads to an isocurvature (non-adiabatic) mode, even when the matter perturbations are assumed adiabatic . We define the total effective dimensionless entropy $`S^{\mathrm{eff}}`$ via the non-adiabatic part of the effective pressure as

$$p^{\mathrm{eff}}S^{\mathrm{eff}}=a^2\text{D}^2p^{\mathrm{eff}}c_{\mathrm{eff}}^2a^2\text{D}^2\rho ^{\mathrm{eff}},$$
(17)
where $`c_{\mathrm{eff}}^2=\dot{p}^{\mathrm{eff}}/\dot{\rho }^{\mathrm{eff}}`$. Then

$`S^{\mathrm{eff}}={\displaystyle \frac{\left[3c_\mathrm{s}^21+\left(2+3w+3c_\mathrm{s}^2\right)\rho /\lambda \right]}{[(1+w)(1+\rho /\lambda )+4\rho ^{}/3\rho ][3w+3(1+2w)\rho /2\lambda +\rho ^{}/\rho ]}}\left[{\displaystyle \frac{4}{3}}{\displaystyle \frac{\rho ^{}}{\rho }}\mathrm{\Delta }(1+w)U\right],`$ (18)
where $`c_\mathrm{s}^2=\dot{p}/\dot{\rho }`$. If $`\rho ^{}=0`$ in the background, then $`U`$ is an isocurvature mode: $`S^{\mathrm{eff}}U`$. If $`\rho ^{}0`$ in the background, then the weighted difference between $`U`$ and $`\mathrm{\Delta }`$ determines the isocurvature mode: $`S^{\mathrm{eff}}(4\rho ^{}/3\rho )\mathrm{\Delta }(1+w)U`$. At very high energies, $`\rho \lambda `$, the entropy is suppressed by the factor $`\lambda /\rho `$. At low energies in the radiation era ($`c_\mathrm{s}^2=\frac{1}{3}`$) it is suppressed by the factor $`\rho /\lambda `$.

The covariant density perturbation equations on the brane  reduce to

$`\dot{\mathrm{\Delta }}`$ $`=`$ $`3wH\mathrm{\Delta }(1+w)Z,`$ (19)
$`\dot{Z}`$ $`=`$ $`2HZ\left({\displaystyle \frac{c_\mathrm{s}^2}{1+w}}\right)\text{D}^2\mathrm{\Delta }\kappa ^2\rho U{\displaystyle \frac{1}{2}}\kappa ^2\rho \left[1+(4+3w){\displaystyle \frac{\rho }{\lambda }}\left({\displaystyle \frac{4c_\mathrm{s}^2}{1+w}}\right){\displaystyle \frac{\rho ^{}}{\rho }}\right]\mathrm{\Delta },`$ (20)
$`\dot{U}`$ $`=`$ $`(3w1)HU+\left({\displaystyle \frac{4c_\mathrm{s}^2}{1+w}}\right)\left({\displaystyle \frac{\rho ^{}}{\rho }}\right)H\mathrm{\Delta }\left({\displaystyle \frac{4\rho ^{}}{3\rho }}\right)Za\text{D}^2Q,`$ (21)
$`\dot{Q}`$ $`=`$ $`(3w1)HQ{\displaystyle \frac{1}{3a}}U{\displaystyle \frac{2}{3}}a\text{D}^2\mathrm{\Pi }+{\displaystyle \frac{1}{3a}}\left[\left({\displaystyle \frac{4c_\mathrm{s}^2}{1+w}}\right){\displaystyle \frac{\rho ^{}}{\rho }}3(1+w){\displaystyle \frac{\rho }{\lambda }}\right]\mathrm{\Delta },`$ (22)
where $`\dot{f}u^a_af`$, and we have corrected the minor errors in the equations given in . The KK anisotropic stress term $`\mathrm{\Pi }`$ occurs only via its Laplacian, $`\text{D}^2\mathrm{\Pi }`$. If we can neglect this term on large scales, then the system of density perturbation equations closes on super-Hubble scales . An equivalent statement applies to the large-scale curvature perturbations . We note that in the low-energy approximation scheme used by Koyama ,

$$\mathrm{\Pi }\text{D}^2U,$$
(23)
so that in particular $`\mathrm{\Pi }0`$ on large scales.

On large scales in the braneworld, neglecting gradient terms, density perturbations are then described by the set $`\{\mathrm{\Delta },Z,U,Q\}`$, where $`Q`$ decouples from the other variables. KK effects introduce two new isocurvature modes on large scales, associated with $`U`$ and $`Q`$ .

We can find a first integral of the large-scale system by following the same approach used in general relativity . A covariant local curvature perturbation is defined by the gradient of $`R^{}`$, the 3-Ricci scalar for $`u^a`$-observers. This scalar is given by the Codazzi equation on the brane ,

$$R^{}=\frac{2}{3}\mathrm{\Theta }^2+2\kappa ^2\rho ^{\mathrm{eff}}.$$
(24)
Then the dimensionless curvature perturbation scalar is defined by $`C=a^4\text{D}^2R^{}`$. This gives

$$C=4a^2HZ+2\kappa ^2a^2\rho \left(1+\frac{\rho }{\lambda }\right)\mathrm{\Delta }+2\kappa ^2a^2\rho U.$$
(25)
Using the large-scale equations (19)–(21), we find that $`C`$ is locally conserved (i.e., constant along fundamental worldlines) on large scales,

$$C=C_0,\dot{C}_0=0.$$
(26)
(The minor errors in  led to the erroneous conclusion that $`C`$ was not in general locally conserved on the brane.)

Thus $`C`$ leads to a first integral for the large-scale system. We also replace $`\mathrm{\Delta }`$ by the dimensionless perturbation variable

$$\mathrm{\Phi }=\kappa ^2a^2\rho \mathrm{\Delta }.$$
(27)
In general relativity, $`\mathrm{\Phi }`$ is the covariant analogue of the Bardeen metric perturbation variable $`\mathrm{\Phi }_H`$ . In the brane-world, high-energy and KK effects mean that the analogue of $`\mathrm{\Phi }_H`$ is a complicated generalisation of Eq. (27) (see ). The key problem in trying to use the generalised variable is that it contains $`\mathrm{\Pi }`$, which remains undetermined. (This is the same reason that the large-scale Sachs-Wolfe effect on the brane cannot be determined from the density or curvature perturbations .) In any case, $`\mathrm{\Phi }`$ as defined above leads to useful simplification of the perturbation equations.

The coupled system of Eqs. (19)–(21), which is closed on large scales (neglecting gradient terms), reduces via Eqs. (26) and (27) to the system

$`\dot{\mathrm{\Phi }}`$ $`=`$ $`H\left[1+(1+w){\displaystyle \frac{\kappa ^2\rho }{2H^2}}\left(1+{\displaystyle \frac{\rho }{\lambda }}\right)\right]\mathrm{\Phi }\left[(1+w){\displaystyle \frac{a^2\kappa ^4\rho ^2}{2H}}\right]U+\left[(1+w){\displaystyle \frac{\kappa ^2\rho }{4H}}\right]C_0,`$ (28)
$`\dot{U}`$ $`=`$ $`H\left[13w+{\displaystyle \frac{2\kappa ^2\rho ^{}}{3H^2}}\right]U{\displaystyle \frac{2\rho ^{}}{3a^2H\rho }}\left[1+{\displaystyle \frac{\rho }{\lambda }}{\displaystyle \frac{6c_\mathrm{s}^2H^2}{(1+w)\kappa ^2\rho }}\right]\mathrm{\Phi }+\left[{\displaystyle \frac{\rho ^{}}{3a^2H\rho }}\right]C_0.`$ (29)
This closed system for $`\mathrm{\Phi },U`$ generalises the equations given in  to the case $`\rho ^{}0`$. Once $`\mathrm{\Phi }`$ and $`U`$ are solved for, $`Q`$ is directly determined by Eq. (22) with $`\text{D}^2\mathrm{\Pi }`$ set to zero.

## IV Solutions

The closed system for $`\mathrm{\Phi }`$ and $`U`$ may be rewritten on using the Friedmann equation (3) and changing the time variable to the number of e-folds, $`tN=\mathrm{ln}(a/a_0)`$. This yields

$`\mathrm{\Phi }^{}`$ $`=`$ $`\left[1+{\displaystyle \frac{3(1+w)(1+\rho /\lambda )}{2(1+\rho /2\lambda +\rho ^{}/\rho )}}\right]\mathrm{\Phi }\left[{\displaystyle \frac{3(1+w)}{2(1+\rho /2\lambda +\rho ^{}/\rho )}}a^2\kappa ^2\rho \right]U+\left[{\displaystyle \frac{3(1+w)}{4(1+\rho /2\lambda +\rho ^{}/\rho )}}\right]C_0,`$ (30)
$`U^{}`$ $`=`$ $`\left[13w+{\displaystyle \frac{2\rho ^{}/\rho }{(1+\rho /2\lambda +\rho ^{}/\rho )}}\right]U\left[{\displaystyle \frac{\rho ^{}}{\rho }}\left\{{\displaystyle \frac{(1+\rho /\lambda )}{(1+\rho /2\lambda +\rho ^{}/\rho )}}{\displaystyle \frac{2c_\mathrm{s}^2}{(1+w)}}\right\}{\displaystyle \frac{2}{a^2\kappa ^2\rho }}\right]\mathrm{\Phi }`$ (32)
$`+\left[{\displaystyle \frac{\rho ^{}}{\rho }}{\displaystyle \frac{1}{(1+\rho /2\lambda +\rho ^{}/\rho )}}{\displaystyle \frac{1}{a^2\kappa ^2\rho }}\right]C_0,`$
where $`f^{}=df/dN`$ and

$$a^2\kappa ^2\rho =\beta \mathrm{exp}(1+3w)𝑑N,\beta =a_0^2\kappa ^2\rho _0.$$
(33)
Here $`\beta `$ is dimensionless and can be given an arbitrary positive value by suitable choice of $`a_0`$.

For comparison, in the general relativity case ($`\rho /\lambda =0=\rho ^{}=U`$) these two equations reduce to the single equation

$$\mathrm{\Phi }_{\mathrm{gr}}^{}=\frac{1}{2}(5+3w)\mathrm{\Phi }_{\mathrm{gr}}+\frac{3}{4}(1+w)C_0.$$
(34)
If $`w`$ is constant, then the non-decaying attractor solution is

$$\mathrm{\Phi }_{\mathrm{gr}}=\frac{3(1+w)}{2(5+3w)}C_0.$$
(35)
Deviations from general relativity arise from early-universe high-energy effects (when $`\rho \lambda `$), and from the KK effects encoded in $`\rho ^{}`$ and $`U`$. At late times, the high-energy corrections become negligible. The background dark radiation decays like ordinary radiation ($`\rho ^{}a^4\rho _{\mathrm{rad}}`$), and is negligible at late times by Eq. (5). However, it can leave a signature that survives at late times, as we show below.

### A No dark radiation in the background

Equation (32) is much simpler in the case when there is no dark radiation in the background ($`\rho ^{}=0`$), i.e. when there is no bulk black hole. (Note also that the errors in  do not affect this case, so that the solutions given there are correct.) When $`\rho ^{}=0`$ we get

$$U=U_0\mathrm{exp}(13w)𝑑N.$$
(36)
Since $`\rho ^{}=0`$, it follows from Eq. (18) that $`U`$ is an isocurvature mode. This isocurvature mode associated with dark radiation perturbations therefore decays for all equations of state less stiff than radiation, i.e. for $`w<\frac{1}{3}`$.

If $`w`$ is constant, then Eq. (30) becomes

$$\mathrm{\Phi }^{}=\left[\frac{(5+3w)+(4+3w)\rho /\lambda }{2+\rho /\lambda }\right]\mathrm{\Phi }\left[\frac{3(1+w)}{2(1+\rho /2\lambda )}\right]e^{2N}\stackrel{~}{U}_0+\left[\frac{3(1+w)}{4(1+\rho /2\lambda )}\right]C_0,$$
(37)
where $`\stackrel{~}{U}_0=\beta U_0`$. The $`U`$ contribution to $`\mathrm{\Phi }`$ is suppressed at low energies, even if $`w>\frac{1}{3}`$ (when $`U`$ grows). Thus the general relativity solution Eq. (35) for constant $`w`$ is a low-energy attractor for $`\rho ^{}=0`$:

$$\mathrm{\Phi }_{\mathrm{low}}\frac{3(1+w)}{2(5+3w)}C_0.$$
(38)
Thus there is no surviving imprint in $`\mathrm{\Phi }`$ of the dark radiation perturbations (the isocurvature mode) when $`\rho ^{}=0`$ in the background.

At high energies/ early times, Eq. (37) becomes

$$\mathrm{\Phi }^{}(4+3w)\mathrm{\Phi }\frac{\lambda }{\rho _0}\left[3(1+w)\right]e^{(1+3w)N}\stackrel{~}{U}_0+\frac{\lambda }{\rho _0}\left[\frac{3}{2}(1+w)\right]e^{3(1+w)N}C_0,$$
(39)
where $`\lambda \rho _0`$. Thus in the early universe, the contribution of the entropy mode $`U`$ is sensitive to the equation of state. This contribution decays for $`w<\frac{1}{3}`$, including the case of inflation. For $`w\frac{1}{3}`$, Eq. (39) shows that the non-decaying solution is

$`\mathrm{\Phi }_{\mathrm{high}}`$ $``$ $`{\displaystyle \frac{\lambda }{\rho _0}}\left[{\displaystyle \frac{3(1+w)}{2(7+6w)}}C_0\right]e^{3(1+w)N}{\displaystyle \frac{\lambda }{\rho _0}}\left[{\displaystyle \frac{3\beta (1+w)}{(5+6w)}}U_0\right]e^{(1+3w)N},`$ (40)
so that in the very high energy limit,

$$\mathrm{\Phi }_{\mathrm{high}}\frac{3}{2}\frac{\lambda }{\rho _0}(1+w)\left[\frac{C_0}{7+6w}\frac{2\stackrel{~}{U}_0}{5+6w}\right].$$
(41)
Hence, for constant $`w`$ ($`\frac{1}{3}`$), high-energy effects initially suppress $`\mathrm{\Phi }`$ (by the factor $`\lambda /\rho _0`$), before it starts to grow. This growth leads to the late-time general relativity value, Eq. (38). By contrast, in general relativity, $`\mathrm{\Phi }`$ remains constant when $`w`$ is constant.

The most relevant case is $`w=\frac{1}{3}`$. If inflation takes place on the brane, reheating will initiate a radiation era. Provided the reheat temperature is high enough, the radiation era will start at high energies. Large-scale perturbations $`\mathrm{\Phi }`$ will be amplified according to Eqs. (41) and (38) with $`w=\frac{1}{3}`$.

For slow-roll inflation, $`w`$ is not constant, but close to -1 and slowly increasing, i.e., $`1+wϵ`$, with $`0<ϵ1`$ and $`H^1|\dot{ϵ}|=|ϵ^{}|1`$. Then Eq. (39) leads to

$$\mathrm{\Phi }_{\mathrm{high}}\frac{3}{2}ϵ\frac{\lambda }{\rho _0}C_0e^{3ϵN},$$
(42)
so that $`\mathrm{\Phi }`$ has a slowly growing mode during high-energy slow-roll inflation. This is different from general relativity, where $`\mathrm{\Phi }`$ is constant during slow-roll inflation. Thus more amplification of $`\mathrm{\Phi }`$ can be achieved than in general relativity, consistent with Eq. (8). The evolution of $`\mathrm{\Phi }`$ is illustrated for a toy model of inflation-to-radiation in Fig. 1 (from ). The early (growing) and late time (constant) attractor solutions are seen explicitly in the plots.

### B Nonzero dark radiation in the background

The presence of dark radiation in the background has interesting qualitative implications, even though the quantitative constraints from observations on the amount of dark radiation are quite severe. The general relativity solution, Eq. (35) is no longer the late-time attractor when $`\rho ^{}0`$.

The dark radiation factor, for constant $`w`$, behaves as

$$\frac{\rho ^{}}{\rho }=\alpha e^{(3w1)N},\alpha =\frac{\rho _0^{}}{\rho _0},$$
(43)
where $`\alpha `$ is a dimensionless constant. For equations of state less stiff than radiation ($`w<\frac{1}{3}`$), the dark radiation factor redshifts away. If $`w>\frac{1}{3}`$, then this factor grows with expansion. In the radiation era, $`w=\frac{1}{3}`$, the dark radiation factor is constant.

If there was a period of inflation on the brane before the radiation era, then dark radiation would be redshifted to negligible levels. However, the reheating era that creates the radiation, could at the same time create dark radiation via 5D graviton emission in high-energy interactions . The fraction of dark radiation is limited by observational constraints, Eq. (5), so that

$$\alpha 0.03.$$
(44)
We can therefore perform a perturbative solution of Eqs. (30) and (32), up to first order in $`\alpha `$. The zero-order solutions are given in the previous subsection.

We find that the non-decaying solutions for radiation in the low-energy regime are

$`\mathrm{\Phi }_{\mathrm{low}}`$ $``$ $`{\displaystyle \frac{C_0}{3}}(1\alpha ),`$ (45)
$`\stackrel{~}{U}_{\mathrm{low}}`$ $``$ $`\stackrel{~}{U}_0+\alpha {\displaystyle \frac{C_0}{3}}e^{2N}.`$ (46)
It follows that the first-order correction to $`\stackrel{~}{U}`$ (where $`\stackrel{~}{U}=\stackrel{~}{U}^{(0)}+\stackrel{~}{U}^{(1)}`$) behaves like the large-scale density perturbations for radiation in general relativity:

$$\stackrel{~}{U}_{\mathrm{low}}^{(1)}\mathrm{\Delta }a^2.$$
(47)
In the high-energy regime, the non-decaying solutions are

$`\mathrm{\Phi }_{\mathrm{high}}`$ $``$ $`{\displaystyle \frac{\lambda }{\rho _0}}\left[{\displaystyle \frac{2}{9}}C_0e^{4N}{\displaystyle \frac{4}{7}}\stackrel{~}{U}_0e^{2N}\right]+4\alpha \left({\displaystyle \frac{\lambda }{\rho _0}}\right)^2\left[\left({\displaystyle \frac{2}{63}}C_0{\displaystyle \frac{4}{49}}\stackrel{~}{U}_0\right)e^{2N}+{\displaystyle \frac{4}{77}}\stackrel{~}{U}_0e^{6N}{\displaystyle \frac{2}{117}}C_0e^{8N}\right],`$ (48)
$`\stackrel{~}{U}_{\mathrm{high}}`$ $``$ $`\stackrel{~}{U}_0+\alpha {\displaystyle \frac{\lambda }{\rho _0}}\left[{\displaystyle \frac{2}{9}}C_0\left(e^{6N}1\right){\displaystyle \frac{4}{7}}\stackrel{~}{U}_0\left(e^{4N}1\right)\right].`$ (49)
In the very high energy limit,

$`\mathrm{\Phi }_{\mathrm{high}}`$ $``$ $`{\displaystyle \frac{\lambda }{\rho _0}}\left[{\displaystyle \frac{2}{9}}C_0{\displaystyle \frac{4}{7}}\stackrel{~}{U}_0\right]+16\alpha \left({\displaystyle \frac{\lambda }{\rho _0}}\right)^2\left[{\displaystyle \frac{C_0}{273}}{\displaystyle \frac{4\stackrel{~}{U}_0}{539}}\right],`$ (50)
$`\stackrel{~}{U}_{\mathrm{high}}`$ $``$ $`\stackrel{~}{U}_0.`$ (51)
Comparing Eq. (50) with Eq. (41) shows that the high energy limit with dark radiation differs only negligibly from the no-dark-radiation limit.

This analysis shows that $`\mathrm{\Phi }`$ is initially suppressed, then begins to grow, as in the no-dark-radiation case, eventually reaching an attractor which is less than the no-dark-radiation attractor. The role of dark radiation is to suppress, by the small fractional amount $`\alpha `$, the final value of $`\mathrm{\Phi }`$, given in Eq. (45). The dark radiation perturbation $`U`$ starts from a constant value at very high energies, and then grows at late times like the density perturbations in radiation. Figures 2 and 3 confirm the qualitative analysis.

The effective total entropy, Eq. (18), has the following behaviour, to $`O(\alpha )`$, in the radiation era:

$`S_{\mathrm{high}}^{\mathrm{eff}}{\displaystyle \frac{8}{5}}{\displaystyle \frac{\lambda }{\rho _0}}\left[U_0\alpha {\displaystyle \frac{\lambda }{\rho _0}}\left({\displaystyle \frac{2C_0}{9\beta }}{\displaystyle \frac{4}{7}}U_0\right)\right],`$ (52)
$`S_{\mathrm{low}}^{\mathrm{eff}}4{\displaystyle \frac{\rho _0}{\lambda }}e^{4N}(12\alpha )U_0.`$ (53)
At high energies $`|S^{\mathrm{eff}}|`$ is suppressed by the factor $`\lambda /\rho _0`$, and the $`\alpha `$-correction is negligible. At low energies, the large factor $`\rho _0/\lambda `$ is overwhelmed by the redshift factor $`e^{4N}`$.

Without dark radiation in the background, $`U`$ in the radiation era remains constant. When there is dark radiation in the background, $`U`$ grows. The growth of $`U`$ leads to a decrease in the final value of $`\mathrm{\Phi }`$.

An extreme case of loss of power in $`\mathrm{\Phi }`$ occurs when the dark radiation redshifts more slowly than the matter, so that it eventually dominates over the matter. This happens when the matter equation of state is stiffer than radiation, $`w>\frac{1}{3}`$. For example, for stiff matter, $`w=1`$, we have

$$\frac{\rho ^{}}{\rho }a^2,$$
(54)
by Eq. (43). As a result, the decrease in $`\mathrm{\Phi }`$ becomes so pronounced that $`\mathrm{\Phi }`$ is driven to 0 at late times. By Eqs. (30) and (32), we find that

$$U_{\mathrm{low}}e^{4N},\mathrm{\Phi }_{\mathrm{low}}e^{2N},$$
(55)
at late times. The suppression of $`\mathrm{\Phi }`$ is confirmed by numerical integration, as shown in Fig. 4.

This example could be relevant to the steep inflation scenario on the braneworld. High-energy effects mean that steep potentials can drive inflation on the brane when they would not give inflation in general relativity. Typically, this inflation ends in a kinetic regime with $`w1`$, and with particle production being gravitational . If dark radiation can be produced by high-energy interactions, then its presence serves to suppress $`\mathrm{\Phi }`$ on large-scales, as in Fig. 4.

## V Conclusion

We have investigated the effects of dark radiation in the background on large-scale density perturbations (neglecting gradient terms) in a Randall-Sundrum 2 type brane-world. Dark radiation significantly complicates the perturbation equations. Since it is limited to be a small fraction of the radiation energy density, its effects can be qualitatively analysed via a perturbative approach. This approach, confirmed by our numerical simulations, shows that in the radiation era, the large-scale density perturbations are suppressed at late times by a small amount. At the same time, the large-scale perturbations in the dark radiation itself grow at late times. By contrast these perturbations have constant amplitude in a universe with no dark radiation in the background. The suppression of density perturbations becomes strong when the matter has a stiff (kinetic-dominated) equation of state.

Acknowledgments

We thank Peter Dunsby for pointing out the minor errors in , which are corrected in Eqs. (20)–(22). BG is supported by the British-Thai Scholarship Scheme and a Royal Thai Government Scholarship. CG and RM are supported by PPARC.


text:
# Hard-Sphere Fluids in Contact with Curved Substrates

## I Introduction

The properties of inhomogeneous fluids with spherical or cylindrical geometries have been investigated by many authors in recent years  because curvature plays an important role in many physical situations, including the following examples. In the context of nucleation phenomena  one considers the spontaneous formation of spherical droplets in a super-saturated vapor. The variation of the local curvature of confinements leads to a more complicated form of the depletion forces in colloidal suspensions . Wetting or critical adsorption of a fluids on curved substrates, such as spherical or cylindrical colloidal particles may lead to the formation of liquid bridges if two colloids come close to each other, and subsequently to flocculation . In studies of the hydrophobic effect one is concerned with the cost in free energy to immerse a hydrophobic particle into a liquid . Despite these efforts, the properties of curved fluid interfaces are not understood as fully as the planar one. Whether one deals with a drop of liquid surrounded by its vapor or with a big spherical particle immersed into a solvent of small particles, the statistical-mechanical description of such systems generates additional difficulties  which call for further studies .

Within this context, here we consider the case of a hard-sphere fluid in the grand canonical ensemble, characterized by the radius $`R`$ of the spheres and the bulk number density $`\rho `$ or equivalently the bulk packing fraction $`\eta =\left(4\pi /3\right)R^3\rho `$, in contact with a hard spherical or cylindrical wall. These curved walls can either be viewed as an external potential exerted on the hard-sphere fluid,

$$V_i^{ext}(r)=\{\begin{array}{cc}0,\hfill & r>R_i\hfill \\ \mathrm{},\hfill & r<R_i,\hfill \end{array}$$
(1)
where the subscript $`i=s`$, $`c`$ denotes the spherical and cylindrical symmetry, respectively, or as the surfaces of a second (big) species of radius $`R_i^{(b)}`$. The radii $`R_i`$ entering the external potential in Eq. (1) are related to the actual radii $`R_i^{(b)}`$ via $`R_i=R_i^{(b)}+R`$ (see Fig. 1).

The grand canonical potential of the system $`\mathrm{\Omega }_i`$ can be divided into a bulk part $`\mathrm{\Omega }_i^{bulk}=pV_i`$, where $`p`$ is the pressure of the hard-sphere fluid and $`V_i`$ is the volume accessible to the system, and an excess (surface) part $`\mathrm{\Omega }_i^{ex}`$: $`\mathrm{\Omega }_i=\mathrm{\Omega }_i^{bulk}+\mathrm{\Omega }_i^{ex}`$. The surface tension $`\gamma _i`$ is readily defined as the excess grand potential per unit area,

$$\gamma _i=\frac{\mathrm{\Omega }_i^{ex}}{A_i^{(b)}}=\frac{\mathrm{\Omega }_i\mathrm{\Omega }_i^{bulk}}{A_i^{(b)}},$$
(2)
where $`A_i^{(b)}`$ is the surface area of the curved wall. While $`\mathrm{\Omega }_i`$ is uniquely defined this division into bulk and surface excess parts is not unique. Depending on the choice of the dividing surface the surface tension can even change sign. As indicated by the superscript $`(b)`$ in Eq. (2) we use the actual wall corresponding to the radius $`R_i^{(b)}`$ as the dividing interface.

Within density functional theory (DFT), the grand potential $`\mathrm{\Omega }_i`$ of the system can be expressed as

$$\mathrm{\Omega }_i[\rho ]=[\rho ]+d^3r\rho (𝐫)(V_i^{ext}(𝐫)\mu ),$$
(3)
where $`[\rho ]=_{ex}[\rho ]+_{id}[\rho ]`$ is the intrinsic free energy functional that can be split into an ideal gas contribution $`_{id}[\rho ]`$ and an excess (over the ideal gas) term $`_{ex}[\rho ]`$. The only explicit dependence of $`\mathrm{\Omega }_i[\rho ]`$ on the external potential originates from the second term in Eq. (3). The equilibrium density profile $`\rho (𝐫)`$ satisfies the Euler-Lagrange equation

$$\frac{\delta \mathrm{\Omega }_i[\rho ]}{\delta \rho (𝐫)}=0.$$
(4)
For the equilibrium density profile $`\rho (𝐫)`$ the density functional reduces to the grand potential of the system, i.e., $`\mathrm{\Omega }_i=\mathrm{\Omega }_i[\rho (𝐫)]`$, and we can rewrite Eq. (2) as

$$\gamma _i=\frac{\mathrm{\Omega }_i[\rho (𝐫)]+pV_i}{A_i^{(b)}}.$$
(5)
From Eqs. (1), (3), and (4) it follows that the change in the grand potential due to an infinitesimal change of the radius $`R_i^{(b)}`$ of the curved wall, at constant chemical potential $`\mu `$ and temperature $`T`$, is given by

$$\left(\frac{\mathrm{\Omega }_i}{R_i^{(b)}}\right)_{\mu ,T}=d^3r\frac{\delta \mathrm{\Omega }_i[\rho ]}{\delta \rho (𝐫)}\frac{\rho (𝐫)}{R_i^{(b)}}+d^3r\rho (𝐫)\frac{V_i^{ext}(𝐫)}{R_i^{(b)}}.$$
(6)
The first term on the right-hand side of Eq. (6) vanishes by virtue of Eq. (4) while the second term, for the external potential given in Eq. (1), is equal to $`\beta ^1A_i\rho (R_i^+)`$, where $`\rho (R_i^+)`$ denotes the contact density of the hard-sphere fluid at the curved wall. This leads to the sum rules

$$\beta \left(\frac{\mathrm{\Omega }_i}{R_i^{(b)}}\right)_{\mu ,T}=\{\begin{array}{cc}4\pi R_s^2\rho (R_s^+),\hfill & \text{for a spherical wall}\hfill \\ 2\pi R_cL\rho (R_c^+),\hfill & \text{for a cylindrical wall.}\hfill \end{array}$$
(7)
In the case of a cylindrical wall we consider the thermodynamic limit $`\underset{L\mathrm{}}{lim}(\mathrm{\Omega }_c/L)`$ so that effects due to a finite length $`L`$ of the cylinder drop out. We note that the sum rules in Eq. (7) are valid for all one-component fluids in contact with hard spherical or cylindrical walls and are satisfied by all density functionals within weighted-density approximation (WDA) .

Using the actual radius $`R_i^{(b)}`$ also as the radius of the dividing interface (see Fig. 1), the grand potential $`\mathrm{\Omega }_i=pV_i+\gamma _iA_i^{(b)}`$ of the system is separated into a bulk and surface term with the accessible volume

$$V_i=V_i(\overline{R})V_i(R_i^{(b)}),$$
(8)
where $`V_s(r)=\frac{4\pi }{3}r^3`$, $`V_c(r)=\pi r^2L`$, and $`\overline{R}`$ is a macroscopicly large radius considered in the thermodynamical limit, and the surface areas

$$A_i^{(b)}=\{\begin{array}{cc}4\pi (R_s^{(b)})^2,\hfill & \text{for a spherical wall}\hfill \\ 2\pi R_c^{(b)}L,\hfill & \text{for a cylindrical wall}.\hfill \end{array}$$
(9)
The sum rule in Eq. (7) can then be expressed as

$$\rho (R_s^+)\left(\frac{R_s}{R_s^{(b)}}\right)^2=\beta p+\frac{2\beta \gamma _s(R_s^{(b)})}{R_s^{(b)}}+\beta \left(\frac{\gamma _s(R_s^{(b)})}{R_s^{(b)}}\right)_{\mu ,T}$$
(10)
for a spherical wall and

$$\rho (R_c^+)\frac{R_c}{R_c^{(b)}}=\beta p+\frac{\beta \gamma _c(R_c^{(b)})}{R_c^{(b)}}+\beta \left(\frac{\gamma _c(R_c^{(b)})}{R_c^{(b)}}\right)_{\mu ,T}$$
(11)
for a cylindrical wall. Equations (10) and (11) indicate that different versions of DFTs might lead to different contact values near curved walls even if the underlying equation of state is the same, which leads to identical contact values at a hard planar wall. As pointed out by J.R. Henderson , these sum rules have interesting implications. If $`\rho (R_i^+)`$ is analytic as a function of $`R_i^1=(R_i^{(b)}+R)^1`$ one has

$$\rho (R_i^+)=c_0+\frac{c_i^{(1)}}{R_i^{(b)}+R}+\frac{c_i^{(2)}}{(R_i^{(b)}+R)^2}+\frac{c_i^{(3)}}{(R_i^{(b)}+R)^3}+\mathrm{}$$
(12)
with $`c_0=\beta p`$ for both the spherical and the cylindrical wall and expansion coefficients $`c_i^{(j)}`$ for $`j=1,2,\mathrm{}`$. Note that the expansion in Eq. (12) is assumed to converge for arbitrary non-zero values of $`R_i`$. This assumption will be corroborated by our numerical results in Sec. III. With Eq. (12) the solution of the differential equation (10) can be written as

$$\beta \gamma _s(R_s^{(b)})=\frac{C_s^{(1)}}{2}\left[1+\frac{2C_s^{(2)}}{C_s^{(1)}}\frac{1}{R_s^{(b)}}+\frac{2c_s^{(3)}}{C_s^{(1)}}\frac{\mathrm{ln}(R_s^{(b)}/R)}{(R_s^{(b)})^2}+\frac{2D_s}{C_s^{(1)}}\frac{1}{(R_s^{(b)})^2}+𝒪\left(\frac{1}{(R_s^{(b)})^3}\right)\right],$$
(13)
where $`D_s`$ is an integration constant, $`C_s^{(1)}=c_s^{(1)}+2\beta pR`$ and $`C_s^{(2)}=c_s^{(2)}+c_s^{(1)}R+\beta pR^2`$. Likewise, with Eq. (12) the solution of the differential equation (11) can be written as

$$\beta \gamma _c(R_c^{(b)})=C_c^{(1)}\left[1+\frac{c_c^{(2)}}{C_c^{(1)}}\frac{\mathrm{ln}(R_c^{(b)}/R)}{R_c^{(b)}}+\frac{D_c}{C_c^{(1)}}\frac{1}{R_c^{(b)}}+𝒪\left(\frac{1}{(R_c^{(b)})^2}\right)\right],$$
(14)
where $`D_c`$ is an integration constant and $`C_c^{(1)}=c_c^{(1)}+\beta pR`$. In the case of the spherical wall $`\gamma _s(R_s^{(b)})`$ has the form

$$\beta \gamma _s(R_s^{(b)})=\beta \gamma _s(\mathrm{})\left[1\frac{2\delta _T}{R_s^{(b)}}+𝒪\left(\frac{\mathrm{ln}(R_s^{(b)}/R)}{(R_s^{(b)})^2}\right)\right]$$
(15)
where $`\gamma _s(\mathrm{})=\beta ^1\frac{c_s^{(1)}}{2}+pR`$ is the surface tension of the planar wall and

$$\delta _T=\frac{c_s^{(2)}+c_s^{(1)}R+\beta pR^2}{c_s^{(1)}+2\beta pR}$$
(16)
is the so-called Tolman length , which plays an important role in nucleation theory . We note that the value of the Tolman length depends on the choice of the dividing surface; here the dividing surface has the radius $`R_s^{(b)}`$. If these logarithmic terms in the expansion of the surface tensions do not vanish, the concept of the Tolman length apparently breaks down for cylindrically curved and generally shaped walls .

The analyticity of $`\gamma (R)`$ in terms of $`\frac{1}{R}`$ is of considerable importance because it provides the basis for the so-called Helfrich expansion of the surface free energy of arbitrarily curved surfaces in powers of the principal curvatures . Although this approach appears to be very useful for describing membranes  it has been shown recently that the presence of long-ranged dispersion forces destroys the aforementioned analyticity for fluid interfaces  and prevents the use of a Helfrich-type interface Hamiltonian for such systems. The possible occurrence of logarithmic terms in Eqs. (13) and (14) following from the simple expansion (12) raises the issue whether the absence of analyticity is not only due the presence of long-ranged forces but is an intrinsically geometrical effect which shows up even for surface free energies of hard-sphere fluids. The present study aims at clarifying this conceptual point.

To this end we investigate the curvature dependence of the surface tension as well as relationships between thermodynamic and local properties of the hard-sphere fluids following two different routes. First, in Sec. II the surface tension is determined from a bulk theory in which the curved wall is considered as the surface of an additional component at infinite dilution. Second, a numerical approach based on the minimization of the Rosenfeld fundamental measure theory (FMT) free energy functional  is presented in Sec. III. In this approach, the hard-sphere fluid is considered to be exposed to the external potential of Eq. (1).

## II Bulk Theory

We start by considering the bulk of a one component hard-sphere fluid in a volume $`V_{tot}`$ with $`V_{tot}^3`$ in the thermodynamic limit. The grand potential of this system is given by $`\mathrm{\Omega }_0=pV_{tot}`$, where $`p`$ is the pressure of the hard-sphere fluid. A single particle of a second component $`b`$ is inserted into this system. In view of the following discussion the only restriction for this particle is that it must be hard and convex. The grand potential of the new system is

$$\mathrm{\Omega }=\mathrm{\Omega }_0+\mathrm{\Delta }\mathrm{\Omega },$$
(17)
where $`\mathrm{\Delta }\mathrm{\Omega }`$ measures the change of the grand potential due to the insertion of a single particle. $`\mathrm{\Delta }\mathrm{\Omega }`$ equals the one-particle direct correlation function  which in the bulk limit considered here is simply the excess chemical potential of species $`b`$ in the dilute limit $`\rho _b0`$, i.e.,

$$\beta \mathrm{\Delta }\mathrm{\Omega }=c_b^{(1)}.$$
(18)
Although Eqs. (17) and (18) are formally exact it is important to notice that this approach within a bulk theory is expected to only work reliably if the considered fluid is in a single phase. In the case of phase separation, which could lead to wetting or drying of the big sphere, this approach most likely would fail .

We evaluate Eq. (18) within the framework of Rosenfeld’s Fundamental Measure Theory (FMT). Within this approach  the excess free energy functional (over the ideal gas) is given by

$$\beta _{ex}=d^3r\mathrm{\Phi }(\{n_\alpha \}),$$
(19)
where $`n_\alpha `$ denote weighted densities

$$n_\alpha (𝐫)=\underset{j=1}{\overset{N}{}}d^3r^{}\rho _j(𝐫^{})w_\alpha ^{(j)}(𝐫𝐫^{}),$$
(20)
with geometrical weight functions $`w_\alpha ^{(j)}`$ of species $`j=1,\mathrm{},N`$. In Rosenfeld’s approach there are six different weighted densities, four are scalar and two are vector-like. For a general hard body $`j`$ one has

$`w_3^{(j)}(𝐫)`$ $`=`$ $`\mathrm{\Theta }(|𝐫𝐑^{(j)}|),`$ (21)
$`w_2^{(j)}(𝐫)`$ $`=`$ $`\delta (|𝐫𝐑^{(j)}|),`$ (22)
$`𝐰_{V2}^{(j)}(𝐫)`$ $`=`$ $`\widehat{𝐧}^{(j)}\delta (|𝐫𝐑^{(j)}|),`$ (23)
$`w_1^{(j)}(𝐫)`$ $`=`$ $`{\displaystyle \frac{H^{(j)}}{4\pi }}w_2^{(j)}(𝐫),`$ (24)
$`w_0^{(j)}(𝐫)`$ $`=`$ $`{\displaystyle \frac{K^{(j)}}{4\pi }}w_2^{(j)}(𝐫),`$ (25)
$`𝐰_{V1}^{(j)}(𝐫)`$ $`=`$ $`{\displaystyle \frac{H^{(j)}}{4\pi }}𝐰_{V2}^{(j)}(𝐫),`$ (26)
where $`𝐑^{(j)}=𝐑^{(j)}(\theta ,\phi )`$ is the radius vector to the surface and $`\widehat{𝐧}^{(j)}`$ is its outward unit vector normal to the surface at point $`𝐑^{(j)}`$ of particle $`j`$. $`H^{(j)}`$ and $`K^{(j)}`$ are the integrated mean curvature and integrated Gaussian curvature of the particle $`j`$, respectively .

In the bulk limit for which the density profiles $`\rho _j(𝐫)`$ are constant, as considered in this section, the vector weighted densities vanish, and Eq. (18) takes the form

$$c_b^{(1)}=\frac{\mathrm{\Phi }}{n_3}\zeta _3+\frac{\mathrm{\Phi }}{n_2}\zeta _2+\frac{\mathrm{\Phi }}{n_1}\zeta _1+\frac{\mathrm{\Phi }}{n_0}\zeta _0,$$
(27)
with characteristic functions $`\zeta _i`$, $`i=0,\mathrm{},3`$ of the shape of the particle of species $`b`$.

Assuming the general functional of the given by Eqs. (19)-(26), one has for a general convex hard particle

$`\zeta _3`$ $`=`$ $`V^{(b)},`$ (28)
$`\zeta _2`$ $`=`$ $`A^{(b)},`$ (29)
$`\zeta _1`$ $`=`$ $`A^{(b)}{\displaystyle \frac{H^{(b)}}{4\pi }},`$ (30)
and

$$\zeta _0=A^{(b)}\frac{K^{(b)}}{4\pi },$$
(31)
where the shape of the particle is specified by its volume $`V^{(b)}`$, its surface area $`A^{(b)}`$, and integrated curvatures $`H^{(b)}`$ and $`K^{(b)}`$.

Within FMT there are several expressions for the excess free energy density $`\mathrm{\Phi }`$. For the given problem we have chosen the original Rosenfeld form

$$\mathrm{\Phi }(\{n_\alpha \})=n_0\mathrm{ln}(1n_3)+\frac{n_1n_2𝐧_1𝐧_2}{1n_3}+\frac{n_2^33n_2𝐧_2𝐧_2}{24\pi (1n_3)^2}.$$
(32)
The partial derivatives of the excess free energy density with respect to the weighted densities in the dilute limit of species $`b`$ depend only on the packing fraction $`\eta `$ and the radius $`R`$ of the hard-sphere fluid:

$`{\displaystyle \frac{\mathrm{\Phi }}{n_3}}`$ $`=`$ $`{\displaystyle \frac{3\eta }{4\pi R^3}}{\displaystyle \frac{1+\eta +\eta ^2}{(1\eta )^3}}\beta p_{PY},`$ (33)
$`{\displaystyle \frac{\mathrm{\Phi }}{n_2}}`$ $`=`$ $`{\displaystyle \frac{3\eta (2+\eta )}{8\pi R^2(1\eta )^2}}\beta \gamma _{HW}^{SPT},`$ (34)
$`{\displaystyle \frac{\mathrm{\Phi }}{n_1}}`$ $`=`$ $`{\displaystyle \frac{3\eta }{R(1\eta )}}\varphi _1,`$ (35)
and

$$\frac{\mathrm{\Phi }}{n_0}=\mathrm{ln}(1\eta )\varphi _0.$$
(36)
In the above equations $`p_{PY}`$ and $`\gamma _{HW}^{SPT}`$ are the Percus-Yevick compressibility equation of state and the scaled-particle theory (SPT) expression for the surface tension of a hard-sphere fluid at a planar hard wall , respectively.

With the characteristic functions in Eqs. (28)-(31) and the partial derivatives given by Eqs. (33)-(36) one finds within this approach for the grand potential of a hard-sphere fluid surrounding any convex hard particle b

$$\beta \mathrm{\Omega }=\beta p_{PY}(V_{tot}V^{(b)})+A^{(b)}\left(\beta \gamma _{HW}^{SPT}+\frac{H^{(b)}}{4\pi }\varphi _1+\frac{K^{(b)}}{4\pi }\varphi _0\right).$$
(37)
By choosing the surface of the convex particle as the dividing interface, the bulk contribution to the grand potential is simply $`\mathrm{\Omega }_0=p^{PY}(V_{tot}V^{(b)})`$ and the surface tension follows directly as

$$\beta \gamma (H^{(b)},K^{(b)})=\beta \gamma _{HW}^{SPT}+\frac{H^{(b)}}{4\pi }\varphi _1+\frac{K^{(b)}}{4\pi }\varphi _0.$$
(38)
The particular choice of the dividing interface underlying Eq. (38) is motivated by the simplicity of the resulting expression for the surface tension. Other choices of the dividing interface would lead to significantly more complex expressions and the elegant, simple form of Eq. (38) would be lost. Obviously within this bulk approach it is impossible to pick up a contribution to the surface tension that is proportional to the logarithm of the radius of curvature.

The surface tension in Eq. (38) is always positive for a non-zero packing fraction $`\eta `$ and its change for small curvature due to an increase in curvature is described by the Tolman length $`\delta _T`$ which within this approach is defined by

$$\gamma (H^{(b)},K^{(b)})=\gamma _{HW}^{SPT}\left(12H^{(b)}\delta _T+𝒪(K^{(b)})\right),$$
(39)
so that it is independent of the shape of the convex particle with

$$\delta _T=\frac{\varphi _1}{8\pi \beta \gamma _{HW}^{SPT}}=\frac{\eta 1}{\eta +2}R.$$
(40)
Thus within this approach we find that the generalized Tolman length for a general hard convex cavity is negative, provided the surface of this hard cavity is chosen as the dividing surface. The grand potential (Eq. (37)) and the surface tension (Eq. (38)) have been derived for a general convex hard wall and are related to the Helfrich Hamiltonian . In order to be able to calculate the contact density of a hard-sphere fluid on the basis of the sum rule in Eq. (7) one has to evaluate the expression in Eq. (38) for special cases. For a big spherical particle the mean curvature $`H_s^{(b)}=1/R_s^{(b)}`$ and the Gaussian curvature $`K_s^{(b)}=1/(R_s^{(b)})^2`$. By combining Eqs. (38) and (10) within the present bulk theory one finds for the contact value of the density profile of a hard-sphere fluid at a big hard sphere:

$$\rho (R_s^+)=\beta p_{PY}\frac{9\eta ^2}{4\pi (1\eta )^3R^3}\left\{(1+\eta )\frac{R}{R_s}\eta \frac{R^2}{R_s^2}\right\}.$$
(41)
For finite curvatures the contact density is smaller than in the planar limit, which in the present approach is given by the Percus-Yevick compressibility pressure $`\beta p_{PY}`$. The difference between the contact density at finite curvature and its planar wall limit is a quadratic function in $`1/R_s`$.

For a cylindrical wall we have to proceed more carefully. Since Eq. (38) was derived assuming a finite convex body we take the limit $`L\mathrm{}`$ only after the cylinder is inserted into the bulk system. In this limit the mean curvature $`H_c^{(b)}=1/(2R_c^{(b)})`$ and the Gaussian curvature $`K_c^{(b)}=0`$. The contact density of a hard-sphere fluid close to a cylindrical wall follows from Eq. (38) and the sum rule in Eq. (11) as

$$\rho (R_c^+)=\beta p_{PY}\frac{9\eta ^2}{4\pi (1\eta )^3R^3}\left\{\frac{(1+\eta )}{2}\frac{R}{R_c}\right\},$$
(42)
which is linear in $`1/R_c`$. For large radii of curvature $`R_i\mathrm{}`$ the contact densities at a big sphere $`\rho (R_s^+)`$ and at a cylinder $`\rho (R_c^+)`$ are related by $`2\rho (R_c^+)\rho (R_s^+)=\beta p_{PY}+𝒪(\frac{1}{R_s^2})`$.

## III Numerical results

As mentioned in the Introduction an alternative approach for calculating the surface tension and the contact density of a hard-sphere fluid at a spherical or cylindrical wall consists of minimizing the density functional, Eq. (3), numerically. Within this approach the curved wall enters the problem via the external potential specified in Eq. (1). The solution of the corresponding Euler-Lagrange equation yields the inhomogeneous equilibrium density profile $`\rho (r)`$. From $`\rho (r)`$ both the contact densities $`\rho (R_i^+)=\rho (r=R_i^+)`$ and the surface tensions $`\gamma _i(R_i^{(b)})`$, Eq. (5), follow directly. In order to verify that the numerical DFT results for the contact density and the surface tension (from minimization) are consistent and to probe the numerical accuracy of our DFT calculation, we confirmed numerically that the sum rule in Eq. (7) is satisfied very accurately (within the range of four to six significant digits depending on the packing fraction).

Although the minimization of the density functional and the bulk route discussed in Sec. II are in principal equivalent approaches, they differ in practice because the actual density functional for the hard-sphere fluid is an approximation. Within an approximate treatment it is expected that the numerical solution of the Euler-Lagrange equation, which is based on a variational principle, is more accurate than the aforementioned bulk approach. Nevertheless, for the hard-sphere fluid which, apart from freezing, exhibits only a single phase both routes are expected to be in qualitative or even semi-quantitative agreement with each other.

First we present in Fig. 2 our numerical results for the surface tension at a spherical wall as obtained by full minimization of the Rosenfeld functional for values of $`\eta `$ in the range between $`\eta =0.1`$ and $`\eta =0.42`$ and with radii of curvature $`R_s^{(b)}`$ in the range $`R_s^{(b)}=R`$ and $`R_s^{(b)}=10000R`$. The DFT results are denoted by symbols and the predictions of the bulk theory by full lines. For small values of $`\eta `$ we find both results to be in excellent, quantitative agreement for all curvatures, see Fig. 2(a). Close to the planar wall limit the level of quantitative agreement reduces significantly for larger packing fractions, while it remains good for large curvatures. This deviation between the two different routes is expected because in the limit $`R_s^{(b)}\mathrm{}`$ the bulk route reduces to the SPT prediction of the surface tension at a planar wall, which is known to overestimate the surface tension of a hard-sphere fluid at a planar hard wall . Since SPT obtains the planar wall limit from extrapolating the behavior at small cavity sizes  this deviation is no surprise. Recently it was suggested that SPT can be re-optimized by taking the planar wall limit into account more accurately .

Even for packing fractions at which the level of quantitative agreement is moderate (see Fig. 2(b)), there is a good qualitative agreement between the DFT results and our bulk approach. This implies that the curvature dependence of our numerical data of $`\gamma _s`$ can be best approximated by a polynomial quadratic in $`1/R_s^{(b)}`$ with no evidence for a logarithmic term.

The DFT results for the surface tension of a hard-sphere fluid in contact with a hard cylindrical wall are shown in Fig. 3. Again the quantitative agreement between the DFT results (symbols) and predictions of the bulk theory (full lines) are excellent for low packing fractions and fair at higher values of $`\eta `$. The curvature dependence of the DFT data is captured very well by a function linear in $`1/R_c^{(b)}`$ as predicted by the bulk route.

A remarkable feature of Eq. (38) is that the ratio of the coefficients corresponding to the term of the leading order in $`1/R_i`$ for spherical and cylindrical symmetry is exactly 2 (since $`H_s^{(b)}=2H_c^{(b)}`$). This behavior can also be found in the numerical results of DFT calculations, i.e., based on the inhomogeneous density profiles $`\rho (r)`$. In Table I we show values of least-square fits assuming that $`\beta \gamma _s(R_s^{(b)})=a_s^{(0)}+\frac{a_s^{(1)}}{R_s^{(b)}}+\frac{a_s^{(2)}}{(R_s^{(b)})^2}`$ and $`\gamma _c(R_c^{(b)})=a_c^{(0)}+\frac{a_c^{(1)}}{R_c^{(b)}}`$, i.e., the simplest formulae to fit our DFT data. To a good approximation the ratio $`a_s^{(1)}/a_c^{(1)}`$ equals 2 for all densities. This ratio agrees with the predictions of the bulk theory; $`a_s^{(0)}=a_c^{(0)}`$ is the planar wall surface tension. The Tolman length (compare Eq. (39)) is given by $`\delta _T^{(s)}=a_s^{(1)}/(2a_s^{(0)})`$ and $`\delta _T^{(c)}=a_c^{(1)}/a_c^{(0)}`$. In Fig. 4 we plot $`\delta _T^{(s)}`$ (diamonds) and $`\delta _T^{(c)}`$ (circles) together with the prediction of the bulk theory Eq. (40) (full line) as a function of the packing fraction $`\eta `$.

Thus neither the bulk theory nor the numerical DFT results provide any hints for the occurrence of logarithmic terms in the expression of the surface tension of hard-sphere fluids in contact with curved hard surfaces.

For the contact density one knows a priori that for all values of $`\eta `$ both the DFT route and the bulk approach yield the same contact value $`\rho (R_i^+\mathrm{})=\beta p`$ of the density profile of a hard-sphere fluid in the limit of the planar wall. For finite curvature we find that the contact densities obtained numerically from DFT (symbols) and analytically from our bulk theory for $`\gamma (R_s^{(b)})`$ (full lines) are in very good quantitative agreement for all values of $`\eta `$, as shown in Fig. 5 for spherical walls and in Fig. 6 for cylindrical walls. Only for very large curvatures the actual DFT contact densities are slightly overestimated by the bulk theory.

For the spherical wall in Fig. 5 we also plot the contact density obtained from an empirical fit to simulation data (dotted lines) by Degreve et al. . In terms of $`R_s^1`$ this fit interpolates linearly between the planar wall limit within the Carnahan-Starling theory and the contact value of the pair correlation function $`g(r)`$ corresponding to $`R/R_s=0.5`$. At high packing fractions this fit deviates from our results in the planar wall limit due to the difference between the Carnahan-Starling equation of state and the Percus-Yevick compressibility pressure which underlies both our approaches. At very high packing fractions there are additional deviations because the simple linear fit does not capture the actual higher order terms in $`R_s^1`$. Nonetheless, the overall agreement between our results and this fit is good.

## IV Conclusions and Summary

Our analysis of the local structure and the thermodynamics of the hard-sphere fluid near hard curved substrates has lead to the following main results:

1. We have studied the curvature dependence of the surface tension of a hard-sphere fluid at hard spherical and cylindrical walls (Fig. 1) obtained by minimizing the Rosenfeld fundamental measure density functional. We have found no indications of logarithmic singularities in the expansion of the surface tension in terms of curvatures of spherical or cylindrical walls. The surface tension appears to be an analytic function of curvature.
2. Based on the Rosenfeld fundamental measure theory we have derived an analytical expression for the surface tension of a hard-sphere fluid close to hard, arbitrarily curved, convex walls in terms of their integrated mean and Gaussian curvatures, $`H^{(b)}`$ and $`K^{(b)}`$ (see Eq. (38)). This approach also does not render logarithmic singularities of the surface tension as a function of the radii of curvature. There is good agreement between the results of the density functional theory and the bulk theory which treats the curved wall as the surface of a particle of a second component in infinite dilution (Figs. 2 and 3). The Tolman length as a function of the packing fraction of the hard-sphere fluid is shown in Fig. 4.
3. Based on the sum rule for curved substrates we have obtained expressions for the contact values of the density profile of a hard-sphere fluid close to spherical (Eq. (41)) and cylindrical (Eq.(42)) walls. A comparison with the density functional results indicates, that the analytical expressions of the contact values are very reliable; only at rather high packing fractions they slightly overestimate the density functional values (Figs. 5 and 6) for large curvatures.
4. The ratio of the leading order terms in the expansion of the surface tension as a function of curvature, obtained from full minimization of the fundamental measure functional, for cylindrical and spherical symmetry equals 2 (see Table I). This is in agreement with the aforementioned bulk theory (Eq. (38)). This also agrees with the general feature of the Helfrich theory , in which for a sphere the contribution to the surface tension linear in the curvature is twice as large as for a cylinder.

In Ref.  the arguments in favor of and against the existence of logarithmic singularities in the curvature dependence of the surface tension of fluids have been reviewed. Our present results indicate that for a hard-sphere fluid near hard curved substrates such logarithmic singularities either do not exist or their influence is negligible.

It would be of interest to investigate the issue of analyticity of the surface tension via a diagrammatic expansion of the contact value of the density profile . Within this approach our preliminary calculations show, that in the lowest order in density the diagrammatic expansion of the contact value of the density profile indeed leads to a surface tension analytical in curvature. However, the extension of this analysis to higher orders in densities is technically involved and is postponed to future work.

Thus our findings indicate that non-analyticities of the curvature dependence of surface tensions arise only via dispersion forces acting on the fluid particles  or via the onset of drying or wetting transitions on curved substrates

###### Acknowledgements.

It is a pleasure to thank Bob Evans for stimulating discussions.