text:
# 1 Introduction and summary

## 1 Introduction and summary

In view of the observation  that gravity may be localized on a brane embedded in space with one extra dimension of infinite size, it is of interest to study whether there exist mechanisms of (quasi-)localization of gravity in spaces with more than one infinite extra dimensions. One proposal of this sort has been put forward in Refs. . The basic idea  is that radiative effects due to matter residing on the brane may induce new terms in the effective action of multi-dimensional gravity (cf. Ref. ), which concentrate on the brane and dominate the gravitational interactions of brane matter. Thus, the effective action has the form

$$S_{tot}=S_{bulk}^{eff}+S_{brane}$$
(1)
Here the bulk term involves $`(4+N)`$-dimensional metric $`g_{AB}`$ ($`N>1`$ is the numer of extra dimensions) and at low energy reduces to the $`(4+N)`$-dimensional Einstein–Hilbert action111Leaving aside the issue of the cosmological constant. with the fundamental scale $`M_{}`$. The brane Einsten-Hilbert term, on the other hand, involves induced $`4`$-dimensional metric $`g_{\mu \nu }`$ on the brane and has its own mass scale $`M_{Pl}`$, which supposedly is determined by dynamics on the brane. It has been argued in Ref.  that the two scales may be completely different, and, in particular, that the relation

$$M_{}M_{Pl}$$
(2)
may hold.

For more than one extra dimensions, $`N>1`$, the model exhibits a potentially interesting UV-IR mixing. Naively, one would expect that at large distances along the brane, the relevant terms in $`S_{bulk}`$ and $`S_{brane}`$ are the multi-dimensional and four-dimensional Einstein-Hilbert terms respectively, while the brane may be treated as $`\delta `$-functional in transverse directions. This is not the case, however, because of the singularity of the $`N`$-dimensional propagator . Hence, the behavior of the model at large distances along the brane depends on how the singularity in transverse dimensions is resolved.

One way to resolve this singularity would be to smear the $`\delta `$-function in the brane action. This proposal, however, suffers the strong coupling problem at unacceptably low energies . Hence, we will not consider this option any longer.

Another proposal  which is the subject of this paper222Just for brevity, we will call this proposal as “brane induced gravity” in what follows., is that the bulk gravity is “soft” at distances shorter than $`M_{}^1`$. Under this assumption, matter on the brane experiences four-dimensional gravity at intermediate distances

$$\frac{1}{M_{}}rr_c\frac{M_{Pl}}{M_{}^2}$$
(3)
while four-dimensional Newton’s law ceases to hold at both short and ultra-large distances. It is worth noting that this multi-dimensional “brane-induced gravity” model, linearized about flat background, leads to pure tensor  four-dimensional gravity on the brane at intermediate distances (3), without an extra scalar inherent in the linearized brane-induced gravity in one extra dimension .

These features make brane-induced gravity with $`N>1`$ potentially interesting, in particular, from the viewpoint of the cosmological constant problem . The violation of four-dimensional Newton’s law at ultra-large distances, combined with the absence of extra scalar interaction on the brane at intermediate scales is alarming, however, as the same property was present in the model of Ref.  which has been found to have a ghost . Hence, brane-induced gravity in more than one extra dimensions is worth studying in some detail.

In this paper we consider brane-induced gravity, linearized about flat multi-dimensional space and flat brane, mostly at $`N>2`$; we discuss somewhat special case $`N=2`$ towards the end. In Section 2 we neglect complications due to tensor structure, and study a scalar counterpart of the model. We find that once the brane has finite thickness, the equivalence between gravitational and inertial masses is generally violated for matter on the brane, even in theory restricted to intermediate scales (3). This is again alarming, since in other models333Leaving aside models with extra light four-dimensional degrees of freedom., violation of “charge universality” (in the gravitational context, equivalence principle) is a signal for inconsistency .

We then proceed in Section 3 to brane-induced gravity itself. We study the linearized field equations, assuming first that the bulk term has the tensor structure of General Relativity. We begin with the study of low-mass states which are localized or quasi-localized near the brane. We find in Section 3.1 that one such state is a four-dimensional scalar; it is exactly localized on the brane and has negative (tachyonic) mass squared. Another state is a massive four-dimensional graviton444The graviton has finite, though very small width, $`\mathrm{\Gamma }_{graviton}m_{graviton}`$, i.e., it is, strictly speaking, quasi-localized.. Both masses are of order

$$|m_{tachyon}||m_{graviton}|r_c^1\frac{M_{}^2}{M_{Pl}}$$

In Section 3.2 we proceed to show that the tachyon is actually a ghost. This can be seen in two ways. One is to study the propagator of the full linearized theory near the tachyon pole and show that the residue has negative sign. Another way is to evaluate the propagator from brane to brane, which describes gravitational interaction of matter on the brane. We find that the brane-to-brane propagator is a sum of two terms, one of which has a pole at $`p^2=m_{graviton}^2`$ with tensor structure appropriate for massive graviton, while another is a scalar ghost term (of overall negative sign) with a pole at $`p^2=m_{tachyon}^2`$. This situation is completely analogous to that in the model of Ref. : at intermediate scales (3), the ghost term cancels out the extra  scalar part of the massive graviton propagator, so that the brane-to-brane propagator at intermediate distances has massless tensor form.

We comment on the case of two extra dimensions, $`N=2`$, in Section 3.3. There are peculiarities, but the outcome is the same: the model has tachyonic ghost.

In Section 4 we generalize by allowing for most general tensor structure of the linearized bulk equations (in fact, there are only two terms consistent with $`(4+N)`$-dimensional general covariance). We again study the case $`N>2`$, and evaluate the brane-to-brane propagator. We find that it again has ghost term, although the mass of the ghost is no longer necessarily tachyonic.

Our overall conclusion is that the linearized brane-induced gravity as it stands has a ghost, if the number of extra dimensions is larger than one. We interprete this property as an indication that this version of induced gravity cannot emerge as a low energy limit of any consistent microscopic theory.

It is worth noting that the low energy action of the general form (1) emerges in string theory context . Furthemore, the hierarchy (2) is also possible in string theory framework . It would be of interest to understand how string theory resolves the UV-IR ambiguity inherent in the case of more than one extra dimensions.

## 2 Scalar model

We begin with a counterpart of the brane-induced gravity with metric perturbations mimicked by a single scalar field $`\mathrm{\Phi }`$. In what follows it will be convenient to consider thick brane, and take the limit of delta-function brane in the end of calculations, if desired. It has been argued in Ref.  that the loops (and/or non-perturbative effects) involving matter on the brane induce non-local terms in the effective action, with the scale of non-locality set by the brane thickness $`\mathrm{\Delta }`$. At quadratic level, this effect is modelled by the induced action of the following form ,

$$S_{brane}^{(2)}=\frac{M_{Pl}^2}{2}d^4xd^Nyd^Ny^{}f^2(y)_\mu \mathrm{\Phi }(x,y)f^2(y^{})^\mu \mathrm{\Phi }(x,y^{})$$
(4)
where $`f(y)`$ is a smooth function localized near the brane; it accounts for thickness of the brane. $`N`$ is the number of extra dimensions; we concentrate on the case $`N>2`$. Without loss of generality, $`f`$ is normalized to unity,

$$d^Nyf^2(y)=1$$
(5)
and is non-zero in a region of size of order $`\mathrm{\Delta }`$. Hereafter $`X_A=(x_\mu ,y_a)`$ are coordinates in $`(4+N)`$ dimensions, $`\mu =0,\mathrm{},3`$; $`a=4,\mathrm{}N+3`$; signature of metric is mostly negative.

Let the bulk theory have the effective action $`S_{bulk}^{eff}[\mathrm{\Phi }]`$. There are two more assumptions in the model : (i) The mass scale entering $`S_{bulk}^{eff}`$ is $`M_{}`$ which is much smaller than $`M_{Pl}`$; (ii) The bulk theory is “soft” at length scales below $`1/M_{}`$, which we understand as the assumption that the Green’s functions of the bulk theory rapidly vanish at high Euclidean momenta.

For considering linearized theory (weak sources), let us neglect the non-linear terms in the bulk effective action. Then the only relevant term in $`S_{bulk}^{eff}`$ is quadratic in $`\mathrm{\Phi }`$, and has the form

$$S_{bulk}^{eff,(2)}=\frac{1}{2}d^{N+4}X\mathrm{\Phi }(X)(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}\mathrm{\Phi }(X)$$

where $`(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}`$ is the exact inverse propagator of the bulk theory. At low energies, the form factor $``$ is a constant of order $`M_{}^{2+N}`$ (note that the field $`\mathrm{\Phi }`$ is dimensionless). Let us denote the exact propagator of the bulk theory by $`D_{}(XX^{})`$, so that

$$D_{}(P)=\frac{1}{P^2(P^2)}$$

where $`P^2=p^2p_y^2`$, and $`p^2=p_\mu p^\mu `$. At momenta below $`M_{}`$, the propagator $`D_{}`$ coincides with the free propagator,

$$D_{}(P)=\frac{1}{M_{}^{2+N}}\frac{1}{P^2},|P^2|M_{}^2$$
(6)
and, by assumption of softness, $`D_{}(P)`$ rapidly tends to zero at large negative (Euclidean) $`P^2`$, with characteristic scale $`M_{}`$.

### 2.1 Scalar propagator

Let the source on the brane be characterized by a spread function $`g^2(y)`$, and be $`\delta `$-function in $`x`$-coordinates, where, again without loss of generality, $`g`$ is normalized to unity,

$$d^Nyg^2(y)=1$$
(7)
It is convenient to work in mixed representation, momentum in four dimensions and coordinate in extra dimesnions. One has the following equation for the propagator from brane to everywhere for given shape of the source,

$$(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}G_g(p,y^{})+M_{Pl}^2p^2f^2(y)d^Ny^{}f^2(y^{})G_g(p,y^{})=g^2(y)$$
(8)
where $`\mathrm{}^{(4+N)}=p^2_y^2`$.

Dvali et.al.  proceed under assumption that $`g^2(y)=f^2(y)`$ with corrections suppressed by $`M_{}/M_{Pl}`$. Let us drop this assumption, and see what happens.

Equation (8) has the following solution,

$$G_g(p,y)=D_g(p,y)\frac{M_{Pl}^2p^2D_{fg}(p)}{1+M_{Pl}^2p^2D_{ff}(p)}D_f(p,y)$$
(9)
where for any function $`u(y)`$ one defines

$$D_u(p,y)=d^Ny^{}D_{}(p,yy^{})u^2(y^{})$$
(10)
and for two functions $`u(y)`$, $`v(y)`$ one writes

$$D_{uv}(p)=D_{vu}(p)=d^Nyd^Ny^{}D_{}(p,yy^{})u^2(y^{})v^2(y)$$
(11)
Let us rewrite the expression (9) in the following suggestive form,

