text:
# Infinite-Order Percolation and Giant Fluctuations in a Protein Interaction Network

## Abstract

We investigate a model protein interaction network whose links represent interactions between individual proteins. This network evolves by the functional duplication of proteins, supplemented by random link addition to account for mutations. When link addition is dominant, an infinite-order percolation transition arises as a function of the addition rate. In the opposite limit of high duplication rate, the network exhibits giant structural fluctuations in different realizations. For biologically-relevant growth rates, the node degree distribution has an algebraic tail with a peculiar rate dependence for the associated exponent.

PACS numbers: 02.50.Cw, 05.40.-a, 05.50.+q, 87.18.Sn

Inter-protein interactions underlie the performance of vital biological functions. Organisms with sequenced genomes, such as the yeast S. cerevisiae, provide important test beds for analyzing protein interaction networks. The number of interactions per protein of S. cerevisiae follows a power-law , a feature common to many complex networks, such as the Internet, the world-wide web, and metabolic networks. Similar behavior is exhibited by protein interaction networks of various bacteria. Based on the observational data, simple proteome growth models have recently been formulated to account for the evolution of this interaction network, where proteins are viewed as the nodes of a graph and links connect functionally related proteins.

In this work, we determine the structure of a minimal protein interaction network model that evolves by the biologically-inspired processes of protein duplication and subsequent mutation. That is, the functionality of a duplicate protein is similar, but not identical, to the original and can gradually evolve with time due to mutations. Within a rate equation approach, we show that: (i) the system undergoes an infinite-order percolation transition as a function of mutation rate, with a rate-dependent power-law cluster-size distribution everywhere below the threshold, (ii) there are giant fluctuations in network structure and no self-averaging for large duplication rate, and (iii) the degree distribution has an algebraic tail with a peculiar rate-dependent exponent when the duplication and mutation rates have biologically realistic values. Some aspects of this last result were recently seen .

In the model, nodes are added sequentially and the new node duplicates a randomly chosen pre-existing “target” node, viz., the new node links to each of the neighbors of the target with probability $`1\delta `$; each new node also links to any previous node with probability $`\beta /N`$, where $`N`$ is the current total number of nodes (Fig. 1). Thus an arbitrary number of clusters can merge when a single node is introduced. As we now discuss, this unusual dynamics appears to be responsible for the unconventional percolation properties of this network in the limit of zero duplication rate but finite mutation rate ($`\delta =0`$, $`\beta >0`$).

Let $`C_s(N)`$ be the expected number of clusters of size $`s1`$. This cluster size distribution obeys the rate equation

$$\frac{dC_s}{dN}=\beta \frac{sC_s}{N}+\underset{n=0}{\overset{\mathrm{}}{}}\frac{\beta ^n}{n!}e^\beta \underset{s_1\mathrm{}s_n}{}\underset{j=1}{\overset{n}{}}\frac{s_jC_{s_j}}{N},$$
(1)
where the sum is over all $`s_11,\mathrm{},s_n1`$ such that $`s_1+\mathrm{}+s_n+1=s`$. The first term on the right-hand side of Eq. (1) accounts for the loss of $`C_s`$ due to the linking of a cluster of size $`s`$ with the newly-introduced node. The gain term accounts for all merging processes of $`n`$ initially separated clusters whose total size is $`s1`$.

Solving for the first few $`C_s(N)`$, we see that they are all proportional to $`N`$. Thus writing $`C_s(N)=Nc_s`$, and introducing the generating function $`g(z)=_{s1}sc_se^{sz}`$, Eq. (1) becomes

$$g=\beta g^{}+(1+\beta g^{})e^{z+\beta (g1)},$$
(2)
where $`g^{}=dg/dz`$. To detect the percolation transition, we use the fact that $`g(0)=sc_s`$ is the fraction of nodes within finite clusters. Thus the size of the infinite cluster (the giant component) is $`NG=N(1g(0))`$. Suppose that we are in the non-percolating phase; this means that $`g(0)=1`$. In this regime, the average cluster size equals $`s=s^2c_s=g^{}(0)`$. To determine $`g^{}(0)`$, we substitute the expansion $`g(z)=1+zg^{}(0)+\mathrm{}`$ into Eq. (2) and take the $`z0`$ limit. This yields a quadratic equation for $`g^{}(0)`$ with solution

$$g^{}(0)=s=\frac{12\beta \sqrt{14\beta }}{2\beta ^2}.$$
(3)
This has a real solution only for $`\beta 1/4`$, thus identifying the percolation threshold as $`\beta _c=1/4`$. For $`\beta >\beta _c`$, we express $`g^{}(0)`$ in terms of the size of the giant component by setting $`z=0`$ in Eq. (2) to give

$$g^{}(0)=\frac{e^{\beta G}+G1}{\beta \left(1e^{\beta G}\right)}.$$
(4)
When $`\beta \beta _c`$, we use $`G0`$ to simplify Eq. (4) and find $`s(1\beta _c)\beta _c^2=12`$. On the other hand, Eq. (3) shows that $`s4`$ when $`\beta \beta _c`$ from below. Thus the average size of the finite clusters jumps discontinuously from 4 to 12 as $`\beta `$ passes through $`\beta _c=\frac{1}{4}`$.

The cluster size distribution $`c_s`$ exhibits distinct behaviors below, at, and above the percolation transition. For $`\beta <\beta _c`$, the asymptotic behavior of $`c_s`$ can be read off from the behavior of the generating function as $`z0`$. If $`c_s`$ has the power-law behavior

$$c_sBs^\tau \mathrm{as}s\mathrm{},$$
(5)
then the corresponding generating function $`g(z)`$ has the following small-$`z`$ expansion

$`g(z)=1+g^{}(0)z+B\mathrm{\Gamma }(2\tau )(z)^{\tau 2}+\mathrm{}.`$ (6)
The regular terms are needed to reproduce the known zeroth and first derivatives of the generating function, while the asymptotic behavior is controlled by the dominant singular term $`(z)^{\tau 2}`$. Higher-order regular terms are asymptotically irrelevant. Substituting this expansion into Eq. (2) we find that the dominant terms are of the order of $`(z)^{\tau 3}`$. Balancing all contributions of this order gives

$$\tau =1+\frac{2}{1\sqrt{14\beta }}.$$
(7)
Intriguingly, a power-law cluster size distribution with a non-universal exponent arises for all $`\beta <\beta _c`$. In contrast to ordinary critical phenomena, the entire range $`\beta <\beta _c`$ is critical.

The power-law tail implies that the size of the largest cluster $`s_{\mathrm{max}}`$ grows as a power law of the system size. From the extreme statistics criterion $`_{ss_{\mathrm{max}}}Nc_s=1`$ and the asymptotics of Eq. (5), we find $`s_{\mathrm{max}}N^{1/(\tau 1)}`$, or $`s_{\mathrm{max}}N^{\frac{1}{2}\sqrt{\beta _c\beta }}`$. In contrast, for conventional percolation below threshold, the largest cluster has size $`s_{\mathrm{max}}\mathrm{ln}N`$, reflecting the exponential tail of the cluster size distribution .

At the transition, Eq. (7) gives $`\tau =3`$. However, the naive asymptotics $`c_ss^3`$ cannot be correct as it implies that $`g^{}(0)`$ diverges. Similarly, we cannot expand the generating function as in Eq. (6) with $`\tau =3`$, since the singular term $`\mathrm{\Gamma }(1)\times (z)`$ has an infinite prefactor. As in other situations where the order of a singular term coincides with a regular term, we anticipate a logarithmic correction. Thus consider the modified expansion $`g(z)=1+4z+zu(z)+\mathrm{}`$, where $`u(z)`$ vanishes slower than any power of $`z`$, as $`z0`$. Substituting this into Eq. (2), setting $`\beta =\beta _c`$, and equating singular terms yields $`(8+u)zu^{}+u^2=0`$. Solving this differential equation asymptotically we obtain the leading behavior $`u8/\mathrm{ln}(z)`$; this indeed vanishes slower than any power of $`z`$ for $`z0`$. Substituting this form for $`u(z)`$ in the modified expansion for $`g(z)`$ and inverting yields

$$c_s\frac{8}{s^3(\mathrm{ln}s)^2}\mathrm{as}s\mathrm{}.$$
(8)
Thus exactly at the transition, the cluster size distribution acquires a logarithmic correction. This result also implies that the size of the largest component scales as $`s_{\mathrm{max}}N^{1/2}/\mathrm{ln}N`$.