$`G_g(p,y)`$ $`=`$ $`{\displaystyle \frac{D_g(p,y)}{1+M_{Pl}^2p^2D_{ff}(p)}}`$ (12)
$`+`$ $`{\displaystyle \frac{M_{Pl}^2p^2[D_{ff}(p)D_g(p,y)D_{fg}(p)D_f(p,y)]}{1+M_{Pl}^2p^2D_{ff}(p)}}`$
Now, recall that

$$D_{}(p,yy^{})=d^Np_yD_{}(p^2p_y^2)\text{e}^{ip_y(yy^{})}$$

To evaluate the integral of the form $`D_{uv}(p)`$, we assume that the brane thickness $`\mathrm{\Delta }`$ is much smaller than $`1/M_{}`$, and write for small $`y`$ and $`y^{}`$

$$D_{}(p,yy^{})=D_{}^{(0)}(p)+D_{}^{(2)}(p)(yy^{})^2+\mathrm{}$$

Clearly,

$$D_{}^{(0)}(p)=d^Np_yD_{}(p^2p_y^2)$$

and

$$D_{}^{(2)}(p)=\frac{1}{2N}d^Np_yp_y^2D_{}(p^2p_y^2)$$
(13)
We assume that the latter integrals are convergent at negative (Euclidean) four-momenta, $`p^20`$, because of softness of the propagator $`D_{}`$ at short distances. On dimensional grounds

$$D_{}^{(0)}(|p|M_{})\frac{1}{M_{}^4}$$
(14)
and

$$D_{}^{(2)}(|p|M_{})\frac{1}{M_{}^2}$$
(15)
To the first non-trivial order in brane thickness, one has (assuming that $`u`$ and $`v`$ are normalized to unity, see eqs. (5) and (7))

$$D_{uv}(p)=D_{}^{(0)}(p)+D_{}^{(2)}(p)\mathrm{\Delta }_{uv}^2$$
(16)
where

$$\mathrm{\Delta }_{uv}^2=d^Nyd^Ny^{}(yy^{})^2u^2(y)v^2(y^{})$$

explicitly depends on the shapes of the functions $`u(y)`$ and $`v(y)`$ and is generically of the order of $`\mathrm{\Delta }^2`$.

At low momenta, $`|p|M_{}`$, one can set $`D_{ff}=\text{const}M_{}^4`$ in the denominators in eq. (12). Then at intermediate distances (3), virtuality $`p^2`$ is large enough, and one has

$`G_g(p,y)`$ $`=`$ $`{\displaystyle \frac{D_g(p,y)}{M_{Pl}^2p^2D_{ff}(p)}}`$ (17)
$`+`$ $`{\displaystyle \frac{D_{ff}(p)D_g(p,y)D_{fg}(p)D_f(p,y)}{D_{ff}(p)}}`$
This propagator determines the field induced by a weak source of shape $`g^2(y)`$ in transverse directions, in the theory restricted to intermediate scales (3).

### 2.2 Potential on the brane

The interaction between sources with spread functions $`g^2(y)`$ and $`h^2(y)`$ is described by the effective four-dimenstional propagator, which is the convolution of $`G_g(p,y)`$ and $`h^2(y)`$. At intermediate values of momenta, $`M_{}^2p^2r_c^2`$, one has from eq. (17)

$$G_{eff,4d}(p)=\frac{D_{gh}(p)}{M_{Pl}^2p^2D_{ff}(p)}+\frac{D_{gh}(p)D_{ff}(p)D_{fg}(p)D_{fh}(p)}{D_{ff}(p)}$$
(18)
Keeping terms of order $`\mathrm{\Delta }^2`$, one finds

$`G_{eff,4d}(p)`$ $`=`$ $`{\displaystyle \frac{D_{}^{(0)}(p)+D_{}^{(2)}(p)(\mathrm{\Delta }_{gh}^2\mathrm{\Delta }_{ff}^2)}{M_{Pl}^2p^2D_{}^{(0)}(p)}}`$ (19)
$`+`$ $`D_{}^{(2)}(p)(\mathrm{\Delta }_{gh}^2+\mathrm{\Delta }_{ff}^2\mathrm{\Delta }_{fg}^2\mathrm{\Delta }_{fh}^2)`$
Consider the first term. Because of explicit $`p^2`$ in denominator, one can replace $`D_{}^{(0)}`$ and $`D_{}^{(2)}`$ by constants at $`pM_{}`$, i.e. at distances larger than $`M_{}^1`$. This leads to four-dimensional Newton’s potential with non-universal gravitational constant

$$G_{Newton,eff}=\frac{1}{M_{Pl}^2}\left[1+\frac{D_{}^{(2)}(0)}{D_{}^{(0)}(0)}(\mathrm{\Delta }_{gh}^2\mathrm{\Delta }_{ff}^2)\right]$$

According to eqs. (14) and (15), the non-universal correction is of order $`\mathrm{\Delta }^2M_{}^2`$.

This is the main result of this Section: the model allows for (weak) violation of the equivalence principle, since the spread functions $`g^2(y)`$ and $`h^2(y)`$ may have different shapes, depending on the type of matter residing on the brane.

The second term in eq. (19) corresponds to short-ranged force. According to eq. (13), one has, in coordinate representation,

$$D_{}^{(2)}(x)=\frac{1}{2N}_y^2D_{}(x^2y^2)|_{y=0}$$

At relatively large distances, $`rM_{}^1`$, the propagator $`D_{}`$ is a free propagator in $`(4+N)`$ dimensions, up to a factor $`1/M_{}^{2+N}`$. This gives

$$𝑑x^0D_{}^{(2)}(x)=\frac{1}{M_{}^{2+N}|𝐱|^{3+N}},|𝐱|M_{}^1$$

up to numerical constant of order one. Hence, the correction to Newton’s potential is

$$\mathrm{\Delta }V(r)=\frac{1}{M_{}^{2+N}r^{3+N}}(\mathrm{\Delta }_{gh}^2+\mathrm{\Delta }_{ff}^2\mathrm{\Delta }_{fg}^2\mathrm{\Delta }_{fh}^2)$$

This is a short-ranged potential, “fifth force”, which again depends on the composition of matter (functions $`g^2`$ and $`h^2`$). It is worth noting that the latter non-universality exists also at $`N=1`$ , where the brane-induced gravity does not show any inconsistency.

## 3 Tachyonic ghost

Let us now consider the linearized brane-induced gravity and keep track of its tensor structure. In this Section we assume for simplicity that the tensor structure of the linearized equations in the bulk coincides with that in the linearized Einstein theory in $`(4+N)`$ dimensions. Then the linearized field equation takes the following form,

$$(\mathrm{}^{(4+N)})G_{AB}(x,y)+M_{Pl}^2f^2(y)𝑑y^{}f^2(y^{})G_{AB}^{(4)}(x,y^{})=T_{AB}(x,y)$$
(20)
where $`G_{AB}=R_{AB}(1/2)g_{AB}R`$ is the linearized Einstein tensor in $`(4+N)`$ dimensions, $`G_{aB}^{(4)}=0`$, the four-dimensional Einstein tensor $`G_{\mu \nu }^{(4)}`$ is constructed in terms of four-dimensional components of the metric. The form-factor $``$ has the same properties as above. The function $`f^2(y)`$ is again the spread function for induced term.

Let us impose the harmonic gauge

$$_Ah_B^A=\frac{1}{2}_Bh_A^A$$
(21)
where $`h_{AB}`$ are perturbations about Minkowski metric $`\eta _{AB}`$; indices are raised and lowered by Minkowski metric. Then one has

$$G_{AB}=\frac{1}{2}\mathrm{}^{(4+N)}(h_{AB}\frac{1}{2}\eta _{AB}h_D^D)$$
(22)
while $`G_{\mu \nu }^{(4)}`$ remains in its general form

$$G_{\mu \nu }^{(4)}=\frac{1}{2}[_\mu _\lambda h_\nu ^\lambda +_\nu _\lambda h_\mu ^\lambda \mathrm{}^{(4)}h_{\mu \nu }_\mu _\nu h_\lambda ^\lambda \eta _{\mu \nu }(_\lambda _\rho h^{\lambda \rho }\mathrm{}^{(4)}h_\lambda ^\lambda )]$$
(23)
Hereafter $`\mathrm{}^{(4)}=_\mu ^\mu `$.

### 3.1 (Quasi-)localized states: tachyon and massive graviton

Let us consider the sourceless field equation, i.e., eq. (20) with $`T_{AB}=0`$, to see whether there exist modes which are (quasi-)localized near the brane. The $`(ab)`$ and $`(a\mu )`$-components of this equation in the gauge (21) read

$$\frac{1}{2}(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}(h_{ab}\frac{1}{2}\eta _{ab}h_D^D)=0$$

$$\frac{1}{2}(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}h_{a\mu }=0$$

These are the equations of the bulk theory for the corresponding combinations of metrics, and they do not have localized solutions. Hence,

$$h_{ab}=\frac{1}{2}\eta _{ab}h_C^C$$
(24)
and

$$h_{a\mu }=0$$
(25)
After taking trace of eq. (24), one expresses the $`(ab)`$-components of the metric in terms of the trace of the four-components,

$$h_{ab}=\frac{1}{N2}\eta _{ab}h_\mu ^\mu $$
(26)
(at this point we specialize to $`N>2`$). Then the gauge condition (21) with $`B=\mu `$, together with eq. (25) give

$$_\mu h_\nu ^\mu =\frac{1}{N2}_\nu h_\lambda ^\lambda $$

Making use of the above relations, one obtains for the remaining, $`(\mu \nu )`$-components of the field equations

$``$ $`{\displaystyle \frac{1}{2}}(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}\left(h_{\mu \nu }+{\displaystyle \frac{1}{N2}}\eta _{\mu \nu }h_\lambda ^\lambda \right)`$ (27)
$`+`$ $`{\displaystyle \frac{M_{Pl}}{2}}f^2(y){\displaystyle 𝑑y^{}f^2\left[\frac{N}{N2}_\mu _\nu h_\lambda ^\lambda \mathrm{}^{(4)}h_{\mu \nu }+\frac{N1}{N2}\eta _{\mu \nu }\mathrm{}^{(4)}h_\lambda ^\lambda \right]}=0`$
Trace of this equation gives

$$(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}h_\mu ^\mu +\stackrel{~}{M}_{Pl}^2f^2(y)𝑑y^{}f^2\mathrm{}^{(4)}h_\mu ^\mu =0$$
(28)
where

$$\stackrel{~}{M}_{Pl}^2=\frac{2(N1)}{N+2}M_{Pl}^2$$
(29)
The latter is a scalar equation, and we are interested in its solution localized near the brane. This solution is expressed in terms of functions $`D_f(p,y)`$ and $`D_{ff}(p)`$ introduced in eqs. (10) and (11). In the mixed representation the solution is

$$h_\mu ^\mu (y)=cD_f(p^2=m_{}^2,y)$$
(30)
where $`c`$ is a normalization constant and the mass is determined by the “eigenvalue equation”

$$m_{}^2=\frac{1}{\stackrel{~}{M}_{Pl}^2D_{ff}(m_{}^2)}$$
(31)
Let us see that the mass squared, $`m_{}^2`$, is, in fact, negative and real,

$$m_{tachyon}^2m_{}^2<0$$
(32)
$$\text{Im}(m_{}^2)=0$$
(33)
so the mode we consider is a tachyon localized near the brane. We first note, that

$$|m_{}|\frac{M_{}^2}{M_{Pl}}r_c^1$$

which is small compared to $`M_{}`$. Now, one has

$$D_{ff}(p^2)=d^Np_y\frac{|f^2(p_y)|^2}{(P^2)(P^2)}$$
(34)
where $`P^2=p^2p_y^2`$, as before. Since one assumes that the form-factor $``$ rapidly grows at large negative $`P^2`$ (the propagator $`D_{}`$ rapidly decays), this integral is convergent in the ulraviolet, and the integrand does not have singularities at $`p^2<0`$ (a zero of $`(P^2)`$ at negative $`P^2`$ would imply that there is a tachyon in bulk theory). For $`N>2`$ the integral is infrared-convergent even at $`p^2=0`$, since for $`|P^2|M_{}^2`$ the form-factor $``$ is constant. For small $`p^2<0`$ the integral here is a real positive constant, which is of order $`M_{}^4`$ on dimensional grounds, so $`D_{ff}(p^2=|m_{}^2|)`$ is a negative constant at small $`|m_{}|`$. One concludes that, as long as scales lower than $`M_{}`$ are concerned, there exists a single solution to eq. (31) which indeed obeys (32), (33).

Finally, we have to show that the wave function (30) decays as $`|y|\mathrm{}`$. One writes

$$D_f(p,y)=d^Np_y\frac{\text{e}^{ip_yy}f^2(p_y)}{P^2(P^2)}$$

Large $`|y|`$ correspond to small $`p_y`$, so at large $`|y|`$ one has

$$D_f(p^2=m_{}^2,y)d^Np_y\frac{\text{e}^{ip_yy}f^2(p_y)}{p_y^2+|m_{}^2|}$$
(35)
Recalling that $`f^2(p_y=0)=d^Nyf^2(y)=1`$, one obtains that the wave function (35) has the shape of $`N`$-dimensional Yukawa potential with (small) mass $`|m_{}|`$. Hence the wave function indeed decays as $`|y|\mathrm{}`$.

To obtain the complete tensor structure of the tachyon mode, one plugs the solution for the trace, eq. (30), back into eq. (27), and obtains, in mixed representation,

$``$ $`(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}h_{\mu \nu }+M_{Pl}^2f^2(y)p^2{\displaystyle 𝑑y^{}f^2(y^{})h_{\mu \nu }(p,y^{})}`$ (36)
$`=`$ $`cM_{Pl}^2D_{ff}p^2\left[{\displaystyle \frac{N}{N2}}{\displaystyle \frac{p_\mu p_\nu }{p^2}}{\displaystyle \frac{N(N1)}{(N2)(N+2)}}\eta _{\mu \nu }\right]D_f(p,y)`$
where $`p^2=m_{}^2`$ for the mode we discuss. This inhomogeneous equation is readily solved. The $`ab`$-components are found from eq. (26). In this way one finds the complete expression for the (unnormalized) tachyon mode

$`h_{\mu \nu }^{(m_{})}`$ $`=`$ $`{\displaystyle \frac{1}{3}}\left(\eta _{\mu \nu }{\displaystyle \frac{N+2}{N1}}{\displaystyle \frac{p_\mu p_\nu }{p^2}}\right)D_f(y)`$
$`h_{\mu a}^{(m_{})}`$ $`=`$ $`0`$
$`h_{ab}^{(m_{})}`$ $`=`$ $`{\displaystyle \frac{1}{N1}}\eta _{ab}D_f(y)`$ (37)
where $`p^2=m_{}^2<0`$ and $`D_f(y)=D_f(p^2=m_{}^2;y)`$.

For completeness, let us consider (quasi-)localized traceless modes, for which $`h_\mu ^\mu =0`$. For these modes, one obtains from eq. (27) the following equation,

$$(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}h_{\mu \nu }+M_{Pl}^2f^2(y)p^2𝑑y^{}f^2(y^{})h_{\mu \nu }(p,y^{})=0$$

The solution to this equation is again of the form

$$h_{\mu \nu }(y)=c_{\mu \nu }D_f(p^2=m^2,y)$$

where $`c_{\mu \nu }`$ are independent of $`y`$, and the mass now obeys

$$m^2=\frac{1}{M_{Pl}^2D_{ff}(m^2)}$$

We are interested in solutions with $`|m|M_{}`$, which are relevant at low energies. According to eq. (34), for such a solution the real part of $`m^2`$ is positive, and is of order $`r_c^2`$. Now, for small positive $`p^2`$, the function $`D_{ff}(p^2)`$ has even smaller imaginary part, which may be estimated as follows. The integrand in eq. (34) is a smooth positive function at $`p_y^2p^2`$, so this region does not contribute to the imaginary part. The imaginary part comes from the infrared region, and is proportional to

$$_0^ϵ\frac{p_y^{N1}dp_y}{p_y^2p^2i0}$$

The imaginary part of the latter integral is proportional to

$$ip^{N2}$$
(38)
So, for $`N>2`$ the quasi-localized graviton has a small mass $`m_{graviton}m`$, where

$$\text{Re}(m)=\frac{1}{M_{Pl}\sqrt{|D_{ff}(0)|}}\frac{M_{}^2}{M_{Pl}}r_c^1$$

and even smaller width,

$$\frac{\mathrm{\Gamma }_{graviton}}{m}\frac{m^{N2}}{M_{}^{N2}}$$

We conclude that in this model, there is a massive four-dimensional graviton with a tiny width. The violation of Newton’s law at distances of order $`r_c`$ is due to the graviton mass, not width, in clear contrast to the five-dimensional case .

### 3.2 Propagators at low energies: the tachyon is a ghost

One way to see that the tachyon is in fact a ghost, is to calculate the full propagator $`D_{AB,CD}(p;y,y^{})`$ near $`p^2=m_{}^2`$, i.e., extract its pole term. This is done in Appendix A. The outcome is

$$D_{AB,CD}^{(pole)}(p;y,y^{})=\frac{3}{M_{Pl}^2[D_{ff}(m_{}^2)]^2}\frac{h_{AB}^{(m_{})}(y)h_{CD}^{(m_{})}(y)}{p^2m_{}^2}$$
(39)
where $`h_{AB}^{(m_{})}(y)`$ is the (unnormalized) tachyon wave function (37). The overall negative sign here means that the tachyon is indeed a ghost.

The structure of the pole term (39) is precisely what one expects for the contribution of a mode localized near the brane. From eq. (39) one deduces also that the properly normalized tachyon-ghost mode is

$$h_{AB}^{normalized}(y)=\frac{\sqrt{3}}{M_{Pl}|D_{ff}|}h_{AB}^{(m_{})}$$

One observes from the latter formula and eq. (37) that the tachyonic ghost couples to matter on the brane at gravitational strength.

It is perhaps more instructive to study the propagator with both end-points on the brane. More precisely, let us consider the source on the brane with the only non-vanishing components $`T_{\mu \nu }`$, which is distributed in the transverse directions with the same555The analysis of the general case of a source with spread function $`g^2`$ different from $`f^2`$ proceeds along the lines of Section 2. This analysis is straightforward but not illuminating. spread function $`f^2(y)`$ as in eq. (20),

$$T_{\mu \nu }(x,y)=\theta _{\mu \nu }(x)f^2(y),$$
(40)
where $`\theta _{\mu \nu }(x)`$ is conserved in the four-dimensional sense. The point is to calculate $`(\mu \nu )`$-components of the metric due to this source. This is done in Appendix B, with the result, in mixed representation,

$`h_{\mu \nu }(p,y)`$ $`=`$ $`{\displaystyle \frac{2}{M_{Pl}^2}}{\displaystyle \frac{D_f(p,y)}{D_{ff}(p)}}\left[{\displaystyle \frac{1}{p^2m^2(p)}}\left(\theta _{\mu \nu }{\displaystyle \frac{1}{3}}\eta _{\mu \nu }\theta _\lambda ^\lambda \right){\displaystyle \frac{1}{6}}{\displaystyle \frac{1}{p^2m_{}^2(p)}}\eta _{\mu \nu }\theta _\lambda ^\lambda \right]`$ (41)
$`+`$ longitudinal part
where

$$m^2(p)=\frac{1}{M_{Pl}^2D_{ff}(p)}$$
(42)
$$m_{}^2(p)=\frac{1}{\stackrel{~}{M}_{Pl}^2D_{ff}(p)}$$
(43)
and the longitudinal part is proportional to $`p_\mu p_\nu `$ and vanishes when contracted with conserved stress-energy (the overall factor $`2`$ in (41) is due to our definition of $`M_{Pl}`$, see eq. (20)). Now, the interaction between two sources of the form (40) may be written in terms of the effective four-dimensional propagator $`D_{\mu \nu ,\lambda \rho }^{(4)}(p)`$, so that one has

$$\theta _{\mu \nu }^{}(p)D_{\mu \nu ,\lambda \rho }^{(4)}(p)\theta _{\lambda \rho }(p)=\theta _{\mu \nu }^{}(p)d^Nyf^2(y)h_{\mu \nu }(y,p),$$

where $`h_{\mu \nu }`$ is given by eq. (41). Hence, the effective brane-to-brane propagator is

$`D_{\mu \nu ,\lambda \rho }^{(4)}`$ $`=`$ $`{\displaystyle \frac{2}{M_{Pl}^2}}\left[{\displaystyle \frac{1}{p^2m^2(p)}}\left({\displaystyle \frac{1}{2}}(\eta _{\mu \lambda }\eta _{\nu \rho }+\eta _{\mu \rho }\eta _{\nu \lambda }){\displaystyle \frac{1}{3}}\eta _{\mu \nu }\eta _{\lambda \rho }\right){\displaystyle \frac{1}{6}}{\displaystyle \frac{1}{p^2m_{}^2(p)}}\eta _{\mu \nu }\eta _{\lambda \rho }\right]`$ (44)
$`+`$ longitudinal part
At low energies the “masses” $`m^2(p)`$ and $`m_{}^2(p)`$ are constants (up to tiny $`p`$-dependent imaginary part, see eq. (38)). Thus, at low energies the propagator (44) corresponds to a massive graviton of mass $`m`$ (note that the Van Dam–Veltman–Zakharov property indeed holds) and tachyonic ghost with negative $`m_{}^2`$. This ghost cancels the contribution of the extra graviton polarization at intermediate momenta $`M_{}|p|(m,m_{})r_c^1`$, so that at these scales the brane-to-brane propagator has the same structure as in General Relativity,

$$D_{\mu \nu ,\lambda \rho }^{(4)}=\frac{1}{M_{Pl}^2}\frac{1}{p^2}(\eta _{\mu \lambda }\eta _{\nu \rho }+\eta _{\mu \rho }\eta _{\nu \lambda }\eta _{\mu \nu }\eta _{\lambda \rho })+\text{longitudinal  part}$$