Above the percolation transition, both $`g(0)=1G`$ and $`g^{}(0)`$ (Eq. (4)) are finite, so that the expansion for $`g(z)`$ has the form $`g(z)=1G+g^{}(0)z+\mathrm{}`$. Substituting this into Eq. (2) one can show that: (i) the full expansion of $`g(z)`$ is regular in $`z`$, and (ii) the generating function diverges at $`z_{}=1/s_{}`$. This latter fact implies that $`c_se^{s/s_{}}`$ as $`s\mathrm{}`$. The location of the singularity is determined by the condition $`e^{z+\beta (g1)}=1`$. This gives $`s_{}16/G`$ as $`\beta \beta _c`$. Realistic protein interaction networks are always above the percolation transition, e.g., for yeast the giant component includes 54% of all nodes and 68% of the links of the system; thus a giant component always exists and the cluster-size distribution has an exponential tail.

The size of the giant component $`G(\beta )`$ is obtained by solving Eq. (2) near $`z=0`$. A lengthy analysis shows that near the percolation threshold:

$$G(\beta )\mathrm{exp}\left(\frac{\pi }{\sqrt{4\beta 1}}\right),$$
(9)
so that all derivatives of $`G(\beta )`$ vanish as $`\beta \beta _c`$. Thus the transition is of infinite order. Similar behavior has been recently observed  for several growing network models where single nodes and links were introduced independently. This generic growth mechanism seems to give rise to fundamentally new percolation phenomena.

We now examine the complementary limit of no mutations ($`\beta =0`$) and show that individual realizations of the evolution lead to widely differing results. Consider first the limit of deterministic duplication of $`\delta =0`$ where all the links of the duplicated protein are completed. There is still a stochastic element in this growth, as the node to be duplicated is chosen randomly. When $`\delta =0`$, the rate equation approach \[Eqs. (14)–(15) below\] predicts that the degree distribution $`N_k`$ (defined as the number of nodes that are linked to $`k`$ other nodes) is given by $`N_k=2(12/N)^{k1}`$.

However, this “solution” does not correspond to the outcome of any single realization of the duplication process. To appreciate this, consider the simple and generic initial state of two nodes that are joined by a single link. We denote this graph as $`K_{1,1}`$, following the graph theoretic terminology that $`K_{n,m}`$ denotes a complete bipartite graph in which every node in the subgraph of size $`n`$ is linked to every node in the subgraph of size $`m`$. Duplicating one of the nodes in $`K_{1,1}`$ gives $`K_{2,1}`$ or $`K_{1,2}`$, equiprobably. By continuing to duplicate nodes, one finds that at every stage the network always remains a complete bipartite graph, say $`K_{k,Nk}`$, and that every value of $`k=1,\mathrm{},N1`$ occurs with equal probability (Fig. 2). Thus the degree distribution remains singular – it is always the sum of two delta functions!

For fixed $`N`$, we average over all realizations of the evolution to obtain the average degree distribution

$$N_k=2\left(1\frac{k1}{N1}\right).$$
(10)
Computing $`N_k`$ for other generic initial conditions, e.g., complete $`m`$-partite graphs and ring graphs, we find that the initial condition dependence persists throughout the evolution. More importantly, self-averaging breaks down: different realizations of the growth lead to statistically distinguishable networks. Similar giant fluctuations arise in the general case of imperfect duplication where $`\beta =0`$ and $`\delta >0`$ . To illustrate the origin of these macroscopic fluctuations, consider the network growth in the limit $`\delta 1`$. The probability that the first few duplication steps are complete (all eligible links are created) is close to one. For this initial development, the degrees of each node increase and the probability to create isolated nodes becomes very small as the network grows. On the other hand, if the first duplication event was totally incomplete, an isolated node would be created. The creation of isolated nodes necessarily leads to more isolated nodes but subsequent duplication events. Thus the number of isolated nodes is a non-self-averaging quantity. In a similar fashion, the number of nodes of degree $`k`$ for any finite $`k>0`$ is also non-self-averaging.

Finally, we investigate to the evolution of the network when both incomplete duplication and mutation occur ($`\delta <1`$ and $`\beta >0`$). Let us first determine the average node degree of the network, $`𝒟`$, for such general rates. In each growth step, the average number of links $`L`$ increases by $`\beta +(1\delta )𝒟`$. Therefore, $`L=[\beta +(1\delta )𝒟]N`$. Combining this with $`𝒟=2L/N`$ gives

$$𝒟=\frac{2\beta }{2\delta 1},$$
(11)
a result that applies only when $`\delta >\delta _c=1/2`$. Below this threshold, the number of links grows as

$$\frac{dL}{dN}=\beta +2(1\delta )\frac{L}{N},$$
(12)
and combining with $`𝒟(N)=2L(N)/N`$, we find

$$𝒟(N)=\{\begin{array}{cc}\mathrm{finite}\hfill & \delta >1/2\text{,}\hfill \\ \beta \mathrm{ln}N\hfill & \delta =1/2\text{,}\hfill \\ \mathrm{const}.\times N^{12\delta }\hfill & \delta <1/2\text{.}\hfill \end{array}$$
(13)
Without mutation ($`\beta =0`$) the average node degree always scales as $`N^{12\delta }`$, so that a realistic finite average degree is recovered only when $`\delta =1/2`$. Thus mutations play a constructive role, as a finite average degree arises for any duplication rate $`\delta >1/2`$.

We now consider this case of $`\delta >1/2`$ and $`\beta >0`$ and apply the rate equation approach  to study the degree distribution $`N_k(N)`$. The degree $`k`$ of a node increases by one at a rate $`A_k=(1\delta )k+\beta `$. The first term arises because of the contribution from duplication, while mutation leads to the $`k`$-independent contribution. The rate equations for the degree distribution are therefore

$$\frac{dN_k}{dN}=\frac{A_{k1}N_{k1}A_kN_k}{N}+G_k.$$
(14)
The first two terms account for processes in which the node degree increases by one. The source term $`G_k`$ describes the introduction of a new node of $`k`$ links, with $`a`$ of these links created by duplication and $`b=ka`$ created by mutation. The probability of the former is $`_{sa}n_s\left(\genfrac{}{}{0pt}{}{s}{a}\right)(1\delta )^a\delta ^{sa}`$, where $`n_s=N_s/N`$ is the probability that a node of degree $`s`$ is chosen for duplication, while the probability of the latter is $`\beta ^be^\beta /b!`$. Since duplication and random attachment are independent processes, the source term is

$$G_k=\underset{a+b=k}{}\underset{s=a}{\overset{\mathrm{}}{}}n_s\left(\genfrac{}{}{0pt}{}{s}{a}\right)(1\delta )^a\delta ^{sa}\frac{\beta ^b}{b!}e^\beta .$$
(15)
From Eq. (14), the $`N_k`$ grow linearly with $`N`$. Substituting $`N_k(N)=Nn_k`$ in the rate equations yields

$$\left(k+\frac{\beta +1}{1\delta }\right)n_k=\left(k1+\frac{\beta }{1\delta }\right)n_{k1}+\frac{G_k}{1\delta }.$$
(16)
Since $`G_k`$ depends on $`n_s`$ for all $`sk`$, the above equation is not a recursion. However, for large $`k`$, we can reduce it to a recursion by simple approximations. As $`k\mathrm{}`$, the main contribution to the sum in Eq. (15) arises when $`b`$ is small, so that $`a`$ is close to $`k`$, and the summand is sharply peaked around $`sk/(1\delta )`$. This simplifies the sum, as we may replace the lower limit by $`s=k`$, and $`n_s`$ by its value at $`s=k/(1\delta )`$. Further, if $`n_k`$ decays as $`k^\gamma `$, we write $`n_s=(1\delta )^\gamma n_k`$ and simplify $`G_k`$ to

$`G_k`$ $``$ $`(1\delta )^\gamma n_k{\displaystyle \underset{s=k}{\overset{\mathrm{}}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{s}{k}}\right)(1\delta )^k\delta ^{sk}{\displaystyle \underset{b=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\beta ^b}{b!}}e^\beta `$ (17)
$`=`$ $`(1\delta )^{\gamma 1}n_k,`$ (18)
since the former binomial sum equals $`(1\delta )^1`$.

Thus for $`k\mathrm{}`$, Eq. (16) reduces to a recursion relation, from which we deduce that $`n_k`$ has the power-law behavior $`k^\gamma `$, with $`\gamma `$ determined from the relation

$$\gamma (\delta )=1+\frac{1}{1\delta }(1\delta )^{\gamma 2}.$$
(19)
Notice that the replacement of $`n_s`$ by $`(1\delta )^\gamma n_k`$ is valid only asymptotically. This explains the slow convergence of the degree distribution to the predicted power law form (Fig. 3). Intriguingly, the exponent $`\gamma (\delta )`$ is independent of the mutation rate $`\beta `$ . Nevertheless, the presence of mutations ($`\beta >0`$) is vital to suppress the non-self-averaging as the network evolves and thus make possible a smooth degree distribution. If we adopt $`\delta =0.53`$, as suggested by observations, we obtain $`\gamma =2.373\mathrm{}`$, compared to the numerical simulation result of $`\gamma =2.5\pm 0.1`$.

In summary, network growth by duplication and mutation leads to rich behavior with an infinite-order percolation transition and no self-averaging in the absence of mutations. Without mutation, different realizations of the network lead to drastically different outcomes and each outcome is itself singular. Mutations are needed to form networks that are statistically similar to observed protein interaction networks. Thus mutations seem to play a constructive role in forming robust networks whose functioning realizes the primary purpose of mutations.

We thank NSF grant DMR9978902, KOSEF grant 2002-2-11200-002-3 in the BRP program, and a travel grant from the BK21 project.


format: text

text:
# Forbidden activation levels in a non-stationary tunneling process

\[

## Abstract

> Tunneling in the presence of an opaque barrier, part of which varies in time, is investigated numerically and analytically in one dimension. Clearly, due to the varying barrier a tunneling particle experiences spectral widening. However, in the case of strong perturbations, the particles’ activation to certain energies is avoided. We show that this effect occurs only when the perturbation decays faster than $`t^2`$

\]

Tunneling, which is a very common phenomenon in the quantum world, is rarely a stationary process. In most cases, it is involved with temporal changes of the medium through which the particle is tunneling. In some cases, such variations are mostly due to transport of the tunneling particles themselves. In these cases, the problems are intrinsically nonlinear. However, when the source of the changes is external, the problem can be treated as a linear one. In this category, one can even find stationary many-body transport problems where a single particle can experience the influence of its surroundings (in this case its neighbors) as an alternating external force. Many works (see, for example, refs. ) done in this field have revealed an abundance of physical effects: resonances and oscillations in energy, exponential increase of the tunneling current, activation-assisted tunneling and elevator resonance activation.

This paper discusses the influence of a very strong perturbation (compared with the incoming particles’ energy) on the tunneling process. It is assumed that the perturbation potential is always positive (see Fig.1), i.e., it increases the barrier (locally) and never creates a local well (inside the barrier). Therefore, effects like elevator resonance activation do not occur (see ref. ). On the contrary, in the adiabatic approximation, the perturbation - being very strong - should reduce the transmission dramatically. However, when the characteristic time of the perturbation decreases, the particle may absorb energy from the external perturbation, which can assist it in its due course of tunneling. It is clear, for example, that in the extreme case, where the potential barrier is smaller than $`\mathrm{}/\tau `$ ($`\tau `$ is the characteristic time of the perturbation), the particle can easily (with high probability) absorb enough energy to leap over the barrier. This paper discusses the intermediate case, where the process has a tunneling nature but allows for energy quanta to be absorbed by the tunneling particle. A very strong perturbation, which is activated for a finite period, permits propagation only at the very onset and at the final decay of the perturbation. When the wave function in these two instances has the same phase, maximum transmission is expected, but when in these two moments the wave function is out of phase, one should expect minimum transmission. These are the forbidden activation energies we are looking for.

This system, which is presented in Fig.1, can be described in terms of the Schrödinger equation in the following way

$$\frac{d^2\psi }{dx^2}+U(x)\psi +f(t/\tau )\delta (x)\psi =i\frac{d\psi }{dt}$$
(1)
where we have used the units $`\mathrm{}=1`$ and $`2m=1`$ ($`m`$ is the electron’s mass), and $`f(x)`$ is a function, which vanishes very rapidly for $`x\mathrm{}`$, and $`U(x)`$ is the barrier’s potential.

We seek the solution

$$\psi (x,t)=\psi _\mathrm{\Omega }(x,t)+𝑑\omega e^{i\omega t}a_\omega G_\omega ^+(x)$$
(2)
where $`\psi _\mathrm{\Omega }(x,t)`$ is the homogeneous solution of eq.1, $`\psi _\mathrm{\Omega }(x,t)\phi _\mathrm{\Omega }(x)e^{i\mathrm{\Omega }t}`$, while $`\phi `$ solves the stationary-state equation

$$\frac{d^2\phi _\mathrm{\Omega }}{dx^2}+\left[U(x)\mathrm{\Omega }\right]\phi _\mathrm{\Omega }=0;$$
(3)
$`G_\omega ^+`$ is the ”outgoing” Green function of equation 3 (but with a general $`\omega `$ instead of $`\mathrm{\Omega }`$).

The perturbation can be regarded as the force that is activated upon the particle for a finite period $`\tau `$ ($`f`$ has momentum dimensions, therefore $`f/\tau `$ has dimensions of force).

A straightforward substitution of eq.2 in eq.1 reveals the integral equation

$$\stackrel{~}{f}(\omega \mathrm{\Omega })\phi _\mathrm{\Omega }(0)+a_\omega 𝑑\omega ^{}\stackrel{~}{f}(\omega \omega ^{})a_\omega ^{}G_\omega ^{}^+(0)=0$$
(4)
where

$$\stackrel{~}{f}(\omega )(2\pi )^1𝑑tf(t/\tau )\mathrm{exp}(i\omega t)$$
(5)
and some tedious calculations  reveal that, for a rectangular barrier with potential height $`U`$ and width $`L`$,

$$G_\omega (0)=\frac{\mathrm{coth}\left[\kappa L+i\mathrm{arctan}(k/\kappa )\right]}{2\kappa }(2\kappa )^1$$
(6)
where $`k\sqrt{\omega }`$ and $`\kappa \sqrt{U\omega }`$ and the approximation on the right hand side takes place for an opaque barrier $`\left(\sqrt{U}L1\right)`$.

It should be noted that the results of this paper are valid for any opaque barrier. We have chosen a rectangular barrier since it has a relatively simple expression. In fact, the approximation in eq.6 (which is a good estimation for any opaque barrier) yields, for all practical purposes, the same results.

In general, eq.4 does not have an analytical solution. Thus, we will begin with the numerical results. In Fig. 2 we present the results of the exact solution of eq.4 for the temporally confined perturbation

$$f(t/\tau )=(f_0/\tau )\mathrm{exp}(|t/\tau |).$$
(7)
The $`x`$-axis represents the frequency $`\omega `$ (i.e., the frequency of the transmitted particles), and the y-axis represents the logarithm of $`\tau `$ (the perturbation’s characteristic time). Within this matrix, the gray levels indicate the logarithm of the transmission coefficient, i.e.,$`\mathrm{ln}|a_\omega G_\omega ^+(x>L)|`$. In the figure, dark spots indicate a high probability of a particle being emitted with the corresponding frequency, or energy ($`\omega `$).

From the figure we see that, for very large $`\tau `$ (i.e., the adiabatic case), the outgoing frequency is almost identical to the incoming one, $`\mathrm{\Omega }`$. That is, the particles simply experience spectral widening.

When $`\tau `$ decreases (and thus the process becomes more energetic), the spectral width of the outgoing particles increases, as it should by the uncertainty principle. The strange part comes later, when $`\tau `$ reaches the value $`\tau _0`$ , after which the increase in the spectral width becomes more moderate. This is the point at which the activated particles prefer to be activated to the first resonant energy $`U+(\pi /2L)^2`$. These resonant energies *are* barrier-dependent; therefore, one should expect differences in these resonances from one barrier to the next. However, these resonances have nothing to do with the following effect, which is the main issue in this paper.

When $`\tau `$ decreases even further, the spectral width is split into discrete probable frequencies-like a fan. That is, the transmitting (through the barrier) particles are activated to specific preferred frequencies, while *activation to others frequencies is forbidden*.

To explain this behavior, we can use the approximation that $`U\mathrm{\Omega }\tau ^1`$. Since, $`U\mathrm{\Omega }`$ can be arbitrarily large, this approximation still enables us to discuss relatively rapid changes. Within this temporal regime, only the spectral vicinity of $`\mathrm{\Omega }`$ affects the amplitude $`a_\omega `$. Thus, we can approximate the Green function in eq. 4 according to: $`G_\omega ^+(0)G_\mathrm{\Omega }^+(0)`$. The important thing, which allows for this exchange, is that when $`\mathrm{\Omega }`$ is very far from $`U`$, $`G_\omega ^+(0)(2\sqrt{U\omega })^1`$ does not change much within the spectral range (again, recall that this approximation is almost barrier-independent). Hence, with this assumption in hand eq.4 can readily be solved with a Fourier transform:

$$A(t)=\frac{f(t/\tau )\phi _\mathrm{\Omega }(0)\mathrm{exp}(i\mathrm{\Omega }t)}{1f(t/\tau )G_\mathrm{\Omega }^+(0)}$$
(8)
where $`A(t)`$ is the inverse Fourier counterpart of $`a_\omega `$. For convenience, we use the dimensionless function $`g`$: $`f(t/\tau )(f_0/\tau )g(t/\tau )`$ (where $`f_0f(0)`$). The results in Fig.2, for example, were taken for $`g(x)=\mathrm{exp}(|x|)`$. It is assumed that $`f_0f(0)`$ is the maximal value of the perturbation. Therefore,

$$A(t)=\frac{f_0g(t/\tau )\phi _\mathrm{\Omega }(0)\mathrm{exp}(i\mathrm{\Omega }t)}{\tau f_0g(t/\tau )G_\mathrm{\Omega }^+(0)}$$
(9)
Eq.9 determines two time scales

$$\tau _1f_0G_\mathrm{\Omega }$$
(10)
and

$$\tau _2\tau g^1(f_0G_\mathrm{\Omega }/\tau ).$$
(11)
When $`\tau \tau _1`$, the denominator of eq. 9 can be ignored, and the broadening of the spectral width is exactly equal to the spectral width of the perturbation. Note that since we are dealing with a strong perturbation (large $`f_0`$), this characteristic time is very long.

When $`\tau <\tau _1`$ (and this is the more interesting regime), the denominator cannot be ignored, and here $`\tau _2`$ comes into play: when $`t\tau _2`$, $`A(t)\phi _\mathrm{\Omega }(0)\mathrm{exp}(i\mathrm{\Omega }t)/G_\mathrm{\Omega }^+(0)`$; except for the oscillating term the absolute value of $`A(t)`$ is almost a constant. However, when $`t\tau _2`$, the denominator of eq. 9 can again be ignored (because of the rapid decay of $`g`$). Therefore, $`A(t)`$ has a kind of rectangular shape: it is almost a constant for a time $`\tau _2\tau `$ and then it abruptly decays with the time scale $`\tau `$ For such a ”rectangular” shape, the Fourier transform is something like the sinc function, which is oscillating and decaying at the same time. This oscillating ”frequency” of $`a_\omega `$ is $`\tau _2`$ (note that we are dealing here with the frequency space, hence the frequency has dimensions of time). Thus, one can write $`a_\omega a[(\omega \mathrm{\Omega })\tau _2]`$. Multiplication by the Green function $`G_\omega ^+(x>L)`$ distorts the picture a little (it increases the high frequencies’ contributions at the expense of those of the low frequencies), but it does not change the fact that the particles can be activated (”act”) only by a product of an integer and the quantity $`\delta \omega \tau _2^1`$. That is,

$$\omega _{act}^n\mathrm{\Omega }+n\delta \omega ,$$
(12)
where $`n`$ is an integer. Similarly, no activation occurs to the forbidden (”for”) activation levels

$$\omega _{for}^n\mathrm{\Omega }+\left(n+\frac{1}{2}\right)\delta \omega .$$
(13)
These are the forbidden frequencies we were looking for. From the definition of $`\tau _2`$ (eq. 11), one learns that it depends not only on the characteristic time $`\tau `$ and on the incoming frequency $`\mathrm{\Omega }`$, but also on the strength of the perturbation $`f_0`$ and, maybe more surprisingly, on the *functional behavior* of $`f`$. For example, suppose

$$g(x)=\mathrm{exp}(|x|^m)$$
(14)
then

$$\tau _2=\tau \left[\mathrm{ln}\left(\frac{f_0G_\mathrm{\Omega }}{\tau }\right)\right]^{1/m}$$
(15)
It is clear that when the perturbation is strong ($`f_0`$ is large) the dependence on $`m`$ can be very significant.

Let us take the following temporal dependence that corresponds to a milder decay:

$$g_n(x)=\frac{1}{1+|x|^n}$$
(16)
Such temporal behavior leads to a particularly simple solution, which allows tracing the onset of the ”discrete activation level effect”. In this case

$$a_\omega =\frac{1}{2\pi }𝑑tA(t)\mathrm{exp}(i\omega t)=\frac{f_0}{2\pi }𝑑x\frac{\mathrm{exp}(ixs)}{|x|^n+\alpha ^n}$$
(17)
where we adopt the dimensionless parameters $`s(\omega \mathrm{\Omega })\tau `$, $`\alpha ^nG(0)f_0/\tau +1`$, and of course $`xt/\tau `$.

But eq.17 can also be rewritten

$$a_\omega =\frac{f_0}{\alpha ^{n1}}\stackrel{~}{g}_n(s\alpha )$$
(18)
where $`\stackrel{~}{g}_n`$ is the inverse Fourier transform of $`g_n`$ (see eq.5).

For $`y\mathrm{}`$, an asymptotic expression can be written

$$\stackrel{~}{g}_n(y)\mathrm{cos}[y\mathrm{cos}(\pi /n)]\mathrm{exp}[y\mathrm{sin}(\pi /n)],$$
(19)
(for $`n=2`$ and $`n=4`$ this is the only term in the asymptotic expansion).

Thus, the oscillating part of $`a_\omega `$ reads

$$\mathrm{cos}\left[(\omega \mathrm{\Omega })\tau \left(\frac{Gf_0}{\tau }+1\right)^{1/n}\mathrm{cos}\left(\frac{\pi }{n}\right)\right]$$
(20)
Similar to the exponential case (eq. 14) we can define a characteristic time

$$\tau _2\tau (Gf_0/\tau +1)^{1/n}\mathrm{cos}(\pi /n)$$
(21)
which determines the activation energies according to eq. 12. Similar to eq. 15, $`\tau _2\tau `$ when the *functional* behavior of the perturbation dies out very quickly ($`n\mathrm{}`$). Moreover, eq. 21 also reveals the onset of the discrete activation pattern. The parameter $`n`$ *must be larger than 2*, otherwise no oscillations occur.

Therefore, we conclude that in order to obtain the discrete activation level effect, the perturbation must die out faster than $`(t/\tau )^2`$.

It should be noted, that because of the point-like nature of the alternating impurity, the whole derivation of the effect depends on the barrier via the Green function, however, as approximation (6) suggests, any opaque barrier would yield similar results.

To summarise, the activation energy in the case of a strongly varying potential barrier was investigated. It was shown that when the temporal perturbation amplitude is very large, the particles cannot be activated to certain discrete energies $`\omega _{for}^m\mathrm{\Omega }+(m+1/2)\delta \omega `$, where $`\delta \omega [\tau g^1(f_0G_\mathrm{\Omega }/\tau )]^1`$, $`\tau `$ is the characteristic time of the perturbation, $`f_0`$ is its strength and $`g^1`$ is the inverse function of the perturbation’s temporal dependence. Moreover, it was demonstrated that in order to see this forbidden activation effect, the perturbation must die out faster than $`t^2`$ .

I would like to thank Prof. Mark Azbel for stimulating discussions.


text:
# Correlations in proton-proton collisions with ALICE

## 1 Particle correlations in pseudo-rapidity

Particle correlations and particle multiplicity distributions (MD’s) cannot be approached independently: a unified description of correlations and multiplicity distributions is always needed in order to understand the underlying dynamics in high energy collisions. The warning comes from charged particle $`H_q`$ moments oscillations and shoulder structure in charged particle multiplicity distributions both in $`\mathrm{p}\overline{\mathrm{p}}`$ collisions and e+e- annihilation, which have been shown to have a common origin in the weighted superposition of different classes of events or topologies, each one described by the same multiplicity distribution, i.e, a negative binomial (NB) (Pascal) MD with characteristic parameters for each substructure . Accordingly, this Section should be considered the natural continuation of Ref.  on global event properties.

Attention here will be on differential variables whose integrated versions have been already discussed in the just mentioned Reference. Indeed, $`n`$-particle exclusive cross-section density

$$\frac{1}{\sigma _{\text{inel}}}\frac{d^n\sigma _n}{d\eta _1\mathrm{}d\eta _n}$$
(1)
integrated over $`\eta _1,\mathrm{},\eta _n`$, corresponds to $`n`$-particle multiplicity distribution $`P_n=\sigma _n/\sigma _{\text{inel}}`$ , with $`\sigma _n`$ the $`n`$-particle topological cross-section and $`_n\sigma _n=\sigma _{\text{inel}}`$; and $`n`$-particle inclusive cross-section density

$$\frac{1}{\sigma _{\text{inel}}}\frac{d^n\sigma }{d\eta _1\mathrm{}d\eta _n}=\rho _n(\eta _1\mathrm{}\eta _n)$$
(2)
integrated over $`\eta _1,\mathrm{},\eta _n`$, corresponds to $`n`$-th order factorial moments, $`F_n`$, of the multiplicity distribution $`P_n`$. In order to avoid inessential contributions due to inclusive lower order combinatorics, the $`n`$-particle correlation functions $`C_n(\eta _1,\mathrm{},\eta _n)`$ are introduced. They are linked, via standard cluster expansion of statistical mechanics, to $`\rho _n(\eta _1\mathrm{}\eta _n)`$. Integrals of $`C_n(\eta _1,\mathrm{},\eta _n)`$ over $`\eta _1\mathrm{}\eta _n`$ define factorial cumulant moments, $`K_n`$, of the multiplicity distribution $`P_n`$. All the above mentioned variables are not normalised, their corresponding normalised version can be obtained by dividing by the product of the single inclusive cross-sections.

It should be pointed out that $`C_n(\eta _1\mathrm{}\eta _n)=0`$ indicates lack of $`n`$-particle correlations, a positive $`C_n(\eta _1\mathrm{}\eta _n)`$ the tendency of the $`n`$ particles to be correlated and to group together and a negative $`C_n(\eta _1\mathrm{}\eta _n)`$ anti-correlations among particles. They prefer to stay far apart one from the other.

Moreover the description of different classes of events in hadron-hadron collisions by NB (Pascal) MD’s would suggest two-particle correlation dominance in rapidity intervals as required by the hierarchical nature of $`n`$-particle correlations in each class of events .

## 2 Two-particle pseudo-rapidity correlations

Two-particle pseudo-rapidity correlations

$$C_2(\eta _1,\eta _2)=\rho _2(\eta _1,\eta _2)\rho _1(\eta _1)\rho _1(\eta _2)$$
(3)
have been shown to be strong at small rapidity distances and to depend on the particles’ charge combination. It should be reminded that by integrating $`C_2(\eta _1,\eta _2)`$ over $`\eta _1`$ and $`\eta _2`$ one gets $`D^2\overline{n}`$, where $`D`$ is the dispersion. In hadron-hadron collisions, the energy dependence of $`C_2(\eta _1=0,\eta )`$ from ISR  to UA5  energies is described in Fig. 1a.

Inclusive two-particle correlation functions $`C_2(\eta _1,\eta _2)`$ can be written in terms of semi-inclusive single- and two-particle densities for particles $`a`$ and $`b`$:

$$\rho _1^{(n)}(\eta )=\frac{1}{\sigma _n}\frac{d^n\sigma _n^a}{d\eta }$$
(4)
and

$$\rho _2^{(n)}(\eta _1,\eta _2)=\frac{1}{\sigma _n}\frac{d^n\sigma _n^{ab}}{d\eta _1d\eta _2}$$
(5)
by mixing events of different multiplicities $`n`$: two new functions can be defined:

$$C_2(\eta _1,\eta _2)=C_S(\eta _1,\eta _2)+C_L(\eta _1,\eta _2)$$
(6)
where

$$C_S(\eta _1,\eta _2)=\underset{n}{}P_nC_2^{(n)}(\eta _1,\eta _2)$$
(7)
and

$$C_L(\eta _1,\eta _2)=\underset{n}{}P_n[\rho _1^{(n)}(\eta _1)\rho _1(\eta _1)]\underset{n}{}P_n[\rho _1^{(n)}(\eta _2)\rho _1(\eta _2)]$$
(8)
with

$$C_2^{(n)}(\eta _1,\eta _2)=\rho _2^{(n)}(\eta _1,\eta _2)\rho _1^{(n)}(\eta _1)\rho _1^{(n)}(\eta _2)$$
(9)
i.e., the two-particle correlation function at fixed multiplicity $`n`$.

$`C_S`$ is the average of the semi-inclusive correlation functions and is more sensitive to dynamical correlations, and $`C_L`$ describes the mixing of different topological single particle densities. As shown in Fig. 1, $`C_S`$ does not depend significantly on energy. This fact should be contrasted with $`C_L`$ behaviour which is strongly energy dependent and is the result of mixing different topological single particle densities. In Fig. 1b and 1c are shown the energy dependence of $`C_S`$ and $`C_L`$: the contribution of $`C_L`$ to $`C_2(\eta _1,\eta _2)`$ is dominant with respect to $`C_S`$, in addition it increases significantly with energy whereas $`C_S`$ does not grow as much with energy around $`\eta =0`$: in addition $`C_S`$ is positive in the just mentioned region and becomes negative for $`|\eta |>1`$. The overall behaviour of $`C_2`$, $`C_L`$ and $`C_S`$ at fixed c.m. energy is summarised in Fig. 2.

It should be pointed out that

$$C_S(\eta _1,\eta _2)𝑑\eta _1𝑑\eta _2=\overline{n}$$
(10)
and

$$C_L(\eta _1,\eta _2)𝑑\eta _1𝑑\eta _2=D^2.$$
(11)
These equations clarify the meaning of the two contributions to total two particle correlations $`C_2(\eta _1,\eta _2)`$ and point out what quantity should be measured with ALICE in the allowed pseudo-rapidity range.

## 3 Bose-Einstein correlations

Among two-particle correlations, Bose-Einstein (BE) correlations are of particular interest in high energy collisions. The production of two identical bosons $`a`$, $`b`$ from two particle sources is controlled by an amplitude which is symmetrised with respect to the interchange of bosons $`a`$, $`b`$ and leads to an enhanced emission probability if bosons have similar momenta.

BE correlations are measured in terms of the second order normalised factorial cumulant $`R_2`$, i.e.,

$$R_2(q_a,q_b)=\frac{\rho _2(q_a,q_b)}{\rho _1(q_a)\rho _1(q_b)}1,$$
(12)
the ratio of two-particle inclusive cross-section $`\sigma ^1d^2\sigma /dq_adq_b`$ over the product of single- $`a`$ and $`b`$ particle inclusive cross-sections $`\sigma ^1d\sigma /dq_a`$, $`\sigma ^1d\sigma /dq_b`$.

$`R_2`$ function is directly related to the Fourier transform of the space-time distribution of particle production points. Accordingly, space distribution and lifetime of boson sources can be measured. In the standard Goldhaber parametrisation

$$R_2(Q)=\lambda \mathrm{exp}\left(G^2Q^2\right),$$
(13)
with $`Q^2=(q_aq_b)^2`$, the square of the four momentum difference of particles $`a`$ and $`b`$, \[$`q_{a,b}(\stackrel{}{p}_{a,b},E_{a,b})`$\]; $`\lambda `$ is the strength of the effect and $`G`$ is a measure of the source size.

The effect depends strongly on the masses of particles used in the analysis. It is well known that the range of the radius of the source from which pions are produced in e+e- annihilation at LEP lies between $`\mathrm{\hspace{0.17em}0.7}`$-$`1.0`$ fm, the range of the radius of the source from which $`\mathrm{K}`$ particles are emitted between $`\mathrm{\hspace{0.17em}0.5}`$-$`0.7`$ fm and that for $`\mathrm{\Lambda }`$ particles between $`\mathrm{\hspace{0.17em}0.1}`$-$`0.2`$ fm. In conclusion the most external shell produces particle pairs of lower mass (pions) and higher mass particles pairs (lambdas) are coming from most internal shells.

Are different particle production shells visible also in pp collisions in the $`\mathrm{TeV}`$ region? does a connection exist between low pair mass scale and short distance scale of QCD mini-jets? Extremely high ALICE statistics should allow to clarify all these points. Bose-Einstein effect depends also on $`p_\mathrm{t}`$ intervals considered and on related multiplicity densities . A sample of 2,400,000 non-single-diffractive $`\mathrm{p}\overline{\mathrm{p}}`$ events at c.m. energy 630 $`\mathrm{GeV}`$ measured by the UA1 central detector have been used. Only vertex associated charged tracks with $`p_\mathrm{t}>0.15`$ GeV$`/c`$, and $`|\eta |<3`$ have been considered. The azimuth angle has been restricted to $`45^\text{o}<|\varphi |<135^\text{o}`$.

Three subsamples for like-sign (ls) and opposite-sign (os) pion pairs have been selected with respect to pion transverse momentum $`p_\mathrm{t}`$:

* the all-$`p_\mathrm{t}`$ sample;
* a low-$`p_\mathrm{t}`$ subsample of charged particles with $`p_\mathrm{t}<0.7`$ GeV$`/c`$ (attention is paid in order to reduce the number of particles originated from jets or mini-jets);
* a high-$`p_\mathrm{t}`$ subsample with $`p_T>0.7`$ GeV$`/c`$ (particles come here predominantly from jets or mini-jets).

Then the behaviour of second order normalised cumulant correlation functions has been studied in the three subsamples for particle multiplicity densities $`dN/d\eta =1.22`$, $`2.72`$ and $`6.85`$. BE correlations show in hadron-hadron collisions a pronounced dependence on multiplicity in the BE strength parameter lambda or equivalently the cumulant of second order $`K_2`$ at $`Q=0`$. In particular the high-$`p_\mathrm{t}`$ sample cumulants exceeds one (which would indicate full coherence in pure BE correlations) and is hence concluded that this data sample is dominated by processes other than BE correlations (Fig. 3).

As far as the low-$`p_\mathrm{t}`$ sample is concerned it is shown that the dependence of the correlation strength and of higher order cumulants on multiplicity is important in order to test different theoretical models (Monte Carlos are totally inadequate here): for this task the low $`p_\mathrm{t}`$ cut-off of ALICE at LHC is required.

## 4 Higher moment correlations in small pseudo-rapidity intervals

The presence of very short range correlations and its connection with BE interference is indeed an important topic to be investigated, as in more general terms it is the search of local fluctuations of multiplicity distributions in momentum space and related scaling properties. Large concentrations of particle number in small pseudo-rapidity intervals for single events have been seen in the JACEE event for a pseudo-rapidity binning $`\delta \eta =0.1`$ (local fluctuations up to $`dn/d\eta `$ $`\mathrm{\hspace{0.17em}300}`$ with a signal-to-background ratio about 1:1)  and in the NA22 event (local fluctuation in rapidity of 60 times the average density) .

A possible explanation of these spikes was related to an intermittent behaviour, i.e., to the guess that there exist a correlations at all scales which implies a power low dependence of normalised factorial moments $`F_q(\delta )`$ on the size $`\delta `$ of the phase space bins:

$$\frac{F_q(\delta )}{n(\delta )^q}=\delta ^{\alpha _q};$$
(14)
the slope $`\alpha _q`$ was shown to increase as $`q`$ increases, see Fig. 4.

Since ALICE allows to study with its high statistics moments properties up to order $`q=14`$, it will be particularly interesting to analyse both intermittent behaviour as bin size decreases and the influence on this phenomenon on BE correlations together with resonance production.

If the above mentioned scaling law will be confirmed in the LHC energy domain a new horizon will be open on self-similar cascading structure and fractal properties of hadron hadron collisions, a fascinating perspective also for more complex collisions, like heavy ions.

## Acknowledgements

We thank J.-P. Revol for stimulating the present work.

## References


text:
# On the Lack of Type I X-ray Bursts in Black Hole X-ray Binaries: Evidence for the Event Horizon?

## 1 Introduction

When gas accretes onto a neutron star (NS) in a low-mass X-ray binary (LMXB), it is compressed and heated as it accumulates on the surface, leading to thermonuclear reactions. In many NS LMXBs, the reactions occur unsteadily and cause Type I X-ray bursts (Grindlay et al. 1976). Type I bursts have been observed in a large number of NS LMXBs (see Lewin, van Paradijs & Taam 1993; Strohmayer, Swank & Zhang 1998 for reviews), and the theory of these bursts is relatively well understood (Hansen & van Horn 1975; Woosley & Taam 1976; Joss 1977; Taam & Picklum 1978; Paczyński 1982, hereafter P82; Fujimoto et al. 1984, 1987; Fushiki & Lamb 1987, hereafter FL87; Taam, Woosley & Lamb 1996; Bildsten 1998).

Among the over 100 LMXBs known in the Galaxy, roughly a dozen systems, most of them transient sources (Tanaka & Shibazaki 1996), have been identified as black hole (BH) candidates. In these BH LMXBs, dynamical measurements give mass estimates for the accreting stars greater than the likely maximum mass $`3M_{}`$ of a NS (Shapiro & Teukolsky 1983, hereafter ST83; Narayan, Garcia & McClintock 2001; and references therein).

No Type I burst has been seen in a BH LMXB, even though, as we show in this Letter, BH LMXBs ought to produce bursts as efficiently as NS LMXBs if the accreting BH candidates possess surfaces. The lack of bursts thus represents possible evidence for the presence of event horizons. In §2 of this paper, we examine the stability of nuclear burning on the surface of a compact star. In §3, we discuss the predictions of the model for accretion onto a $`1.4M_{}`$ NS and a $`10M_{}`$ BH candidate with a hypothetical surface. We conclude in §4 with a discussion.

## 2 The Model

We consider a compact spherical star of mass $`M`$ and radius $`R`$, accreting gas steadily at a rate $`\dot{\mathrm{\Sigma }}(\mathrm{g}\mathrm{cm}^2\mathrm{s}^1)`$. In the local frame, the gravitational acceleration is $`g=GM(1+z)/R^2`$, where the redshift $`z`$ is given by $`1+z=(1R_S/R)^{1/2}`$, and $`R_S=2GM/c^2`$ is the Schwarzschild radius. We assume that the accreting material has mass fractions $`X_0`$, $`Y_0`$ and $`Z_0=1X_0Y_0`$, of hydrogen, helium and heavier elements (mostly CNO).

We consider a layer of accreted material of surface density $`\mathrm{\Sigma }_{max}`$ sitting on top of a substrate of fully burnt material ($`X=Y=0`$, $`Z=1`$). (The particular composition of the substrate is unimportant since we apply the inner boundary condition at its top.) Since the physical thickness of the accreted layer is much less than the radius, we work in plane parallel geometry and take $`g`$ to be independent of depth. We solve for the density $`\rho `$, the temperature $`T`$, the outgoing flux $`F`$, and the hydrogen, helium and heavy element fractions, $`X`$, $`Y`$, $`Z=1XY`$, as functions of the column density $`\mathrm{\Sigma }`$ ($`0\mathrm{\Sigma }\mathrm{\Sigma }_{max}`$).

The evolution equations for the gas in the layer are the standard equations of stellar physics. H- and He-burning give

$$\frac{dX}{dt}=\frac{ϵ_H}{E_H^{}},\frac{dY}{dt}=\frac{ϵ_H}{E_H^{}}\frac{ϵ_{He}}{E_{He}^{}},$$
(1)
where for our problem the total time derivative takes the form $`d/dt/t+\dot{\mathrm{\Sigma }}/\mathrm{\Sigma }`$. Here, $`ϵ_{H,He}`$ are the respective nuclear energy generation rates, and $`E_{H,He}^{}`$ are the corresponding energy release per unit mass of H and He burned (P82). For $`ϵ_H`$, we include the pp chain and the CNO cycle, including fast-CNO burning, saturated CNO burning, and electron capture reactions, as described in Mathews & Dietrich (1984) and Bildsten & Cumming (1998). Since we are not concerned with modelling the bursts themselves, and since our stability criterion does not depend on the detailed treatment of the deep crust, we do not include proton captures onto heavier nuclei; Schatz et al. (1999) illustrate some of the consequences of the $`rp`$-process burning on accreting neutron stars. For He-burning, we include the triple-$`\alpha `$ reaction, but not pycnonuclear reactions (e.g., ST83). We do not correct the reaction rates to include screening (e.g. FL87) since we are concerned only with determining whether nuclear burning of H and He can proceed stably under given conditions; stable burning of H and He utilizes almost exclusively the reactions included.

Hydrostatic equilibrium gives $`P/\mathrm{\Sigma }=g`$. For the pressure $`P`$ we use the expressions given in P82 for the gas, radiation and degeneracy pressure, along with a correction when the degenerate electrons become relativistic. Radiative transfer gives $`T/\mathrm{\Sigma }=(3\kappa F)/(16\sigma T^3)`$, where $`\sigma `$ is the Stefan-Boltzmann constant and we write the opacity $`\kappa `$ in the form $`1/\kappa =1/\kappa _{rad}+1/\kappa _{cond}`$. We employ Iben’s (1975) fitting functions for the radiative opacity $`\kappa _{rad}`$, and an analytical formula from Clayton (1968), suitably modified for relativistic electrons, for the conductive opacity $`\kappa _{cond}`$; the latter expression agrees well with more modern treatments (e.g. Heyl & Hernquist 2001). Finally, the energy equation gives

$$\rho T\frac{ds}{dt}=\rho (ϵ_H+ϵ_{He})+\rho \frac{F}{\mathrm{\Sigma }},$$
(2)
where $`s`$ is the entropy per unit mass. The above five equations form a closed set.

We have four outer boundary conditions at the surface of the star ($`\mathrm{\Sigma }=0`$). Two of these are $`(X,Y)=(X_0,Y_0)`$. We obtain the third boundary condition by equating the accretion luminosity of the infalling gas, $`L_{acc}=4\pi R^2\dot{\mathrm{\Sigma }}c^2z/(1+z)`$, to blackbody emission from the surface: $`4\pi R^2\sigma T_{out}^4=L_{acc}`$. This gives the surface temperature $`T_{out}`$. Then, using $`T_{out}`$ and an assumed value of $`F_{out}`$, we solve for the surface density profile $`\rho (\mathrm{\Sigma })`$ from the radiative transfer equation, thus obtaining the fourth boundary condition.

At the base of the accreted layer we have an inner boundary condition. We assume that the temperature at the top of the substrate, $`T_{in}=T(\mathrm{\Sigma }_{max})`$, is fixed. We examine several values of $`T_{in}`$ for layers with $`\mathrm{\Sigma }_{max}=10^9,10^{10}`$ and $`10^{11}`$ g cm<sup>-2</sup>. Applying the boundary condition at $`\mathrm{\Sigma }_{max}`$ rather than deeper down is an approximation, but the error due to this is not serious. For the high surface densities we consider, the heat transfer is dominated by conduction (Heyl & Hernquist 2001), so the temperature gradient for $`\mathrm{\Sigma }>\mathrm{\Sigma }_{max}`$ is small. Moreover, we examine models for several values of $`T_{in}`$ which further mitigates any error.

The calculations proceed in two stages. First, we solve for the steady state profile of the accretion layer by making the replacement $`d/dt\dot{\mathrm{\Sigma }}d/d\mathrm{\Sigma }`$, $`/\mathrm{\Sigma }d/d\mathrm{\Sigma }`$ in the governing equations. This gives five ordinary differential equations with four outer boundary conditions and one inner boundary condition. The solution gives the profiles of the basic fluid quantities: $`X(\mathrm{\Sigma }),Y(\mathrm{\Sigma }),P(\mathrm{\Sigma })=g\mathrm{\Sigma },\rho (\mathrm{\Sigma }),T(\mathrm{\Sigma })`$.

Having calculated the steady state structure of the accretion layer, we next check its stability. Various local stability criteria have been discussed in the literature (e.g., Bildsten 1998), in which one considers the properties of the gas at a single depth. FL87 proposed a global criterion involving an integral over the entire layer. While an improvement, this approach is still unsatisfactory since the authors assumed a constant temperature perturbation throughout the layer.

We have carried out a full linear stability analysis of the accretion layer. We start with the steady state solution and assume that it is slightly perturbed, $`Q(\mathrm{\Sigma })Q(\mathrm{\Sigma })+Q^{}(\mathrm{\Sigma })\mathrm{exp}(\gamma t)`$, where $`Q`$ corresponds to each of our five variables, and the perturbations $`Q^{}(\mathrm{\Sigma })`$ are taken to be small. We linearize the five equations described earlier (three of which include time derivatives), apply the boundary conditions, and solve for the eigenvalue $`\gamma `$. We obtain a large number of solutions for $`\gamma `$ (technically, there is an infinite number since we are dealing with continuous functions), many of which have both a real and imaginary part. We consider the accretion layer to be unstable if any eigenvalue $`\gamma `$ has a real part (growth rate) greater than the characteristic accretion rate $`\gamma _{acc}=\dot{\mathrm{\Sigma }}/\mathrm{\Sigma }_{max}`$.

If the steady-state model is unstable according to the above criterion, it cannot accrete matter stably with the particular $`\mathrm{\Sigma }_{max}`$ and $`\dot{\mathrm{\Sigma }}`$. Whether this instability manifests itself as a Type I burst depends on how the burning flame, once ignited, envelopes the surface of the star. This is as yet an unsolved problem, although Spitkovsky, Levin & Ushomirsky (2002) present a possible solution. In the following we assume that the instability will grow rapidly to the nonlinear regime and that the system will indeed exhibit Type I bursts.

## 3 Results

Figure 1 shows results for solar composition material ($`X_0=0.7`$, $`Y_0=0.27`$, $`Z_0=0.03`$) accreting on two kinds of compact stars.

The top three panels correspond to an accreting NS of mass $`1.4M_{}`$. We consider a range of accretion rates, parameterized by the ratio $`L_{acc}/L_{Edd}`$, where we take the Eddington luminosity to be $`L_{Edd}=4\pi GMc/\kappa _{es}`$ with $`\kappa _{es}=0.4\mathrm{cm}^2\mathrm{g}^1`$. We also consider a range of radii for the NS: $`\mathrm{log}(R/R_S)=0.20.6`$ ($`R=6.516`$ km). For each choice of $`L_{acc}/L_{Edd}`$ and $`R/R_S`$, we try three values of the surface mass density of the accreted layer: $`\mathrm{\Sigma }_{max}=10^9,10^{10},10^{11}\mathrm{g}\mathrm{cm}^2`$. If any of the three cases is unstable, i.e., if it has any eigenvalue with $`\mathrm{Re}(\gamma )>\gamma _{acc}`$, then we claim that the system will exhibit instabilities that may grow into Type I bursts.

The results presented correspond to three choices of the temperature at the base of the accreted layer: $`T_{in}=10^{8.5},10^8`$ and $`10^{7.5}`$ K from left to right. In this Letter we are primarily interested in transient LMXBs. Since these sources have very low luminosities in quiescence ($`L_X<10^{33}\mathrm{erg}\mathrm{s}^1`$), the core temperatures of the NSs are likely to be $`T_{in}10^8`$ K (Possenti et al. 2001; Brown, Bildsten & Chang 2002).

The calculations shown in Fig. 1 indicate that NSs are unstable to bursts for a wide range of $`T_{in}`$, but that the width of the instability strip (as a function of $`L_{acc}`$) is less for higher values of $`T_{in}`$. The reason for this is clear from the analysis of P82 who showed that when the flux escaping from the stellar core into the accretion layer increases (which happens when $`T_{in}`$ increases), bursting behavior is restricted to a smaller range of $`\dot{\mathrm{\Sigma }}`$.

In the case of transient NS LMXBs, we expect $`T_{in}`$ to be $`10^8`$ K. Fig. 1 indicates that these sources should be unstable to bursts for accretion luminosities up to $`L_{acc}0.3L_{Edd}`$. The predicted luminosity limit is generally consistent with observations; the only NSs that are known not to burst are bright Z sources with $`L_{acc}L_{Edd}`$ (e.g., Matsuba et al. 1995), and X-ray pulsars. Although the latter are significantly less luminous than Eddington, they accrete effectively at close to the Eddington rate since the accreting matter is channeled onto a small area on the NS surface by strong magnetic fields (see Lamb 2000 for a detailed discussion of this argument).

Below $`L_{acc}0.3L_{Edd}`$, the instability is initially of a mixed form in which a He-burning instability triggers a burst in which H and He both burn explosively. At lower luminosities, nearly all the H is burned steadily and the instability corresponds to a pure He burst. These results are consistent with previous work (e.g., Bildsten 1998). The calculations indicate that bursting behavior cuts off below an accretion luminosity $`L_{acc}10^{1.5}L_{Edd}`$. The cutoff is the result of the restriction $`\mathrm{\Sigma }_{max}10^{11}\mathrm{g}\mathrm{cm}^2`$ in our models. Systems with luminosities below our cutoff are still unstable to He bursts, but only at extremely high column density $`\mathrm{\Sigma }_{max}`$. We have ignored such bursts since the recurrence time $`t_{rec}`$ is too long to be of interest for observations of transient X-ray binaries.

For $`T_{in}=10^{7.5}`$ K, there is a second instability strip at low luminosities $`L_{acc}10^2L_{Edd}`$. This strip corresponds to pure H bursts; the possibility of such bursts was first noted by Fujimoto et al. (1987). An interesting difference between the two instability strips is that the strip on the right generally has complex eigenvalues $`\gamma `$ for the unstable modes while that on the left has real eigenvalues. As explained above, the gap between the two strips is real if we restrict ourselves to reasonable values of $`\mathrm{\Sigma }_{max}`$ and $`t_{rec}`$. We are not aware of clear observational evidence for or against the gap, but such evidence could be searched for in future observations.

The bottom three panels in Fig. 1 show results for a $`10M_{}`$ BH candidate with a surface. The three panels correspond to different choices of the base temperature: $`T_{in}=10^{7.5},10^7`$ and $`10^{6.5}`$ K from left to right. The particular choices of $`T_{in}`$ are motivated by the extraordinarily low quiescent luminosities of many transient BH LMXBs ($`L_X<10^{31}\mathrm{erg}\mathrm{s}^1`$, Narayan et al. 2001; a few BH LMXBs are brighter than this limit, but even these are not likely to have $`T_{in}10^{7.5}`$ K). We consider stellar radii in the range $`\mathrm{log}(R/R_S)=0.050.45`$, corresponding to $`R`$ between $`(9/8)R_S`$ and about $`3R_S`$, a factor of nearly twenty in the surface gravity. The choice $`(9/8)R_S`$ corresponds to the smallest radius in general relativity for an object whose density either decreases or remains level with increasing radius (ST83).

The calculations indicate that BH candidates with surfaces are not very different from NSs in their bursting behavior. Except for a modest rightward shift of the positions of the instability zones, the results in the lower three panels are quite similar to those for a NS with a similarly low value of $`T_{in}`$ (upper right panel). As in the case of NSs, bursts are expected for BH candidates even for larger values of $`T_{in}`$ than we have considered, except that the instability strip becomes narrower when $`T_{in}`$ exceeds $`10^8`$ K. We conclude that BH candidates are as prone to the instabilities that lead to Type I bursts as NSs are. The absence of bursts in BH LMXBs is thus highly significant and argues for the lack of surfaces in these systems.

Although our analysis is for a spherical, non-rotating $`10M_{}`$ star, the results are not expected to differ greatly for a rotating object. The effective surface gravity and accretion rate vary by only a factor of 2–3 from pole to equator even for a maximally rotating oblate ellipsoid that is on the verge of the secular triaxial instability; furthermore, the dependence on compactness and equation of state is small (see Gondek-Rosinska & Gourgoulhon 2002; ST83). In comparison, Fig. 1 shows that bursts are present for a wide range of surface gravity $`g`$ (factor of 20) and accretion rate $`\dot{\mathrm{\Sigma }}`$ (factor of 100). The rather modest variation of $`g`$ and $`\dot{\mathrm{\Sigma }}`$ with latitude in a rotating star is thus not likely to have an effect.

Because our model focuses only on the most important physical effects, and neglects many details, the exact positions of the instability strips in Fig. 1 may be uncertain at the level of say a factor of two in accretion luminosity. We believe, however, that the overall pattern of instability we have computed is fairly robust. Also, the calculations we present here do not directly predict what kind of bursts are produced. We leave this for a more detailed paper, but briefly, we find that for NSs with $`L_{acc}0.1L_{Edd}`$, bursts have durations of a few seconds and $`t_{rec}10`$ hours. For lower $`L_{acc}`$, as is well-understood, bursts have longer $`t_{rec}`$ and larger fluences. These results agree qualitatively with observations.

## 4 Discussion

It is clear from the theory of bursts (Bildsten 1998, and references therein) that bursting behavior is largely determined by the surface gravity $`g`$, the mass accretion rate, and the composition of the accreting material. Since these parameters are similar in NS and BH systems, the bursting behavior of the two should be similar. We have quantified this argument with a global linear stability analysis which represents a technical advance over previous calculations.

The results presented in Fig. 1 show that if BH candidates had surfaces they ought to experience thermonuclear instabilities as commonly as NSs do; by inference, they should have frequent Type I bursts. However, no BH candidate has exhibited Type I bursts. The most obvious explanation is that NSs have surfaces and BH candidates have event horizons (Menou 2001). If there is no surface, material cannot accumulate, and therefore cannot become hot or dense enough to trigger a thermonuclear explosion.

Before we can claim that this “proves” the reality of the event horizon, more work is needed. First, we need to show that the model is able to reproduce the more detailed features of Type I bursts as observed in NS LMXBs. The statistics of burst durations and recurrence times (e.g., van Paradijs, Penninx & Lewin 1988) and the occurrence of oscillatory behavior in some systems (Revnivtsev et al. 2000) ought to appear naturally in a realistic model. Also, any NS systems that burst when they should not by our calculations, or do not burst when they should according to the model, need to be explained. Second, the role of the inner boundary condition needs to be explored in detail. In Fig. 1, we see that different choices of the base temperature for a BH candidate give similar results. We have tried other boundary conditions, and also tried changing the composition of the accreting gas. In all cases we find that the accumulating layer is unstable to bursts for a wide range of luminosity, both in $`1.4M_{}`$ NSs and $`10M_{}`$ compact stars with surfaces. Finally, the difficult issue of flame propagation over the surface of the star once the instability has been triggered needs to be addressed (e.g., Lamb 2000, Spitkovsky et al. 2002); the effect of rotation, for instance, is presently unclear.

We should caution that the analysis presented here assumes that the accumulating gas on the surface of a BH candidate behaves like normal matter, with nucleons and electrons; this is the case for standard models of strange stars and Q-stars (Glendenning 1997). Obviously, the argument becomes invalid if the properties of the gas change drastically, e.g., if the nuclei disappeared and were replaced by quarks. Whether such extreme changes are plausible remains to be seen. The density and pressure at the base of the bursting layer do not go above $`\mathrm{few}\times 10^8\mathrm{g}\mathrm{cm}^3`$ and $`10^{26}\mathrm{erg}\mathrm{cm}^3`$ even in the most extreme cases we have considered. It is hard to imagine exotic physics being important under these conditions (Glendenning 1997).

On the observational front, we should check whether some NSs that lie within the unstable regions of Fig. 1 are stable to bursts. Any obvious large-scale disagreement between the observed burst behavior of NS LMXBs and the results presented here would indicate that the model is missing important physics. In the case of BH binaries, we should use observational data to derive quantitative limits on bursting activity. The transient BH LMXBs are particularly important for such work since they vary over a wide range of $`L_{acc}`$ during their accretion outbursts (Tanaka & Shibazaki 1996). If these sources have surfaces it is virtually impossible to arrange matters such that the objects have no bursts at all over the entire range of $`L_{acc}`$. A firm demonstration that BH transients do not have Type I bursts would thus be a strong argument for the presence of event horizons in these systems.

We thank Andrew Cumming, Alex Ene, Kristen Menou, Bohdan Paczyński and Greg Ushomirsky for useful discussions, and two referees for useful comments. RN was supported in part by the W. M. Keck Foundation as a Keck Visiting Professor. RN’s research was supported by NSF grant AST-9820686 and NASA grant NAG5-10780. JSH was supported by the Chandra Postdoctoral Fellowship Award # PF0-10015 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-39073.