This is precisely the same situation as in the model of Ref. : the reason of why the correct tensor structure emerges in the linearized theory at intermediate distances is the existence of a ghost field.

### 3.3 N=2

The case $`N=2`$ is somewhat special. Let us first consider the tachyon mode. Equation (24) implies now

$$h_\mu ^\mu =0$$
(45)
while $`h_a^a`$ is arbitrary at this point. The four-dimensional trace of the sourceless equation (20) gives then

$$(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}h_a^a+\stackrel{~}{M}_{Pl}^2f^2(y)𝑑y^{}f^2\mathrm{}^{(4)}h_a^a=0$$

This equation has the same structure as eq. (28), so there again exists a tachyon. At $`N=2`$, it is the extra-dimensional metric $`h_{ab}`$ and traceless part of $`h_{\mu \nu }`$ that do not vanish in the tachyon mode (in the gauge (21)).

Another point is that the integral (34) is logarithmic at $`N=2`$, so the estimate for the graviton and tachyon masses is now

$$\text{Re}(m^2)|m_{}^2|\frac{M_{}^4}{M_{Pl}^2}\mathrm{log}\frac{M_{Pl}}{M_{}}$$

The imaginary part of the graviton mass is suppressed relative to its real part by logarithm only,

$$\frac{\mathrm{\Gamma }_{grav}}{|m|}\frac{1}{\mathrm{log}\frac{M_{Pl}}{M_{}}}$$

Yet the graviton width is smaller than its mass.

The tachyon is a ghost at $`N=2`$ as well. A simple way to see this is to redo the calculation leading to the brane-to-brane propagator. One finds that the expression (44) remains valid at $`N=2`$, the property (45) being ensured by the appropriate structure of the longitudinal terms. The negative sign of the last term on the right hand side of eq. (44) tells that the tachyon is indeed a ghost.

So, in spite of peculiarities, the conclusion for $`N=2`$ is the same as for $`N>2`$: the model has a tachyonic ghost.

## 4 Generalized model

In this section we drop the assumption that the tensor structure of the linearized bulk equations coincides with that in the linearized Einstein theory and consider the most general tensor structure compatible with the $`(4+N)`$-dimensional general covariance. The linearized field equation in the bulk theory has the following general form,

$$𝒟_{ABCD}h^{CD}=0$$
(46)
with some linear operator $`𝒟_{ABCD}`$. The symmetry of this operator under $`AB`$, $`CD`$ and $`(AB)(CD)`$ implies the following structure of $`𝒟_{ABCD}`$

$`𝒟_{ABCD}=a_A_B_C_D+b\left(_A_B\eta _{CD}+\eta _{AB}_C_D\right)`$
$`+c\left(_A\eta _{BC}_D+_B\eta _{AC}_D+_A\eta _{BD}_C+_B\eta _{AD}_C\right)`$
$`+d\eta _{AB}\eta _{CD}+e\left(\eta _{AC}\eta _{BD}+\eta _{AD}\eta _{BC}\right),`$
where $`a,b,c,d,e`$ are yet arbitrary functions of $`\mathrm{}^{(N+4)}`$. Now, gauge invariance implies

$$_A𝒟_{BCD}^A=0,$$

This leaves only two possible tensor structures which may appear in $`𝒟_{ABCD}`$, namely the usual Einstein structure and the product of two projectors, $`𝒟_{ABCD}P_{AB}P_{CD}`$, where

$$P_{AB}=\mathrm{}^{(4+N)}\eta _{AB}_A_B$$

In the harmonic gauge (21) one has

$$P_{AB}P_{CD}h^{CD}\mathrm{\Pi }_{AB}=\frac{1}{2}(\mathrm{}^{(4+N)}\eta _{AB}_A_B)\mathrm{}^{(4+N)}h_C^C.$$

Consequently, the generalization of eq. (20) is

$$(\mathrm{}^{(4+N)})G_{AB}(x,y)+𝒢(\mathrm{}^{(4+N)})\mathrm{\Pi }_{AB}(x,y)+M_{Pl}^2f^2(y)𝑑y^{}f^2(y^{})G_{AB}^{(4)}(x,y^{})=T_{AB}(x,y),$$
(47)
where $`G_{AB}`$ is given by eq. (22), and a new form-factor $`𝒢`$ is assumed to have the same ultraviolet properties as the form-factor $``$.

We will not repeat all the steps of the analysis of Section 3. To see the existence of the localized ghost and find its wave function in this general setup, it suffices to study the structure of the brane-to-brane propagator. We again consider a source of the form (40) and evaluate the $`(\mu \nu )`$-components of the metric induced by this source. The result is (see Appendix C for calculational details)

$`h_{\mu \nu }(p,y)`$ $`=`$ $`{\displaystyle \frac{2}{M_{Pl}^2}}\left[{\displaystyle \frac{D_f(p,y)}{D_{ff}(p)}}{\displaystyle \frac{1}{p^2m^2(p)}}\left(\theta _{\mu \nu }{\displaystyle \frac{1}{3}}\eta _{\mu \nu }\theta _\lambda ^\lambda \right){\displaystyle \frac{\stackrel{~}{D}_f(p,y)}{\stackrel{~}{D}_{ff}}}{\displaystyle \frac{1}{6}}{\displaystyle \frac{1}{p^2m_{}^2(p)}}\eta _{\mu \nu }\theta _\lambda ^\lambda \right]`$ (48)
$`+`$ longitudinal part
where $`\stackrel{~}{D}_f`$ and $`\stackrel{~}{D}_{ff}`$ are defined in a similar way as $`D_f`$ and $`D_{ff}`$ (see eqs. (10), (11)) but with the a new function $`\stackrel{~}{D}_{}(p;y,y^{})`$ substituted for $`D_{}(p;y,y^{})`$. The function $`\stackrel{~}{D}_{}(p;y,y^{})`$ is a solution of the following equation

$$𝒪(\mathrm{}^{(4+N)})(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}\stackrel{~}{D}_{}(p;y,y^{})=\delta (yy^{}),$$

where the operator $`𝒪(\mathrm{}^{(4+N)})`$ is

$$𝒪(\mathrm{}^{(4+N)})=\frac{N1+N(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}}{N+2(N+3)(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}}$$

with

$$(\mathrm{}^{(4+N)})=\frac{2𝒢(\mathrm{}^{(4+N)})}{(\mathrm{}^{(4+N)})}.$$

The “mass” $`m^2(p)`$ entering eq. (48) is the same as in Section 3, while $`m_{}^2(p)`$ is now

$$m_{}^2(p)=\frac{1}{2M_{Pl}^2\stackrel{~}{D}_{ff}(p^2)}$$

Thus, the brane-to-brane propagator still has the form (44). The second term in eq. (44) again has negative sign, so the model again has ghost field, but now the mass of the ghost is a solution to the following eigenvalue equation

$$m_{}^2=\frac{1}{2M_{Pl}^2\stackrel{~}{D}_{ff}(\stackrel{~}{m_{}^2)}}.$$
(49)
The difference with the case $`𝒢=0`$ studied in Section 3 is that the ghost field in principle need not be a tachyon in the general case and that the wave functions of the ghost and graviton have different profiles in the transverse directions. The graviton wave function is again $`D_f(p;y)`$, while the ghost wave function is $`\stackrel{~}{D}_f(p;y)`$.

## Acknowledgements

The authors are indebted to I. Antoniadis, A. Barvinsky, R. Rattazzi, S. Sibiryakov, M. Shaposhnikov, M. Shifman, P. Tinyakov and F. Viniegra for helpful discussions. We are grateful to G. Gabadadze for stimulating correspondence. This work has been supported in part by Russian Foundation for Basic Research, grant 02-02-17398, Swiss Science Foundation grant 7SUPJ062239, and the CRDF grant RP1-2364-MO-02. The work of S.D. has been supported in part by the INTAS grant YS 2001-2/128.

## Appendix A

Here we calculate the tachyon pole term in the full propagator of the model of Section 3. We do this by solving eq. (20) with conserved right hand side,

$$_AT_B^A=0$$
(50)
Otherwise $`T_{AB}`$ is arbitrary. We still use the gauge conditions (21).

We begin with $`(ab)`$-components of eq. (20) which read

$$\frac{1}{2}\mathrm{}^{(4+N)}(h_{ab}\frac{1}{2}\eta _{ab}h_D^D)=T_{ab}$$
(51)
We decompose $`h_{ab}`$ in the following way,

$$h_{ab}=h_{ab}^T+\frac{1}{N}\eta _{ab}h_c^c$$
(52)
where

$$h_{ab}^T=h_{ab}\frac{1}{N}\eta _{ab}h_c^c$$
(53)
is the traceless part. The traceless part obeys

$$\frac{1}{2}(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}h_{ab}^T=T_{ab}\frac{1}{N}\eta _{ab}T_c^c$$
(54)
while the trace of eq. (51) gives

$$h_c^c=\frac{N}{N2}h_\mu ^\mu +b$$
(55)
where $`b`$ obeys

$$(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}b=\frac{4}{N2}T_c^c$$
(56)
and hence is equal to

$$b(X)=\frac{4}{N2}d^{4+N}X^{}D_{}(XX^{})T_c^c(X^{})$$
(57)
We also have for the overall trace

$$h_C^C=\frac{2}{N2}h_\mu ^\mu +b$$
(58)
Let us now consider $`(a\mu )`$-components of eq. (20). They read

$$\frac{1}{2}(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}h_{a\mu }=T_{a\mu }$$
(59)
Hence,

$$h_{a\mu }(X)=2d^{4+N}X^{}D_{}(XX^{})T_{a\mu }(X^{})$$
(60)
Finally, let us study $`(\mu \nu )`$-components of eq. (20). We need the expression for $`_\lambda h^{\lambda \rho }`$ which enters $`G_{\mu \nu }^{(4)}`$. This expression is obtained by making use of the gauge condition

$$_Ah_\mu ^A_ah_\mu ^a+_\lambda h_\mu ^\lambda =\frac{1}{2}_\mu h_C^C$$
(61)
Now, because of the conservation property (50), one has

$$_ah_\mu ^a=2d^{4+N}X^{}D_{}(XX^{})_aT_\mu ^a(X^{})=2d^{4+N}X^{}D_{}(XX^{})_\lambda T_\mu ^\lambda (X^{})$$
(62)
Hence, $`(\mu \nu )`$-components of eq. (20) may be written in terms of $`h_{\mu \nu }`$, the trace $`h_\lambda ^\lambda `$ and components $`T_\lambda ^\rho `$ and $`T_a^a`$ of the stress-energy tensor. After some algebra one obtains

$``$ $`(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}\left(h_{\mu \nu }+{\displaystyle \frac{1}{N2}}\eta _{\mu \nu }h_\lambda ^\lambda \right)`$ (63)
$`+`$ $`M_{Pl}^2f^2(y){\displaystyle 𝑑y^{}f^2(y^{})\left[\frac{N}{N2}_\mu _\nu h_\lambda ^\lambda \mathrm{}^{(4)}h_{\mu \nu }+\frac{N1}{N2}\eta _{\mu \nu }\mathrm{}^{(4)}h_\lambda ^\lambda \right]}`$
$`=`$ $`2T_{\mu \nu }{\displaystyle \frac{2}{N2}}\eta _{\mu \nu }T_a^a`$
$``$ $`2M_{Pl}^2f^2(y){\displaystyle }d^Ny^{}D_f(y^{})[_\lambda _\mu T_\nu ^\lambda +_\lambda _\nu T_\mu ^\lambda {\displaystyle \frac{2}{N2}}_\mu _\nu T_a^a`$
$``$ $`\eta _{\mu \nu }(_\lambda _\rho T^{\lambda \rho }{\displaystyle \frac{1}{N2}}\mathrm{}^{(4)}T_a^a)]`$
Trace of this equation gives

$``$ $`(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}h_\mu ^\mu +\stackrel{~}{M}_{Pl}^2f^2(y){\displaystyle 𝑑y^{}f^2\mathrm{}^{(4)}h_\mu ^\mu }`$ (64)
$`=`$ $`{\displaystyle \frac{2(N2)}{N+2}}[T_\mu ^\mu {\displaystyle \frac{4}{N2}}T_a^a`$
$`+`$ $`2M_{Pl}^2f^2(y){\displaystyle }d^Ny^{}D_f(y^{})(_\lambda _\rho T^{\lambda \rho }{\displaystyle \frac{1}{N2}}\mathrm{}^{(4)}T_a^a)]`$
The solution to the latter equation is conveniently written in mixed representation,

$`h_\mu ^\mu (p,y)`$ $`=`$ $`{\displaystyle \frac{2(N2)}{N+2}}{\displaystyle d^Ny^{}D_{}(p;y,y^{})\left(T_\mu ^\mu (p,y^{})\frac{4}{N2}T_a^a(p,y^{})\right)}`$ (65)
$`+`$ $`{\displaystyle \frac{4(N2)}{N+2}}{\displaystyle \frac{p^2M_{Pl}^2D_f(y)}{1\stackrel{~}{M}_{Pl}^2p^2D_{ff}(p)}}`$
$`\times `$ $`{\displaystyle d^Ny^{}D_f(p,y^{})\left(\frac{N1}{N+2}T_\mu ^\mu \frac{p_\mu p_\nu }{p^2}T^{\mu \nu }\frac{3}{N+2}T_a^a\right)}`$
where we made use of the relation (29). This expression is still exact. Clearly, it has a pole at $`p^2=m_{}^2`$ (in the second term in the right hand side of eq. (65)).

To find the tachyon pole term in $`h_{\mu \nu }`$, one plugs the solution for $`h_\lambda ^\lambda `$ back into eq. (63). In terms of eq. (63), the tachyon pole in $`h_{\mu \nu }`$ comes entirely from the pole part in $`h_\lambda ^\lambda `$. One makes use of eq. (64) and writes eq. (63) in the following form

$``$ $`(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}h_{\mu \nu }+M_{Pl}^2f^2(y)p^2{\displaystyle 𝑑y^{}f^2(y^{})h_{\mu \nu }(p,y^{})}`$ (66)
$`=`$ $`M_{Pl}^2p^2f^2(y){\displaystyle d^Ny^{}f^2\left(\frac{N(N1)}{(N2)(N+2)}\eta _{\mu \nu }\frac{N}{N2}\frac{p_\mu p_\nu }{p^2}\right)h_\lambda ^\lambda }+\mathrm{}`$
where dots denote terms that do not contain a pole at $`p^2=m_{}^2`$. This equation is straightforwardly solved, and after some algebra one obtains that the tachyon pole part of $`h_{\mu \nu }`$ is

$`h_{\mu \nu }^{(pole)}(p,y)`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \frac{1}{p^2m_{}^2}}{\displaystyle \frac{D_f(p,y)}{M_{Pl}^2D_{ff}^2}}\left(\eta _{\mu \nu }{\displaystyle \frac{N+2}{N1}}{\displaystyle \frac{p_\mu p_\nu }{m_{}^2}}\right)`$ (67)
$`\times `$ $`{\displaystyle d^Ny^{}D_f(p,y^{})\left[\left(\eta _{\lambda \rho }\frac{N+2}{N1}\frac{p_\lambda p_\rho }{m_{}^2}\right)T^{\lambda \rho }(p,y^{})\frac{3}{N1}T_a^a(p,y^{})\right]}`$
It remains to find the tachyon pole parts of other metric components. According to eqs. (54), (56) and (59), the traceless part $`h_{ab}^T`$, the term $`b`$ and the metric components $`h_{a\mu }`$ do not have poles at $`p^2=m_{}^2`$. The pole term in $`h_a^a`$ is determined by the pole term in $`h_\mu ^\mu `$ through eq. (55). Thus, one finds

$`h_{ab}^{(pole)}(p,y)`$ $`=`$ $`{\displaystyle \frac{1}{N2}}\eta _{ab}h_\mu ^{(pole)\mu }(p,y)`$ (68)
$`=`$ $`{\displaystyle \frac{1}{N1}}\eta _{ab}{\displaystyle \frac{1}{p^2m_{}^2}}{\displaystyle \frac{D_f(p,y)}{M_{Pl}^2D_{ff}^2}}`$
$`\times `$ $`{\displaystyle d^Ny^{}D_f(p,y^{})\left[\left(\eta _{\lambda \rho }\frac{N+2}{N1}\frac{p_\lambda p_\rho }{m_{}^2}\right)T^{\lambda \rho }(p,y^{})\frac{3}{N1}T_a^a(p,y^{})\right]}`$
We see from eqs. (67) and (68) that the pole terms in $`h_{AB}`$ may indeed be written in the form

$$h_{AB}^{(pole)}(p,y)=d^Ny^{}D_{AB,CD}^{(pole)}(p;y,y^{})T^{CD}(p,y^{})$$
(69)
where the pole term in the propagator is given by eq. (39).

## Appendix B

Let us calculate the $`(\mu \nu )`$-components of the metric due to the source of the form (40) with conserved $`\theta _{\mu \nu }`$. For this particular type of source, the expression (65) simplifies considerably,

$$h_\mu ^\mu (p,y)=\frac{2(N2)}{N+2}\frac{D_f(p,y)}{\stackrel{~}{M}_{Pl}^2D_{ff}(p)}\frac{1}{p^2m_{}^2(p)}\theta _\mu ^\mu (p)$$
(70)
where $`m_{}^2(p)`$ is given by eq. (43). We plug this expression into eq. (63) and again make use of the properties of the source to obtain

$``$ $`(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}h_{\mu \nu }+M_{Pl}^2f^2(y)p^2{\displaystyle 𝑑y^{}f^2(y^{})h_{\mu \nu }(p,y^{})}`$ (71)
$`=`$ $`f^2(y)\left(2\theta _{\mu \nu }{\displaystyle \frac{2}{N+2}}\eta _{\mu \nu }\theta _\lambda ^\lambda {\displaystyle \frac{N}{N+2}}{\displaystyle \frac{p^2}{p^2m_{}^2(p)}}\eta _{\mu \nu }\theta _\lambda ^\lambda \right)`$
$`+`$ longitudinal part
where the longitudinal part is proportional to $`p_\mu p_\nu `$. After some algebra, one finds that the solution has indeed the form (41).

## Appendix C

Here we sketch the calculations leading to the result (48). The steps are similar to those in Appendices A and B. First, one considers $`(ab)`$-components of eq. (47) and finds the following generalization of eq. (55) (recall that we consider a source of the form (40)),

$$h_c^c=\frac{N+\left(\left(N1\right)\mathrm{}^{(4+N)}+\mathrm{}^{(4)}\right)}{N2+\left(\left(N1\right)\mathrm{}^{(4+N)}+\mathrm{}^{(4)}\right)}h_\mu ^\mu .$$
(72)
For the overall trace one has

$$h_C^C=\frac{2}{N2+\left(\left(N1\right)\mathrm{}^{(4+N)}+\mathrm{}^{(4)}\right)}h_\mu ^\mu .$$
(73)
From the $`(a\mu )`$-components of eq. (47) one finds

$$h_{a\mu }=\frac{}{N2+\left(\left(N1\right)\mathrm{}^{(4+N)}+\mathrm{}^{(4)}\right)}_a_\mu h_\nu ^\nu .$$
(74)
By making use of the gauge condition (61) one obtains the following expression for the longitudinal components of the four-dimensional part of the metric,

$$_\nu h_\mu ^\nu =\frac{1\mathrm{}^{(N)}}{N2+\left(\left(N1\right)\mathrm{}^{(4+N)}+\mathrm{}^{(4)}\right)}_\mu h_\nu ^\nu ,$$
(75)
where $`\mathrm{}^{(N)}=\delta ^{ab}_a_b`$. Plugging expressions (73), (75) into the $`(\mu \nu )`$-components of eq. (47) one arrives at the following analog of eq. (63),

$``$ $`(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}\left(h_{\mu \nu }+{\displaystyle \frac{1+\mathrm{}^{(N+4)}}{N2+\left(\left(N1\right)\mathrm{}^{(4+N)}+\mathrm{}^{(4)}\right)}}\eta _{\mu \nu }h_\lambda ^\lambda \right)`$ (76)
$`+`$ $`(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}{\displaystyle \frac{}{N2+\left(\left(N1\right)\mathrm{}^{(4+N)}+\mathrm{}^{(4)}\right)}}_\mu _\nu h_\lambda ^\lambda `$
$`+`$ $`M_{Pl}^2f^2(y){\displaystyle }dy^{}f^2(y^{})[{\displaystyle \frac{N+\left(\left(N+1\right)\mathrm{}^{(4+N)}\mathrm{}^{(4)}\right)}{N2+\left(\left(N1\right)\mathrm{}^{(4+N)}+\mathrm{}^{(4)}\right)}}_\mu _\nu h_\lambda ^\lambda `$
$``$ $`\mathrm{}^{(4)}h_{\mu \nu }+{\displaystyle \frac{N1+N\mathrm{}^{(4+N)}}{N2+\left(\left(N1\right)\mathrm{}^{(4+N)}+\mathrm{}^{(4)}\right)}}\eta _{\mu \nu }\mathrm{}^{(4)}h_\lambda ^\lambda ]=2\theta _{\mu \nu }f^2(y).`$
Trace of this equation gives

$``$ $`(\mathrm{}^{(4+N)})\mathrm{}^{(4+N)}{\displaystyle \frac{N+2+\left(N+3\right)\mathrm{}^{(4+N)}}{N2+\left(\left(N1\right)\mathrm{}^{(4+N)}+\mathrm{}^{(4)}\right)}}h_\mu ^\mu `$ (77)
$`+`$ $`2M_{Pl}^2f^2(y){\displaystyle 𝑑y^{}f^2\mathrm{}^{(4)}\frac{N1+N\mathrm{}^{(4+N)}}{N2+\left(\left(N1\right)\mathrm{}^{(4+N)}+\mathrm{}^{(4)}\right)}h_\mu ^\mu }=2\theta _\mu ^\mu f^2(y).`$
The solution of this equation is

$$h_\mu ^\mu =\frac{2\theta _\mu ^\mu }{12p^2D_{ff}(p)}\frac{N2+(\left(N1\right)\mathrm{}^{(4+N)}+\mathrm{}^{(4)})}{N+2+\left(N+3\right)\mathrm{}^{(4+N)}}\stackrel{~}{D}_f.$$
(78)
Plugging this result back into eq. (76) one obtains after some algebra the desired expression (48).


format: text

text:
# Transport on percolation clusters with power-law distributed bond strengths: when do blobs matter?

## I Introduction

The transport properties of percolation clusters have been a subject of interest for many years SA94 ; SN98 . A natural problem to study is e.g. conductivity, and one often complicates it further by using random bond “strengths” $`\sigma `$ with a power-law tail of the form $`\sigma ^\alpha `$ KSD79 ; SN82 ; HFSD85 ; MGMC86 ; FHST87 ; MC88 ; LTR88 ; OLN91 ; SJC01 . In the first place, because this allows to represent continuum percolation HFSD85 ; FHST87  and thus get closer to some actual physical realizations of percolation. A second reason why these systems are interesting is the equivalence AHLH71  between transport on strongly disordered systems and percolative transport 111 Although this equivalence is strict only in the $`\alpha 1`$ limit, it holds with corrections for $`1>\alpha `$ as well..
Transport critical exponents on these systems are found to depend on $`\alpha `$, which means that strict universality is lost. The original observation that transport exponents become nonuniversal is due to Kogut and Straley KSD79 , who used mean-field type arguments. Later Straley SN82 , with the help of the nodes-links-blobs SST75 ; DR76 ; SC77 ; CC82  picture of the backbone, concluded that the conductivity exponent $`t`$, such that $`\mathrm{\Sigma }(pp_c)^t`$, is the maximum of the universal exponent $`t_0=(d2)\nu +\zeta `$ and the $`\alpha `$-dependent exponent $`\overline{t}(\alpha )=(d2)\nu +1/(1\alpha )`$. Here $`\nu `$ is the correlation length exponent and $`\zeta `$ measures the contribution of blobs to the resistance between two points on the backbone, for the case of constant conductances. For the conductivity problem thus there is a crossover from the universal exponent $`t_0`$ for $`\alpha <\alpha _c`$ to $`\overline{t}(\alpha )`$ in the “anomalous regime” $`\alpha >\alpha _c`$. Although not without some controversy initially BBN81 ; MGMC86 ; MC88 ; LTR88 , this result is by now well established MC88 ; LTR88 ; OLN91 ; SJC01 .

It is somehow surprising that $`\overline{t}(\alpha )`$ can be analytically calculated in the “anomalous regime”, given that the universal exponent $`t_0`$, which applies to the arguably simpler case of constant conductance, has not been analytically derived up to now. The difficulty in deriving $`t_0`$ resides in that $`\zeta `$ is determined by the blobs, and thus one would require detailed information BHF96  about the internal structure of the blobs.

On the contrary, it has been argued by several authors MGMC86 ; BHF96 ; SJC01  that $`\overline{t}(\alpha )`$ in the anomalous regime is determined by the cutting-bonds alone. Since these form linear chains of typically $`L^{1/\nu }`$ bonds at $`p_c`$ CC82 , the resulting conductivity exponent is easily derived. The argument to support the belief that blobs are irrelevant in the anomalous regime seems to be roughly the following: an exceedingly small conductivity falling on a blob has little effect on the overall conductance, because there are many alternative parallel paths. On the other hand, if this small conductivity is located on a cutting bond it will certainly dominate the system conductance. While this argument is true in principle, this reasoning misses the fact that the number of cutting bonds is itself a fluctuating quantity. The issue of blob-irrelevance has been considered by Machta and collaborators MGMC86  using a hierarchical model for the backbone, to reach similar conclusions. However, as noted by the authors, their model does not include structural fluctuations. We will demonstrate in this work that it is in fact the blobs and not the cutting bonds, that determine the critical transport properties, even in the anomalous regime. This as we will see, is due to structural fluctuations. However, the resulting transport properties turn out to be the same as those given by the most simplistic red bond estimate.

A related critical transport problem, which is relevant for disordered superconductors, is that of determining the critical current density $`J_c=I_cL^{(d1)}`$ that a percolation network can sustain, and which above $`p_c`$ behaves as $`J_c(pp_c)^v`$ DRC79 ; LHSN87 ; OOAN88 ; GGWC88 ; HOAN94 ; MDHN99 . This problem has a simple geometrical interpretation. Finding the maximum flow of current, or *maxflow*, is equivalent to finding the surface across the system, on which the sum of critical currents of the bonds is maximized. As we will draw advantage of this analogy later on, we note that this surface is called a *mincut* in computer science language PSC98 ; CLRI93 ; RBL89 .

In this paper our goal is to present the first comprehensive study of the maxflow problem on percolation clusters. This is motivated by the following observations. First, this is the simplest transport problem that one can think about, and has not been as such discussed much in the literature. Second, we are able to use to our advantage recent developments ADME01  on combinatorial optimization algorithms, in the context of disordered systems. Here one can use a three-step approach, in which first a critical spanning cluster is set, its backbone is pruned out, and finally that is used for the max-flow/min-cut problem. Each stage is solved with one of the powerful graph optimization algorithms for the particular problem, as discussed later.

In the simplest version of the maxflow problem, all present bonds have the same critical current, or *capacity* $`i_c`$ and absent bonds have $`i_c=0`$. At criticality, a typical percolating cluster is a linear chain of cutting bonds and thus $`I_c=i_c`$. From this observation plus the usual scaling relation $`J_c(L)L^{v/\nu }`$, one concludes that DRC79  $`v=(d1)\nu `$. This result is consistent with experiments DRC79 ; GGWC88  and numerical simulation OOAN88 ; MDHN99 . In a more realistic model, each present bond has a random capacity $`i_c`$ with power-law distribution $`P(i_c)i_c^\alpha `$. This is for example the case for continuum percolation models LHSN87 ; OOAN88 ; HOAN94 . A simple extension of the “typical cutting bond string” argument gives $`v(\alpha )=(d1)\nu +1/(1\alpha )`$ as we show later in Section II.

In the following we will find it useful to compare the conductivity and critical current problems to each other. This comparison is done by interpreting the random bond variables $`i_c`$ alternatively as bond conductances $`\sigma `$ or as bond capacities $`i_c`$. Consider for example two bonds with $`\sigma _1`$ and $`\sigma _2`$ connected in parallel. The resulting conductance $`\sigma _{par}=\sigma _1+\sigma _2`$ is then the same as the maximum current $`I_{max}`$ that can flow if $`\sigma _j`$ are capacities. If these bonds are instead connected in series, then $`\sigma _{series}=(\sigma _1^1+\sigma _2^1)^1`$ and $`I_{max}=min(\sigma _1,\sigma _2)`$ are no longer equal. However the series conductance can be written as LP89  $`\sigma _{series}=min(\sigma _1,\sigma _2)(1+\beta )^1`$, with $`\beta =min(\sigma _1,\sigma _2)/max(\sigma _1,\sigma _2)`$. In the limit of strong disorder ($`\alpha 1`$), $`\beta `$ is typically negligible. We conclude that, in this limit, also in the series case the conductance equals exactly the maximum current obtained by interpreting $`\sigma _j`$ as capacities $`i_c`$. Therefore in the $`\alpha 1`$ limit, the resistive current problem and the superconducting current problem (maxflow) are equivalent, at least for all structures that can be solved by a combination of series and parallel bond reductions 222Similar ideas were put forward in MGMC86 , for hierarchical structures which are reducible by series and parallel transformations.. Moreover as shown in Section III.2, we find that the equivalence noticed above is valid not only in the $`\alpha 1`$ limit but for a range of $`\alpha `$ values, for strings of bonds in series.

In deriving $`\alpha `$-dependent exponents, both for $`v(\alpha )`$ and $`\overline{t}(\alpha )`$, the assumption is made that the backbone always contains $`L^{1/\nu }`$ cutting bonds. While this is true typically, the number of cutting bonds is in fact a fluctuating variable whose distribution may extend down to zero in the form of a power law (see later). The existence of such fluctuations has been noted by some works previously KS86b ; HAR97 , but their role in transport properties has not been considered. These number fluctuations, we will show in Section II, do modify the transport exponent that results from a string of cutting bonds. Then by analyzing the conceptually and numerically simple maxflow problem, we will be able to show that in fact blobs cannot be neglected. The net outcome, which we justify by a heuristic hierarchical picture, is that although the simplest cutting bond scaling (without fluctuations) is restored, it is in fact the blobs that set this scaling behavior.

The structure of the rest of the paper is as follows. Section II presents the analytical discussion, based on a “fluctuating number of cutting-bonds” picture. In Section III we go through one-by-one the numerical methods employed, the findings about structural fluctuations, and some further numerical analysis of the extremal statistics aspects. Section IV contains the results concerning the maxflow problem, and some details of interest that can be determined from analyzing large statistics. Section V finishes the paper with a discussion.

## II Critical current density

We consider diluted lattices where the maximum supercurrent $`i_c`$ that a present bond can sustain is a random variable distributed between $`0`$ and $`1`$ according to

$$P(i_c)=(1\alpha )c^\alpha ,$$
(1)
with $`\alpha <1`$.

Let $`I_c`$ be the maximum supercurrent (or *maxflow*) that the whole system, given a set of values $`\{i_c\}`$, can sustain. The average current density $`J_c`$ is then $`J_c=<I_c>/L^{(d1)}`$, and goes to zero at $`p_c`$ as

$$J_c(pp_c)^v$$
(2)
Right at $`p_c`$, and for a system of finite linear size $`L`$, usual finite-size scaling arguments SAI94 ; BHF96  imply that

$$J_c(p_c,L)L^{v/\nu },$$
(3)
where $`\nu `$ is the percolation correlation length exponent. The nodes-links-blobs picture of the percolation cluster SST75 ; DR76 ; SC77 ; CC82  tells us that, right at $`p_c`$, there is typically a single connected path through the sample. This path is a sequence of multiply connected regions (blobs) connected by strings of singly connected bonds, also called *cutting bonds*. The average number of cutting bonds is of order $`L^{1/\nu }`$ at $`p_c`$ CC82 .

We now start by considering the maximum flow $`f^{}`$ allowed by a string of $`n`$ cutting bonds, and which obviously equals the least capacity among the $`n`$ bonds. The typical least value $`f_n^{}`$ among a collection of $`n>>1`$ random numbers $`i_c`$ with probability $`P(i_c)`$ satisfies

$$_0^{f_n^{}}P(i)𝑑i=1/n$$
(4)
Thus

$$f_n^{}=n^{1/(1\alpha )}$$
(5)
On a system of linear size $`L`$ at $`p_c`$, the average number of cutting bonds is $`L^{1/\nu }`$ CC82 . In replacing this one obtains $`f_n^{}L^{1/\nu (1\alpha )}L^{d1}J_cL^{(d1)v(\alpha )/\nu }`$ and thus from Eq. (3),

$$v(\alpha )=(d1)\nu +1/(1\alpha )$$
(6)
as advanced in the introduction.

This typical-$`n`$ argument however neglects the fact that $`n`$ is a fluctuating number. Since $`𝒫(n)`$ actually has a power law tail extending down to $`n=0`$ HAR97 , this neglect turns out not to be correct for quantities that depend on $`1/n`$ as Eq. (5).

We now present a more careful treatment, which takes into account the fluctuations in $`n`$. It is known HAR97  that $`𝒫_L(n)=(n_L^{})^1\widehat{𝒫}(n/n_L^{})`$, where $`\widehat{𝒫}(\widehat{n})`$ is a size-independent function, and $`n_L^{}L^{1/\nu }`$ CC82 . Since for the purpose of our discussion all that matters is the behavior of $`𝒫(n)`$ as $`\widehat{n}0`$, we take for simplicity $`\widehat{𝒫}(\widehat{n})=(1+a)\widehat{n}^a`$, for $`0<\widehat{n}1`$. Thus,

$$𝒫(n)=(1+a)(n_L^{})^{(1+a)}n^a,$$
(7)
for $`1nn_L^{}`$. We will for the moment assume that $`n`$ cannot be zero.

Let now $`f`$ be the minimum among $`n`$ numbers $`x`$ distributed with probability $`P(x)`$. The distribution $`m_n(f)`$ of $`f`$ is determined as

$`m_n(f)`$ $`=`$ $`nP(f)\left\{1{\displaystyle _0^f}P(x)𝑑x\right\}^{n1}`$ (8)
$``$ $`nP(f)e^{n_0^fP(x)𝑑x}.`$
Because of the strong exponential suppression that occurs for $`f`$ larger than $`f_n^{}`$ defined by Eq. (4), we can approximate $`m_n(f)`$ by

$$m_n(f)\{\begin{array}{cccc}& nP(f)\hfill & \text{if}\hfill & 0<ff_n^{}\hfill \\ & & & \\ & 0\hfill & \text{if}\hfill & f>f_n^{}\hfill \end{array}$$
(9)
Now allowing for the fact that $`n`$ fluctuates, the probability distribution function (PDF) of the maxflow $`f`$ through a string of cutting bonds is

$`m(f)`$ $`=`$ $`{\displaystyle _1^{\mathrm{}}}𝑑n𝒫(n)m_n(f)`$ (10)
$`=`$ $`(1+a)(n_L^{})^{(1+a)}{\displaystyle _1^{n_L^{}}}𝑑nn^am_n(f),`$
for $`0<f<1`$. From Eq. (5) and Eq. (9) we conclude that, for a given value of $`f`$, the only nonzero contributions in Eq. (10) come from $`n`$ values which are smaller than $`\eta (f)=f^{(1\alpha )}`$. Thus

$`m(f)`$ $`=`$ $`{\displaystyle \frac{P(f)(1+a)}{(n_L^{})^{1+a}}}{\displaystyle _1^{\text{min}(n_L^{},\eta (f))}}n^{a+1}𝑑n`$ (11)
Defining $`f_{typ}=(n_L^{})^{1/(1\alpha )}`$, $`\kappa =(a+1)(1\alpha )`$ and $`\lambda =\kappa /(a+2)`$, this last expression can be written as

$$m(f)=\{\begin{array}{cccc}& \lambda f_{typ}^1\left(\frac{f}{f_{typ}}\right)^\alpha \hfill & \text{if}\hfill & 0<f<f_{typ}\hfill \\ & & & \\ & \lambda f_{typ}^1\left(\frac{f}{f_{typ}}\right)^{(1+\kappa )}\hfill & \text{if}\hfill & f_{typ}<f<1.\hfill \end{array}$$
(12)
This gives the PDF for the maxflow $`f`$ through a string of cutting bonds on a system of size $`L`$, allowing for fluctuations in the number $`n`$ of bonds on the string. The strength of the fluctuations of $`1/n`$ is characterized by the exponent $`\kappa `$, which in turn depends on $`a`$. If $`a\mathrm{}`$ (non-fluctuating limit) $`m(f)`$ is nonzero only for $`f<f_{typ}`$. Thus $`<f>f_{typ}`$ and Eq. (6) is recovered in this case. However it is known that $`a0.22`$ in two dimensions HAR97 .

For general $`a`$ and $`\alpha `$, $`m(f)`$ has a power-law tail with exponent $`(1+\kappa )`$ for $`f>>f_{typ}`$. The importance of this power-law tail is evidenced by considering the average flow

$`<f>`$ $`=`$ $`{\displaystyle \frac{\lambda }{2\alpha }}f_{typ}+{\displaystyle \frac{\lambda }{1\kappa }}\left(f_{typ}\right)^\kappa `$ (13)
When $`\kappa >1`$, $`<f>f_{typ}(n_L^{})^{1/(1\alpha )}L^{1/\nu (1\alpha )}`$ and Eq. (6) is recovered. However if $`\alpha >a/(a+1)`$ ($`\kappa <1`$), the power-law tail dominates the average. In this case $`<f>f_{typ}^\kappa >>f_{typ}`$. Therefore $`<f>L^{(a+1)/\nu }`$, and Eq. (3) implies that in this case,

$$v=(d1)\nu +a+1$$
(14)
The meaning of this is clear. If $`\alpha `$ is large, typical cases with $`𝒪(n_L^{})`$ cutting bonds will only allow an exceedingly small flux $`f`$. The average flow $`<f>`$ however will be dominated by the very rare cases in which $`n`$ is small and for which $`f𝒪(1)>>f_{typ}`$. So finally we conclude that, if we idealize the backbone at $`p_c`$ as a string of $`n`$ cutting bonds, and if $`𝒫(n)`$ behaves for small $`n`$ as $`n^a`$, one has that

$$v(\alpha )=\{\begin{array}{cccccc}& \frac{1}{1\alpha }+\nu (d1)\hfill & \text{if}\hfill & \alpha <\frac{a}{1+a}\hfill & & \\ & & & & & \\ & a+1+\nu (d1)\hfill & \text{if}\hfill & \alpha >\frac{a}{1+a}\hfill & & \end{array}$$
(15)
## III Numerical Results

### III.1 Algorithms

In this section we test our analytical derivation of $`v(\alpha )`$ of Section II in two and three dimensions on large systems, with the help of powerful combinatorial algorithms ADME01 . Percolation backbones are first generated by means of a matching algorithm MA98 ; ADME01 , for square and cubic lattices. We do this by randomly adding bonds one at a time until a percolation path is first present. At this point the matching algorithm identifies the conducting backbone, exactly at the percolation point for each sample. Alternatively one could fix the density of present bonds to a value close to the infinite system critical density $`p_c`$, and then identify the percolating backbone with the same algorithm. However, our procedure has the advantage that no separate estimate is necessary for $`p_c`$.

For each percolating backbone, capacities are drawn from the given distribution, and the maxflow is calculated by means of a flow augmentation algorithm (see ADME01  for a review of the max-flow problem). The efficiency of the maxflow algorithm is highly increased when working on the backbone only, so we are able to analyze thousands of samples for each value of $`\alpha `$. In this way we estimate numerically the average flow at $`p_c`$ for several linear sizes $`L`$, and from its scaling properties $`v(\alpha )`$ is derived. The largest sample sizes studied were $`L=4000`$ in two dimensions and $`L=120`$ in three dimensions. These are mostly set by the CPU usage of the combination of the matching and flow algorithms, which in turn is dominated for $`L`$ large by the scaling of the matching part. The maxflow code is actually sub-linear in $`n=L^2`$ in CPU time, since the mass of the backbone scales with its fractal dimension. Notice that once the backbone of a sample has been established, it can be used for several consequent maxflow determinations for different $`\alpha `$ to save CPU-time. In the Appendix we present an idea for an optimal algorithm for this problem.

### III.2 Results

Results are shown in Fig. 1.

Our numerical simulation results confirm Eq. (6) nicely. However the saturation of $`v(\alpha )`$ predicted by Eq. (15) for $`\alpha >a/(a+1)`$ does not occur. Notice that $`a`$ is not a universal exponent but depends on the ensemble.

For example if the ensemble is determined by fixing $`p=p_c`$, numerical measurements and RG calculations HAR97  give $`a0.22`$. Additionally, roughly $`20\%`$ of the connected samples have zero cutting bonds KS86b ; HAR97 , that is, $`𝒫(n)0.20\delta (n)+cn^{0.22}`$ for small $`n`$.

However other ensembles can be considered. Consider for example the percolation cluster defined by Ambegaokar *et al* construction, in which conductances are laid down on the lattice in order of increasing conductivity until a percolating path is created AHLH71 . At least the last conductance to be laid down is a cutting bond, so $`𝒫(0)=0`$. Experimentally this situation is realized when superconductive samples grow percolatively by deposition HOAN94 . In this case the point at which the supercurrent is nonzero for the first time is defined by the first appearance of a connected path, not by a fixed density of occupied bonds.

As we add bonds one at a time, our numerical simulations correspond to this case rather than to fixing $`p=p_c`$. Our measurements of the distribution of the number of cutting bonds indicate (Fig. 2) that $`𝒫(n)n^a`$ for small $`n`$, with $`a1.25`$ in two and three dimensions.

Different ensembles give rise to different distributions of the number of cutting bonds, and specifically to different values of $`a`$ so, if Eq. (15) were to hold for percolation clusters, the resulting transport exponent would be ensemble-dependent. However, our maxflow measurements on percolation clusters are consistent with Eq. (6) for all $`\alpha `$, without any sign of saturation.

In view of the failure of percolation clusters to show the predicted exponent saturation, we first confirmed the validity of Eq. (15) for strings of cutting-bonds. We did so by numerically studying strings of bonds whose number $`n`$ is distributed according to (7), and whose conductances (or capacities, for the maxflow problem) are distributed according to (1). The maximum flow is simply the least critical current $`i_c`$. Alternatively, bond capacities $`i_c`$ may be interpreted as conductances, in which case the resulting conductance for the whole string is simply $`\sigma =1/_{j=1}^n1/i_c(j)`$.

We find for these strings of cutting bonds (Fig. 3) that Eq. (15) is satisfied very accurately. Fig. 3 also shows that the conductivity and maxflow exponents are the same for $`\alpha >0`$, indicating that the conductance is dominated by the least $`i_c`$ value in that regime. We conclude that Eq. (15) is exact for strings of cutting-bonds. Thus the failure of Eq. (15) for percolation clusters simply means that these *do not* behave as strings of cutting bonds do. In other words, for $`\alpha `$ near $`1`$, it is not correct to approximate a percolation cluster as a string of cutting bonds.

## IV The role of blobs

### IV.1 Structural fluctuations

Our results (Fig. 1) show that the maxflow exponent $`v(\alpha )`$ follows Eq. (6), although fluctuations in the number $`n`$ of cutting bonds, which exist and are relevant in real percolation clusters, were disregarded in its derivation. So we face a somehow paradoxical situation, since a naive calculation gives the correct result (Eq. (6)), while a seemingly more careful calculation that takes into account the fluctuations in $`n`$ (Eq. (15)) does not. As mentioned in the previous section, this means that our assumption that the maximum flow is determined by the cutting bonds alone needs to be revised. In order to test this assumption, we separately measure the maximum flow allowed by cutting bonds and by blobs, which we call $`m_c`$ and $`m_b`$ respectively, for each percolation cluster. The overall maximum flow is the minimum of these. The procedure works such that one picks first the smallest of the cutting bond capacities, and then assigns to them an infinite capacity. Then the maxflow is found, which is now given by the minimal blob min-cut (configuration).

Figures 4 a to c show how the PDF of $`m_c`$ (cutting-bond flow) and the total PDF vary with $`\alpha `$. For non-anomalous values $`\alpha =0`$ (Fig. 4a), the distribution is centered around a well-defined mean value. With increasing $`\alpha `$ one enters the anomalous regime, and the PDF develops a power-law tail. This would be expected to result from the cutting bonds, while the blob flows $`m_b`$ have a much narrower distribution, decaying roughly exponentially for large flows. This means that, when $`m_c`$ is large, most probably $`m_b`$ will be much smaller and thus the overall flow will be determined by $`m_b`$. Thus, although our derivation of Eq. (15) is correct for strings of cutting bonds, it is the blobs that determine the flow in those rare cases in which $`m_c`$ is large. Therefore the power-law tail in $`P(m_c)`$, which is responsible for the saturation of $`v(\alpha )`$ at large values of $`\alpha `$ in Eq. (15), is suppressed by blobs on percolation clusters. Figure 5 illustrates this by showing that the fraction of cases - for a given maxflow $`m`$ \- that are dominated by the blob contribution follow a separate PDF. The collapse is not completely perfect, since there may be a very slight trend in the total fraction of blob-dominated cases with increasing $`L`$. On the other hand, the variances of the maxflow distributions scale as expected (as the mean). It is worth mentioning that the distribution of $`k`$-cuts (number of bonds in the min-cut) is roughly exponential, so that $`k`$ is of the order of 1.4 $`\mathrm{}`$ 1.5 for $`\alpha =0.7`$.

It is also worth pointing out that there are cross-correlations between the structural quantities on one hand, and between the structure and the maxflow on the other hand. These are illustrated in Figures 6 and 7. In a system with a given $`L`$ it is after a moment’s deliberation rather clear that there may be an inverse correlation between the  *number* of cutting bonds and the sample-to-sample weight of the backbone. We have not tried to measure this relationship quantitatively, but given such a relation it is no surprise (Fig. 7) that the mass of the backbone correlates strongly with the maxflow value.

### IV.2 Blob dominance

In order to prove that our hypothesis, namely that the blobs set the maxflow scale, is correct, we still have to show that the blob flow $`m_b`$ has the right scaling properties, i.e. $`\overline{m}_bL^{1/\nu (1\alpha )}`$. A complete calculation of the maximum flow allowed by blobs would require detailed information about the blob’s internal structure. However an estimate can be obtained from the following arguments. It is known that the backbone at $`p_c`$ has a hierarchical, or self-similar, structure CC82 . At the top level of this hierarchy, the backbone itself can be thought of as a string of singly connected (cutting) bonds interspersed with blobs. Blobs in turn are loops made of doubly connected bonds interspersed with smaller blobs and so on, as depicted in Fig. 8. This hierarchical structure has its counterpart in a similar classification of surfaces which separate the backbone into two pieces (cuts). At the top level of this hierarchy are the surfaces $`\{S_1`$} that cut the backbone at just one bond, next come those surfaces $`\{S_2\}`$ that cut the backbone at exactly two bonds, etc. The capacity $`C(S)`$ of a cut $`S`$ is defined as the sum of the capacities $`i_c`$ of the bonds crossed by it. Because of the maxflow-mincut theorem, the maximum flow equals the minimum of the cuts’ capacities. Our assumption that cutting bonds alone determine the maximum flow is equivalent to minimizing the capacities among the $`S_1`$ alone. We now describe how the next level $`S_2`$ in this hierarchy can be analyzed.

Coniglio CC82  has shown that the derivative of the spanning probability $`p^{}`$ with respect to $`p`$ is proportional to the average number of cutting bonds $`<n>`$. An extension of his reasoning, due to Kantor KS86b , allows one to write the second derivative of $`p^{}(p)`$ with respect to $`p`$ at $`p_c`$ as $`^2p^{}/p^2|_{p_c}<n(n1)2N_2>|_{p_c}`$, where $`n`$ is the number of cutting bonds and $`N_2`$ is the number of pairs of doubly connected bonds. Because by definition $`^2p^{}/p^2=0`$ at $`p_c`$ 333Here we depart slightly from Kantor’s reasoning, who justifies this last point as due to duality. Duality is not necessary and therefore $`2<N_2>=<n^2><n>`$ for any lattice and dimension, at $`p_c`$., one finds that $`2<N_2>=<n^2><n>`$ at $`p_c`$. Since $`<n^k>L^{k/\nu }`$ 444 This results from the fact that the distribution of $`n`$ decays exponentially fast for large $`n`$. See Fig. 2 and ref. HAR97 . , we conclude that the typical number of pairs of doubly connected bonds at $`p_c`$ is $`N_2L^{2/\nu }`$. However this alone is not enough to estimate the typical maximum flow allowed by doubly connected bonds, for they might be grouped into blobs in different ways. Fortunately the total number $`n_2`$ of doubly connected bonds at $`p_c`$ can also be calculated CC82 , and it turns out to be $`<n_2>L^{1/\nu }`$. This means that the blob statistics is dominated by one large ring of roughly $`L^{1/\nu }`$ bonds and therefore containing a number of pairs of doubly connected bonds which is of order $`L^{2/\nu }`$. Using this information we can now estimate the maximum flow allowed by blobs at the level of doubly connected bonds. This large blob dominates the maximum flow since lesser blobs, located somewhere else along the backbone, will allow a larger flow. Thus one has to find the maximum flow for two parallel strings, each containing $`L^{1/\nu }`$ cutting bonds. The typical flow allowed by each string is of order $`L^{1/\nu (1\alpha )}`$ and therefore the typical maximum flow allowed by doubly-connected blobs, which is twice this, is of the right order.

Our reasoning for doubly connected bonds only considers typical cases, i.e. fluctuations in the number of doubly connected bonds are disregarded. If a particular cluster, in addition to having a small number of cutting bonds, also has a small number of doubly connected bonds, then the next levels in this hierarchy would be relevant. The same sort of reasoning can be used at all levels in the hierarchy of cuts, but because the algebra becomes too complicate for triply connected bonds already, we did not test this in detail. However it seems safe to assume that $`min\{C(S_k)\}L^{1/\nu (1\alpha )}`$ for all $`k>2`$ as well. Additionally notice that, in order for the mincut to be located at triply connected bonds, it is necessary that the numbers of singly and doubly connected bonds be simultaneously small, an occurrence which arguably has a small probability.

We then see that, in those rare cases in which the number of singly connected bonds is small (they allow a large flow), blobs take their role thus limiting the flow to a value which is typically of order $`L^{1/\nu (1\alpha )}`$. This then shows that the correct value of the exponent $`v(\alpha )`$ is given by Eq. (6). Similar behavior is of course expected for other transport properties, e.g. conductivity, in the limit of anomalous distributions of bond strengths. This is so because in this limit the resistance of the whole cluster is dominated by that of the mincut, where conductivities are interpreted as critical currents.

## V Conclusions

In this article we have demonstrated that the transport problem on percolation clusters still holds surprises. Our findings deny the widespread notion that, in the limit of anomalous strength distributions, it is the cutting-bonds alone that determine the transport properties. We show analytically, and confirm numerically that, if blobs could be neglected (because of their allowing a larger maxflow than cutting bonds) then the overall system’s behavior would be strongly dependent on the ensemble (the cutting-bond PDF tail exponent). This ensemble-dependence would come about because the number of cutting bonds has a “broad” distribution extending down to zero. However the predicted ensemble-dependence is not there, as we show numerically on large two and three-dimensional systems. Using scaling arguments we then demonstrate that it is in fact the blobs that finally determine the average maxflow. However, we are forced to finish with the paradoxical conclusion that though the expected mechanism for the maxflow, namely cutting-bond dominance, does not work in the anomalous regime (large $`\alpha `$), the original cutting-bond estimate for the transport exponent is nevertheless restored by the limiting effect of the blobs.

###### Acknowledgements.

This collaboration was partially supported by Väisälän Foundation, Finland. CFM acknowledges finantial support by CONACYT, México, through research project 36256-E. CFM also wishes to acknowledge the kind hospitality of the Laboratory of Physics, HUT, where parts of this work were done. MJA is supported by the Academy of Finland Center of Excellence program. Appendix - an optimal algorithm: We note that the augmenting path -method is better here than in the general maxflow problem so-called push/relabel preflow algorithms GTA88  enjoy the most popularity). This is since the structure of the backbone is essentially one-dimensional, the number of augmentations remains small, of the order of one. To remind the reader, such an algorithm consists of  *flow augmentations*, which are repeated till the min-cut is formed (by a surface of blocked bonds) and maxflow is reached. For each augmentation one needs to establish a path from the “source” to the “sink”, which can be done e. g. . by using shortest-distance path methods ADME01 . The one-dimensional nature means that the backbone can be decomposed into strings of subsequent cutting bonds $`C_i`$ and and blobs separating such strings $`B_i`$. Thus the structure is equivalent to the one-dimensional series $`\mathrm{}C_iB_iC_{i+1}\mathrm{}`$. In principle one may thus write a more efficient algorithm by abandoning the lattice structure, and describing the internal geometry of each $`B_i`$ separately. Thus an optimal version of the algorithm would entail the following steps:  *i):* establish the structure ($`B_i`$, $`C_i`$).  *ii):* Find an augmenting path along the chain, across all $`B_i`$.  *iii):* Augment flow, that is: find the smallest capacity in the $`C_i`$, and the smallest capacity in all the $`B_i`$. This is $`f_1`$.  *iv):* if $`f_1`$ equals the minimal cutting bond capacity stop, otherwise augment (subtract $`f_1`$ from $`C_i`$, and the paths inside $`B_i`$).  *v):* update the paths inside those $`B_j`$, only, where a bond was saturated by $`f_i`$ ($`i=1`$ to begin with). Goto  *iii)*. We have used, instead, an Euclidean background for the maxflow part, since the scaling of the matching program is indeed the bottleneck.