text:
# Fractals and dynamical chaos in a random 2D Lorentz gas with sinks

## 1 Introduction

Many recent works have been devoted to the problem of diffusion in the presence of static absorbing traps, or sinks. In such models, the diffusing particles can be seen as modeling some active chemical species or some kind of excitation. The trapping then consists in a reaction with another chemical species leading to a non-reactive product, or in an excitation trapping converting the energy to some other form . The important quantity to study is the survival probability $`P(c,t)`$, i.e., the probability for a diffusing particle to survive in the system at time $`t`$ for a trap concentration $`c`$. In the case of randomly distributed traps, in a system of infinite spatial extension, the survival probability $`P(c,t)`$ of the diffusing particles obeys asymptotically a stretched exponential law of the form

$$P(c,t)\mathrm{exp}\left[A_du^{2/(d+2)}t^{d/(d+2)}\right],$$
(1)
where $`d`$ is the dimension of the system, $`c`$ is the concentration of static traps, $`u=\mathrm{ln}(1c)`$ and $`A_d`$ is a constant depending on the dimension and on the geometry of the system. This behavior is due to the existence with a non-vanishing probability of very large trap-free regions in which particles can survive for a very long time .

The stretched exponential law (1) has been much studied in stochastic systems although, until now, no study of this phenomenon has been carried out in deterministic chaotic systems. Such a study is of special interest in deterministic chaotic systems where the trapping process is characterized not only by the behavior of the survival probability but also by the fractal set of all the trajectories that never meet a trap. This set is defined in the phase space of the deterministic system and is referred to as a fractal repeller .

The fractal repeller associated with the deterministic chaotic diffusion in the presence of periodically distributed traps or sinks has been studied in a previous paper . The system studied there was a two-dimensional periodic Lorentz gas: a point particle undergoes elastic collisions on hard disks fixed in the plane and forming a regular triangular lattice. Some of these disks, forming a regular superlattice over the disk lattice, play the role of traps: when the point particle collides on one of them, it is annihilated. In this case, since the size of trap-free regions is bounded, the decay is exponential, characterized by an escape rate $`\gamma `$, so

$$P(c,t)e^{\gamma t}.$$
(2)
The dynamics is chaotic and the repeller has a non-integer fractal dimension.

The purpose of the present paper is to introduce a model for the stretched exponential law (1) which is based on a deterministic dynamics in order to be able to investigate the deterministic features of the process and, especially, the properties of the associated fractal repeller. We might suspect that the fractal dimension of the repeller is integer and equal to the phase space dimension $`D_H=3`$ since, in a random system of infinite spatial extension, the survival probability (1) decays more slowly than the exponential (2). However, if the dimension is integer, a finer characterization of the repeller is required and we might even wonder if the repeller is still a fractal.

In order to find a solution to this puzzle, we first apply the escape-rate formalism to finite two-dimensional random Lorentz gases. We confirm that, in such systems, the asymptotic decay is exponential and that the corresponding escape rate can be related to the positive Lyapunov exponent and to the fractal dimension of the fractal repeller.

We next consider infinite 2D random Lorentz gases. Thanks to our fast event-driven numerical algorithm and appropriate simulation method, we are able to extend the study of the survival probability to these deterministic chaotic systems.

Finally, we study the transition from finite to infinite 2D random Lorentz gases, which then allows us to determine the evolution of the scaling properties of the fractal repeller when the system size increases and to solve the aformentioned puzzle.

The paper is organized as follows. Section 2 is devoted to finite two-dimensional random Lorentz gases with sinks. The asymptotic decay is shown to be exponential, described by an escape rate $`\gamma `$. Moreover, we show that this open chaotic system is characterized by the existence of a fractal repeller. In these open systems, the average Lyapunov exponent and the fractal dimensions of the repeller are calculated with a non-equilibrium probability measure. The validity of the relation (2) is confirmed for these finite size systems. The case of infinite two-dimensional random Lorentz gases is presented in section 3. Thanks to our appropriate simulation method, the stretched exponential decay is observed for the first time in a deterministic chaotic system. The transition from finite systems with exponential decay to infinite ones obeying stretched exponential laws is studied in section 4. For a given concentration of traps, the Hausdorff dimension of the repeller is shown to increase with system size, to finally reach an integer value equal to the dimension of phase space, as it would be the case for a closed system. We show that, nevertheless, the repeller remains fractal in the infinite-size limit and that it is no longer characterized by a simple power-law scaling as the classical fractals but by a more complicated scaling. Conclusions are drawn in section 5.

## 2 The two-dimensional random Lorentz gas of finite spatial extension

### 2.1 Description of the model

The 2D random Lorentz gas consists of a point particle colliding elastically on randomly and uniformly distributed non-overlapping hard disks fixed in the plane. In finite (or periodic) geometry we will require in addition that the horizon is finite, i.e., that the point particle, starting from any initial condition, collides on a disk after a finite time. Under these conditions the system is diffusive, with a well-defined diffusion coefficient. Some disks are chosen to be sinks: when the point particle collides on one of them, it is absorbed, according to the following reaction scheme:

$`\mathrm{X}+\mathrm{inert}\mathrm{disk}`$ $``$ $`\mathrm{X}+\mathrm{inert}\mathrm{disk},`$ (3)
$`\mathrm{X}+\mathrm{sink}`$ $``$ $`\mathrm{}+\mathrm{sink}.`$ (4)
In practice, we consider a square unit cell of size $`L\times L`$. In this cell, $`n_d`$ disks of unit radius are placed randomly in the following way. The disks are first placed periodically, on a square or triangular lattice. A random velocity of unit magnitude is assigned to each of them, and the system is evolved as a hard ball system during a relaxation time $`T_{relax}`$. The configuration is then frozen and $`n_s`$ disks are randomly chosen to be sinks. Each disk has the probability $`c`$ to be a sink, independently of the others, so $`n_s`$ is always close to $`cn_d`$. Three examples of such unit cells are shown in Fig. 1.

The motion of the point particle in the system is defined under periodic boundary conditions. The collisions on the disks are elastic so that the energy of the system is conserved. The magnitude of the point particle velocity is therefore a constant, fixed here to one. The particle velocity can be specified by a unique coordinate, for example the angle $`\psi `$ between the velocity and the $`x`$-axis. Its position is specified by two coordinates $`\{x,y\}`$. The three phase-space variables of the point particle are therefore $`\mathrm{\Gamma }=\{x,y,\psi \}`$.

### 2.2 Escape and fractal repeller

#### 2.2.1 Asymptotic behavior

To study the asymptotic time behavior in such a system, we define the survival probability $`P(c,t)`$, which is the probability for a point particle to be still in the system at time $`t`$, for a sink concentration $`c`$ now taken to be exactly $`n_s/n_d`$. For finite systems in fact this survival probability depends on the positions of the scatterers as well as on the distribution of the traps, but as system size increases the fluctuations of $`P(c,t)`$ around its average value will become increasingly smaller. Suppose the long time behavior of $`P(c,t)`$ can be written in the form

$$P(c,t)\mathrm{exp}\left[Ct^{\alpha (t)}\right].$$
(5)
In the case of an exponential decay, $`\alpha `$ would be equal to $`1`$, and for the stretched exponential behavior we expect to hold asymptotically in two dimensions, it would be equal to $`0.5`$. This may be generalized into a time dependent $`\alpha (t)`$ defined as

$$\alpha (t)=\frac{d}{d\mathrm{ln}t}\mathrm{ln}[\mathrm{ln}P(c,t)].$$
(6)
This quantity, obtained numerically, is plotted versus time in Fig. 2.

The first graph, in Fig. 2a, represents the values of $`\alpha `$ for the disk configuration represented in Fig. 1a, but with only 4 traps, corresponding to a trap density $`c0.05`$. An average has been taken over 10 different trap configurations. Figure 2b (respectively 2c) presents the results obtained for the disk configuration of Fig. 1b (respectively 1c), with 5 (respectively 6) traps, corresponding to $`c0.05`$. In all three cases, we observe that, after a transient, $`\alpha `$ slowly tends to 1 .

The escape rate characterizing the exponential decay may be obtained as follows . Let us take $`N_0`$ initial conditions randomly distributed in the phase space of variables $`\mathrm{\Gamma }=(x,y,\psi )`$ according to some initial measure $`\nu _0`$. After a time $`t`$, only a number $`N_t`$ of particles will remain in the system. This number will decrease monotonically with $`t`$. The set of particles remaining in the random Lorentz gas at time $`t`$ consists of those colliding on a trap, and therefore escaping, at an escape time $`𝒯(\mathrm{\Gamma }_0)`$, which is larger than $`t`$:

$$\mathrm{{\rm Y}}(t)=\{\mathrm{\Gamma }_0:t<𝒯(\mathrm{\Gamma }_0)\}.$$
(7)
The decay of the number of particles is expected to be asymptotically exponential due to the finite size of the system. The escape rate is thus defined as

$$\gamma =\underset{t+\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}\frac{N_t}{N_0}=\underset{t+\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}\nu _0[\mathrm{{\rm Y}}(t)].$$
(8)
The escape rate is a characteristic quantity of the system, independent of the choice of the initial measure $`\nu _0`$ as long as this latter is smooth enough. The escape rate is easily accessible by numerical computations, as shown in Fig. 3. The logarithm of the survival probability is plotted versus time in Fig. 3 for the disk and trap configurations presented in Figs. 1a and c. The behavior is clearly exponential after a short transient. The slopes of the curves give escape rates of $`\gamma =0.025`$ and $`\gamma =0.017`$ respectively.

#### 2.2.2 Fractal repeller and nonequilibrium probability measure

A two-dimensional random Lorentz gas of finite spatial extension with traps is an open system, obeying an asymptotically exponential decay. Moreover, its dynamics is chaotic. It can therefore be studied by the methods developed in Refs. , in order to relate its microscopic chaotic dynamics to its macroscopic exponential decay.

An important property of such a system is the existence of a fractal repeller. For a typical initial condition $`\mathrm{\Gamma }_0`$, the particle will collide on a trap and escape after a finite escape time $`𝒯(\mathrm{\Gamma }_0)`$. The same will happen if we follow the trajectory in the negative time direction. Although a forward and a backward escape time exist for almost all trajectories, there do exist trajectories that never collide on a sink and remain trapped forever in the Lorentz gas in both time directions. They can be both periodic and non-periodic. Because of the defocusing character of the collisions on the disks, these non-colliding trajectories are unstable and form a fractal set of zero Lebesgue measure in phase space: this set is called the fractal repeller $``$. It is a typical feature of chaotic-scattering processes.

Evidence of the fractal character of this repeller is given by the forward escape time $`𝒯(\mathrm{\Gamma }_0)`$ as a function of the initial condition, as shown in Fig. 4. This function is finite for almost all initial conditions since the particle typically collides with a sink after a finite time. However, it is infinite for trajectories trapped on the fractal repeller. These trajectories start from initial points on the stable manifold (in the forward direction) of a trajectory that is trapped in both time directions. The infinities of the escape-time function thus establish a fractal set formed by the stable manifold of the fractal repeller, $`W_\mathrm{s}(_n)`$.

Fig. 4b-d shows this escape time as a function of the initial position, given by an angle $`\theta _0`$ around the disk marked by a star in Fig. 4a. The velocity is chosen normal to this disk. Due to the random distribution of disks, the fractal set is not quasi self-similar, as it was the case in Ref. . However, Figs. 4c-d show that a fine structure of escape windows exists at smaller and smaller scale. The characteristics of this structure, such as correlation functions and the distribution of sizes of empty windows, may be expected to show self-similarity on small length scales.

The fractal repeller is the support of an invariant probability measure which can be constructed by considering statistical averages over all the trajectories which have not yet escaped during a time interval $`[T/2,+T/2]`$ and, thereafter, taking the limit $`T\mathrm{}`$. In this limit the trajectories which do not escape during $`[0,+T/2]`$ have initial conditions on the stable manifold of the repeller, $`W_\mathrm{s}()`$. On the other hand, the initial conditions of the trajectories which, in the reverse time direction do not escape during $`[0,T/2]`$ approach the unstable manifold, $`W_\mathrm{u}()`$. In the limit $`T\mathrm{}`$, imposing no escape over the whole interval $`[T/2,+T/2]`$ selects trajectories which approach closer and closer the repeller, given by the intersection $`W_\mathrm{s}()W_\mathrm{u}()=`$. Statistical averages over these selected trajectories define an invariant measure having the fractal repeller for support.

In systems with time-reversal symmetric collision dynamics, such as the present one, statistical averages calculated over the aforementioned invariant measure are equivalent to statistical averages over a conditionally invariant measure defined by selecting the trajectories which do not escape during $`[0,+T/2]`$ only and taking the limit $`T\mathrm{}`$. This conditionally invariant measure has been constructed in Refs.  and we report the reader to these papers for further detail. Let us here add that, in the limit of low densities of sinks, the system tends to the closed finite Lorentz gas, whereupon the invariant and the conditionally invariant measures reduce to the microcanonical equilibrium measure. In this limit, the escape process stops and the fractal repeller fills the whole phase space.

#### 2.2.3 Average Lyapunov exponent of the repeller

The positive Lyapunov of the repeller can be calculated by averaging over the $`N_T`$ particles still in the system after time $`T`$ according to

$$\lambda =\underset{T\mathrm{}}{lim}\underset{N_0\mathrm{}}{lim}\frac{1}{T}\frac{1}{N_T}\underset{j=1}{\overset{N_T}{}}\mathrm{ln}\mathrm{\Lambda }_T(\mathrm{\Gamma }_0^{(j)}),$$
(9)
where $`\mathrm{\Lambda }_T`$ at time $`T`$ is the stretching factor of an initial perturbation on the trajectory $`\mathrm{\Gamma }_0^{(j)}`$ of surviving particles . The average (9) is equivalent to an average based on the aforementioned conditional invariant measure as shown in Refs. .

Numerically, the positive Lyapunov exponent is obtained as the average slope of the logarithm of the stretching factor as a function of time according to Eq. (9). This is illustrated in Figs. 5a and 6a, for the disk and trap configurations depicted respectively in Figs. 1a and 1c. The average Lyapunov exponent respectively equals $`\lambda =0.75`$ and $`\lambda =1.67`$.

#### 2.2.4 Fractal dimensions of the repeller

Important characteristics of the fractal properties of a repeller are its generalized fractal dimensions. Repellers most often are multifractals, with nontrivial generalized fractal dimensions $`D_q`$ . Here the embedding dimension of the repeller is three, since we are working in a 3-dimensional phase space. Therefore, the fractal dimensions of the repeller must satisfy the inequality $`0D_q<3`$. Instead of calculating directly the dimensions $`D_q`$ of the fractal repeller of the full flow we will consider a set of so-called partial fractal dimensions $`d_q`$ defined for the intersection of the stable manifold of the repeller $``$ with a line $``$ across this stable manifold. Like $``$ this mostly is a multifractal set. Obviously the values of $`d_q`$ must be between 0 and 1. Because the system is time-reversal symmetric the partial dimension in the stable direction is equal to the one in the unstable direction. The direction of the flow contributes a partial dimension equal to unity. Hence for given $`q`$, the dimension $`D_q`$ is related to the partial dimension $`d_q`$ by $`D_q=2d_q+1`$.

In the present paper, we focus on the fractal dimensions for $`q=0`$ and $`q=1`$. For $`q=0`$, we have the Hausdorff dimension commonly defined by

$$d_\mathrm{H}=d_0=\underset{\epsilon 0}{lim}\frac{\mathrm{ln}N(\epsilon )}{\mathrm{ln}\epsilon },$$
(10)
where $`N(\epsilon )`$ is the number of intervals $`I_j`$ of equal size $`\epsilon `$, required to cover the fractal $`f`$ defined on a line. For $`q=1`$, we find the so-called information dimension

$$d_\mathrm{I}=d_1=\underset{\epsilon 0}{lim}\frac{1}{\mathrm{ln}\epsilon }\underset{j=1}{\overset{N_l}{}}p_j(\epsilon )\mathrm{ln}p_j(\epsilon ),$$
(11)
where $`p_j(\epsilon )=\mu (I_jf)`$ is the probability weight of the interval $`I_j`$ of size $`\epsilon `$. Here, $`\mu `$ denotes the probability measure describing the distribution of points in the fractal set $`f`$. The Hausdorff and information codimensions are defined respectively as $`c_H=1d_H`$ and $`c_I=1d_I`$.

For fractal repellers of an open hyperbolic system, a relation exists between the partial information codimension $`c_I`$, the (largest) positive Lyapunov exponent $`\lambda `$, and the escape rate $`\gamma `$

$$\gamma =\lambda c_I.$$
(12)
To a good approximation the information codimension can be replaced by the Hausdorff codimension . This last dimension may be calculated with the aid of a numerical algorithm developed by the group of Maryland . The basic idea of this algorithm is to consider an ensemble of pairs of initial conditions, separated by $`\epsilon `$, along the line $``$ defined above. A pair is said to be uncertain if there is a discontinuity of the escape-time function between both initial conditions. On the other hand, when the pair is certain, the initial conditions belong to an interval of continuity of this function. An estimation of the fraction of uncertain pairs may be obtained in the following way: typically for an uncertain pair the time after which the pair becomes uncertain coincides with the time that the distance between the points becomes of order unity. This time may thus be estimated as $`t_u=(1/\lambda )\mathrm{ln}(1/\epsilon )`$. Since trajectories escape at the rate $`\gamma `$, the fraction of uncertain pairs would be proportional to $`\mathrm{exp}(\gamma t_u)\epsilon ^{\gamma /\lambda }=\epsilon ^{c_I}`$. However, such an estimation needs to be refined to take into account the multifractal character of the repeller and the possibility that $`c_Ic_H`$ due to local fluctuations of the Lyapunov exponents. A refined estimation of the fraction of uncertain pairs given in the Appendix A for the simple case of a piecewise-linear one-dimensional map shows that

$$f(\epsilon )\epsilon ^{c_\mathrm{H}},$$
(13)
as proved for the general case in Ref. . We notice that, to leading order in the sink density $`c`$, the Hausdorff codimension equals the information codimension .

In the case of the Lorentz gas, discontinuities in the escape-time function are not very practical to calculate numerically. Instead one may calculate the fraction of pairs that do not collide with exactly the same sequence of disks until escape from the system. For small $`\epsilon `$ this will behave in precisely the same way as the fraction of uncertain pairs, because it is determined by the same criterion of reaching a distance of order unity between both trajectories of a pair before escaping. The details of the algorithm used here are described in Refs. .

In Figs. 5b and 6b, the logarithm of the fraction of uncertain pairs (in fact pairs with different collision sequences before escape) is plotted as a function of the logarithm of the small $`\epsilon `$ separating the two initial conditions of a pair, for the systems defined respectively by the unit cells depicted in Figs 1a and 1c. After a transient regime, the linear behavior confirms the power-law prediction of Eq. (13) . The slope gives a value of the Hausdorff codimension equal to $`c_H=0.034`$ and $`c_H=0.01`$ respectively.

For both configurations, the value of the product $`\lambda c_H`$ is in good agreement with the value of the corresponding escape rate: $`\lambda c_H=0.025=\gamma `$ in the former case and $`\lambda c_H=0.017=\gamma `$ in the latter. Since $`c_H=c_I`$ to leading order in the sink density $`c`$, this confirms the validity of the relation (12) for random Lorentz gases of finite spatial extension. This expression connects directly two quantities characterizing the microscopic chaotic dynamics of the Lorentz gas, the average Lyapunov exponent $`\lambda `$ and the Hausdorff codimension $`c_H`$ of the fractal repeller, and a quantity describing the macroscopic transport behavior of this open system, the escape rate $`\gamma `$.

## 3 The two-dimensional random Lorentz gas of infinite spatial extension

### 3.1 Description of the model and simulation method

We now consider a two-dimensional random Lorentz gas of infinite spatial extension. This supposes that the sinks are randomly distributed among an infinite number of disks. As a consequence, the density of larger and larger trap-free regions decays exponentially with their volume but does not vanish. A suitable method is needed to simulate this situation. It cannot be done with the simple method presented in the previous section: if, in a possibly large, but in practice always finite system, the nature of the disks is fixed a priori, the size of the largest trap-free region is determined. Trap-free regions of a size that is very rare for the given system size in practice just will not be observed. However, a certain probability $`c`$ to be a sink can be assigned to each disk . This leads to an average over the different possible sink configurations, as a result of which the probabilities for finding large trap-free regions will attain just the proper values.

Let us define the corresponding configuration-averaged survival probability $`P(c,t)`$, where $`c`$ may be interpreted alternatively as the average concentration of sinks. The time evolution of $`P(c,t)`$ is given by

$$P(c,t)=\underset{n}{}(1c)^np(n,t)=(1c)^{n_t}$$
(14)
where $`n=n_t`$ is the number of different disks on which the point particle has collided during the time $`t`$ and $`p(n,t)`$ is the probability that this number is equal to $`n`$. $`\mathrm{}`$ is an ensemble average. $`(1c)`$ is the probability that a disk is not a sink, so that $`(1c)^n`$ is the probability that none of the disks met at time $`t`$ is a sink, in other words the probability that the point particle survives.

To compute the probability $`p(n,t)`$ to meet $`n`$ different disks in a time $`t`$, a purely deterministic simulation is used. The two-dimensional random Lorentz gas in which the point particle diffuses is defined as in the previous section, by a square unit cell of size $`L\times L`$. Periodic boundary conditions are applied. However, an index is associated with each square cell, allowing us to know which one is currently visited. This is necessary to distinguish one disk from its copies in the other cells. Indeed, if the point particle meets a disk and, later, one of its copies, this should be considered as collisions with different disks. In other words, although the scatterers are located on a periodic lattice, be it one with a very large unit cell of random filling, the distribution of traps has no periodicity whatsoever.

### 3.2 Theoretical predictions and simulation results

The asymptotic time dependence of the survival probability $`P(c,t)`$ is expected to be of the form (1) . For checking this dependence in our system, the easiest method using the idea of Eq. (14) is a direct simulation, in which $`p(n,t)`$ is estimated as an average over a large number of uniformly distributed initial conditions. The system considered here has as unit cell the one depicted in Fig. 1c, of size $`25\times 25`$ and with 127 disks. The number of initial conditions is of the order of $`10^8`$.

We present our data following a procedure proposed by Barkema et al. in Ref. . At short times and small concentrations, Eq. (14) can be approximated by

$$P(c,t)=(1c)^{n_t}(1c)^{n_t}.$$
(15)
This is the so-called Rosenstock approximation. In two dimensions, $`n`$ scales as $`4\pi \rho Dt/\mathrm{ln}t`$, with $`\rho `$ the average density of scatterers. We therefore expect a time dependence of the survival probability of the form:

$`\mathrm{ln}\left[P(c,t)\right]`$ $``$ $`4\pi \rho Dut/\mathrm{ln}t\text{at short times and small concentrations},`$
$`\mathrm{ln}\left[P(c,t)\right]`$ $``$ $`(ut)^{1/2}\text{asymptotically},`$ (16)
with again $`u=\mathrm{ln}(1c)`$. An efficient way to visualize the crossover from one regime to the other is to plot $`\mathrm{ln}\left[P(c,t)\right]/\mathrm{ln}t`$ as a function of $`\sqrt{ut}/\mathrm{ln}t`$ . In the Rosenstock regime, we expect a quadratic dependence, while the stretched exponential regime should produce a linear behavior. The results are represented in Fig. 7 in a double-logarithmic plot. The data collapse for values of $`c`$ taken in the range $`[0.01,\mathrm{},0.09,0.1,\mathrm{}0.9]`$. The crossover from a slope 2 to a slope 1 is clear.

The intersection point of the two linear fits gives an estimation of the crossover time from the Rosenstock regime to another regime which ends with the Donsker-Varadhan asymptotics. Using this as reference point, we obtain the dependence on the trap density $`c`$ of the crossover time and of the crossover survival probability , as presented in Figs. 8a and b respectively. These results show that the crossover occurs at shorter and shorter times and at larger and larger probabilities as the trap density increases. We observe that the crossover probability changes from $`𝒪(10^{20})`$ to $`𝒪(10^5)`$ as the trap density only varies from $`𝒪(10^2)`$ to $`𝒪(1)`$.

A comment is here in order. The extensive numerical simulations of Ref.  show that the Donsker-Varadhan behavior is approached only very gradually. It is known to provide an exact asymptote for the logarithm of $`P(c,t)`$ including the prefactor of $`(ut)^{1/2}`$. If one inserts the relevant values of the parameters one find this asymptote in Fig. 7 about 20% above the asymptote obtained from the data fit. Given the limited time range over which our numerical data were obtained we may conclude that their behavior is fully compatible with the one observed for stochastic systems in Ref. .

## 4 Transition from finite to infinite random systems

We have considered finite random Lorentz gases with exponential behavior and infinite ones with stretched exponential behavior. How does the transition between these two regimes evolve with increasing system size $`L`$? And what happens to the fractal repeller?

To answer these questions, let us first briefly sketch the arguments of Grassberger and Procaccia for the infinite case . They consider diffusing particles distributed according to a density $`\rho _c(𝐫,t)`$. The average density $`n_c(𝐫,t)`$ is defined as

$$n_c(𝐫,t)=\rho _c(𝐫,t)_s$$
(17)
where $`\mathrm{}_s`$ is an average over all the realizations of trap distributions or equivalently over all choices of the origin from which all particles start diffusing initially. A perfectly absorbing sphere of large radius $`R`$ centered at the origin is added to the system. The value of the new density $`n_c^{}(𝐫,t)`$ obtained in the presence of this sphere is necessarily smaller than $`n_c(𝐫,t)`$

$$n_c(𝐫,t)n_c^{}(𝐫,t).$$
(18)
$`n_c^{}(𝐫,t)`$ depends on the distribution of traps inside the sphere. If $`\{S\}=\{S_1,S_2,\mathrm{},S_k\}`$ are the positions of the traps in the sphere, we have

$$n_c(𝐫,t)n_c^{}(𝐫,t)=\underset{\{S\}}{}n_c^{}(𝐫,t|\{S\})P(\{S\}).$$
(19)
The sum of the right-hand side can be bounded from below by considering only one term, corresponding to the situation of a sphere free of traps:

$$n_c(𝐫,t)>n_c^{}(𝐫,t|\{S\}=0)P(0).$$
(20)
A random distribution of sinks corresponds to a Poisson distribution. The probability of finding $`k`$ traps in a volume $`v`$ is given by

$$P_k=\frac{\overline{N}^k}{k!}\mathrm{exp}(\overline{N})$$
(21)
with average value $`\overline{N}=n_sv`$, where $`n_s`$ is the average density of sinks. The probability to have no trap in the sphere is thus given by $`P(0)=\mathrm{exp}(Cn_sR^d)`$, where $`d`$ is the dimension of the system and $`C`$ the volume of the unit ball in $`d`$ dimensions.

On the other hand, $`n_c^{}(𝐫,t|\{S\}=0)`$ is obtained as solution of a diffusion equation of coefficient $`𝒟_0`$ with absorbing boundary conditions on the border of a sphere of radius $`|𝐫|=R`$. Such a solution is known to decay asymptotically as the exponential $`\mathrm{exp}(\kappa _d\frac{𝒟_0}{R^2}t)`$, where $`\kappa _d`$ is a constant depending on the dimensionality of the system.

Therefore, Eq. (20) becomes

$$n_c(𝐫,t)>n_c^{}(𝐫,t|\{S\}=0)P(0)\mathrm{exp}\left(\kappa _d\frac{𝒟_0}{R^2}tCn_sR^d\right).$$
(22)
The right-hand side may be maximized as a function of $`R`$, with $`R`$ between $`0`$ and $`+\mathrm{}`$ in an infinite system. This gives us an optimal value for $`R`$

$$R=\left(\frac{2\kappa _d}{dC}\frac{𝒟_0}{n_s}t\right)^{1/(d+2)}=\alpha t^{1/(d+2)}.$$
(23)
Inserting this result in Eq. (22), one obtains

$$n_c(𝐫,t)>\mathrm{exp}\left[\beta t^{d/(d+2)}\right].$$
(24)
Let us now consider a finite system of linear size $`2L`$ with a fixed configuration of traps. As long as $`R=\alpha t^{1/d+2}\left(\frac{\mathrm{ln}L}{n_s}\right)^{1/d}`$, the optimal value of $`R`$ is found in a fraction of the volume described by (21) and the stretched exponential behavior is observed. However, for longer times, the maximal value of $`R`$ is proportional to $`l_0(\mathrm{ln}L/n_s)^{1/d}`$. Inserting $`R=l_0`$ in Eq. (22) leads to an exponential decay of the form $`\mathrm{exp}\left(\kappa _d\frac{𝒟_0}{l_0^2}t\right)`$. As $`L`$ increases, the time $`t^{}`$ at which the transition occurs increases.

A logarithmic dependence of $`t^{}`$ on $`L`$ in two dimensions has been proven in Ref. .

We now want to study the properties of the repeller. As a preparation, let us first consider the simple one-dimensional map with escape represented in Fig. 9. Starting from the full unit interval at time zero one finds after $`k`$ discrete time steps an approximate forward repeller consisting of $`N_k=2^k`$ pieces of sizes $`\epsilon _k=(1/3)^k`$. The number of pieces is equal to the exponential of the Kolmogorov-Sinai entropy times time, $`N_k=\mathrm{exp}(h_{KS}k)`$, and the linear size of these pieces is given by $`\epsilon _k=\mathrm{exp}(\lambda k)`$ with the Lyapunov exponent $`\lambda =\mathrm{ln}3`$. This last equality can be intuitively understood as follows: as time goes on, new structures appear on the repeller on smaller and smaller scales. Therefore, long times correspond to small spatial scales. The total length $`L_k`$ of the repeller at time step $`k`$ is the product of $`N_k`$ and $`\epsilon _k`$. It decreases exponentially with time as

$$L_k=N_k\epsilon _k=\mathrm{exp}(h_{KS}k)\mathrm{exp}(\lambda k)=\mathrm{exp}(\gamma k),$$
(25)
where $`\gamma =\lambda h_{KS}`$ is the escape rate.

The fractal (Hausdorff) dimension of the repeller is thus given by the well-known result

$$D_H=\underset{\epsilon 0}{lim}\frac{\mathrm{ln}N(\epsilon )}{\mathrm{ln}\left(\frac{1}{\epsilon }\right)}=\frac{h_{KS}}{\lambda }=\frac{\mathrm{ln}2}{\mathrm{ln}3},$$
(26)
in the illustrative example of Fig. 9. On the other hand, we know that the fraction of uncertain pairs $`f(\epsilon )`$, as defined for the Maryland algorithm, is proportional to the size $`\epsilon `$ to the power $`C_H`$, with $`C_H`$ the codimension of the repeller. From this we obtain

$$f(\epsilon )\epsilon ^{C_H}\epsilon ^{1D_H}\frac{\epsilon }{\epsilon ^{D_H}}\epsilon N(\epsilon ).$$
(27)
Let us now imagine a one-dimensional system in which the escape would obey a stretched exponential law. In this case, the remaining length of the repeller at time $`k`$ would obey

$$L_k=N_k\epsilon _k=\mathrm{exp}(ck^{1/2}).$$
(28)
¿From the relation $`\epsilon =\mathrm{exp}(k\lambda )`$, we can express the time as $`k=\frac{1}{\lambda }\mathrm{ln}\frac{1}{\epsilon }`$. We get

$$N(\epsilon )=\frac{L_k}{\epsilon }=\frac{\mathrm{exp}(ck^{1/2})}{\epsilon }=\frac{1}{\epsilon }\mathrm{exp}\left[c\left(\frac{1}{\lambda }\mathrm{ln}\frac{1}{\epsilon }\right)^{1/2}\right].$$
(29)
¿From Eq. (27), we have

$$f(\epsilon )\epsilon N(\epsilon )\mathrm{exp}\left[c\left(\frac{1}{\lambda }\mathrm{ln}\frac{1}{\epsilon }\right)^{1/2}\right].$$
(30)
This dependence of the fraction of uncertain pairs on $`\epsilon `$ for a system with stretched exponential behavior can be tested numerically for the two-dimensional infinite random Lorentz gas, with the help of the Maryland algorithm. The algorithm has to be slighty modified to simulate an infinite system. As mentioned before, the absorbing or non-absorbing nature of the disks is not fixed a priori, but each disk has a probability $`c`$ to be a sink. At each collision on a disk, the separation of the pair of particles considered is tested. If separation occurs, the fraction of uncertain pairs is increased by $`(1c)^n`$, where $`n`$ is the number of different disks that the pair has met before separating, so $`(1c)^n`$ is the probability that the particles have not escaped before separating. Simultaneously, the fraction of escaped pairs is increased by $`[1(1c)^n]`$. In Fig. 10a, the logarithm of the fraction of uncertain pairs, $`\mathrm{log}_{10}f(\epsilon )`$, is plotted as a function of the logarithm of the distance between the two initial conditions of the pair, $`\mathrm{log}_{10}\epsilon `$. In order to check that this curve can be fitted by the expression (30), we plot in Fig. 10b $`\mathrm{log}[\mathrm{log}_{10}f(\epsilon )]`$ versus $`\mathrm{log}(\mathrm{log}_{10}\epsilon )`$. We obtain a linear behavior with slope $`\frac{1}{2}`$, in excellent agreement with Eq. (30).

In the case of a system of finite spatial extension but with high density of disks, as represented in Fig. 1c, we see in Fig. 3b that the transient behavior persists until $`t50`$. According to the relation $`\epsilon =\mathrm{exp}(t\lambda )`$, this corresponds to $`\mathrm{log}_{10}\epsilon 34`$. This is consistent with the deviation from the linear curve which is observed in Fig. 6b for large values of $`\epsilon `$.

In two-dimensional systems like the Lorentz gases we consider here, the codimension appearing in the relation $`f(\epsilon )\epsilon ^{c_H}`$ is a partial codimension, associated with a partial dimension $`d_H`$. The total dimension is given by $`D_H=1+2d_H`$. The total number of small three-dimensional pieces needed to cover the fractal repeller is given by

$$N(\epsilon )\frac{1}{\epsilon }\left[\frac{f(\epsilon )}{\epsilon }\right]^2\frac{1}{\epsilon ^3}\mathrm{exp}\left[2c\left(\frac{1}{\lambda }\mathrm{ln}\frac{1}{\epsilon }\right)^{1/2}\right].$$
(31)
This type of dependence of $`N(\epsilon )`$ on the size $`\epsilon `$ leads to an integer value of the Hausdorff dimension $`D_H`$. Indeed,

$$D_H=\underset{\epsilon 0}{lim}\frac{\mathrm{ln}N(\epsilon )}{\mathrm{ln}\left(\frac{1}{\epsilon }\right)}=3\underset{\epsilon 0}{lim}2c\left(\frac{1}{\lambda }\frac{1}{\mathrm{ln}\frac{1}{\epsilon }}\right)^{\frac{1}{2}}=3.$$
(32)
This dimension does not contain the information concerning the correction to the power law dependence appearing in Eq. (31). We notice that corrections to a power law also appear for Brownian motion in the plane . In this case, the dependence essentially would take the form

$$N(\epsilon )\frac{1}{\epsilon ^2}\frac{1}{\mathrm{ln}\mathrm{ln}\frac{1}{\epsilon }}$$
(33)
leading to the integer Hausdorff dimension $`D_H=2`$.

When the size $`L`$ of a finite system with fixed traps increases, the time during which stretched exponential behavior is observed increases likewise, and the range of $`\epsilon `$ over which the fraction of uncertain pairs deviates from a power-law behavior in $`\epsilon `$. Correspondingly, as seen from Fig. 10, the slope in the linear regime of $`\mathrm{log}_{10}f(\epsilon )`$ versus $`\mathrm{log}_{10}\epsilon `$ becomes smaller and smaller, i.e., the partial Hausdorff codimension of the repeller will be smaller. This means that with increasing $`L`$ the Hausdorff dimension of the repeller progressively increases towards the asymptotic value $`D_H=3`$. The same holds for the information dimension, as in the dominant regions the trap density approaches zero and we saw before that both dimensions approach each other in this limit.

In the case of the two-dimensional random Lorentz gas of infinite spatial extension, the integer value $`D_H=3`$ of the dimension of the fractal repeller can be understood by intuitive arguments. Indeed, this dimension corresponds to the case of a closed system, without escape. In an infinite system, for a given time $`t`$, there exists with a finite probability a region free of traps of size $`R=\left(\frac{2\kappa _d}{dC}\frac{𝒟_0}{n_s}t\right)^{1/(d+2)}`$. The point particles in this region behave as particles in a closed system until they escape through the “boundary” at $`R`$. As time goes on, larger and larger trap-free regions have to be considered. Even if their density decays exponentially with their size, it does not vanish, so that these regions give the dominant contribution to the measure of the repeller.

Our results show that the repeller is indeed a fractal object in spite of its integer dimension $`D_H=3`$. This fractality is characterized by the complicated scaling behavior of the fraction (30) of uncertain pairs or of the number (31) of pieces needed to cover the fractal. These theoretical predictions are nicely confirmed by our numerical results.

## 5 Conclusions

We have studied in this paper several types of two-dimensional random Lorentz gases with sinks.

The first model considered is a two-dimensional random Lorentz gas with sinks, of finite spatial extension. In this system, the sinks are randomly distributed among the disks. Due to the finite size of the system, the size of the largest trap-free region is bounded. Therefore, the decay of the number of particles is asymptotically exponential. For this open system, we have shown the existence of a fractal repeller formed by the trajectories that never escape from the system. Using a non-equilibrium measure on the repeller, we have characterized its chaotic and fractal properties. We have obtained the escape rate $`\gamma `$, the positive Lyapunov exponent $`\lambda `$, and the fractal dimensions of the repeller. A relation between $`\gamma `$, $`\lambda `$, and the Hausdorff codimension $`c_H`$ of the fractal repeller has been numerically tested in the limit of low sink densities. This relation establishes a link between the macroscopic decay and the underlying microscopic chaos.

The second case considered is the two-dimensional random Lorentz gas with sinks, of infinite spatial extension. This situation can be simulated by giving to each disk a certain probability to be a sink, while distinguishing between different periodic images of the same disk. This leads to an average over the different sink configurations. The decay of the survival probability asymptotically obeys the stretched exponential law (1). Here this behavior has been studied numerically for the first time in a deterministic chaotic system.

For a dense 2D random Lorentz gas of finite size $`L`$, the transient stretched exponential regime persists during a time fixed in first approximation by the typical size of the largest trap-free region in the system \[see Eq. (23)\]. As $`L`$ increases, the transient regime is observed for longer and longer times. Simultaneously, the dimension of the repeller increases, reaching the integer value $`D_H=3`$ in the infinite-size limit. We have here been able to show that the repeller remains fractal in the limit and to answer the puzzle mentioned in the introduction by finding the nontrivial scaling property of this plain fractal. Its fractal property can be characterized equivalently by the fraction of uncertain pairs of trajectories starting from initial conditions separated by $`\epsilon `$ or by the number of cells of size $`\epsilon `$ needed to cover the fractal. Both scaling functions show a correction to a simple power law proving that the repeller remains fractal. The correction depends on the Lyapunov exponent as well as on the exponent $`1/2`$ of the stretched exponential law (1) in $`d=2`$. This nontrivial scaling is remarkably confirmed by our numerical results.

In conclusion, we have been able to show that, in a random deterministic chaotic system of infinite spatial extension with sinks, the repeller composed of all the trajectories that never meet a sink is a plain fractal with special scaling properties related to the stretched exponential decay of the survival probability.

## Acknowledgments

The authors thank Professor G. Nicolis for support and encouragement in this research. IC and PG are supported financially by the National Fund for Scientific Research (F. N. R. S. Belgium). H.v.B. also acknowledges support by the Mathematical physics program of FOM and NWO/GBE.

## Appendix A Direct proof that the uncertainty exponent is the Hausdorff codimension

Let us consider the piecewise-linear one-dimensional map:

$$x_{n+1}=\{\begin{array}{cc}\mathrm{\Lambda }_0x_n0<x_n<A\hfill & \\ \mathrm{\Lambda }_1(1x_n)A<x_n<1\hfill & \end{array}$$
(34)
We assign the symbol $`\omega =0`$ to the interval $`[0,A]`$ and the symbol $`\omega =1`$ to the interval $`[A,1]`$.

If $`\mathrm{\Lambda }_0,\mathrm{\Lambda }_1>2`$, almost all the trajectories escape after a finite number of iterations. The subinterval $`I_{\omega _0\omega _1\mathrm{}\omega _{n1}}`$ is defined as the set of initial conditions for the trajectories which escape in $`n`$ iterations and successively visit the intervals corresponding to the symbolic sequence $`\omega _0\omega _1\mathrm{}\omega _{n1}`$. These subintervals form a partition of the unit interval $`[0,1]`$:

$$_{\omega _0\omega _1\mathrm{}\omega _{n1}}I_{\omega _0\omega _1\mathrm{}\omega _{n1}}=[0,1]$$
(35)
The length of the subinterval $`I_{\omega _0\omega _1\mathrm{}\omega _{n1}}`$ is

$$\mathrm{}_{\omega _0\omega _1\mathrm{}\omega _{n1}}=\frac{K}{\mathrm{\Lambda }_{\omega _0}\mathrm{\Lambda }_{\omega _1}\mathrm{}\mathrm{\Lambda }_{\omega _{n1}}}$$
(36)
Since the sum of lengths should be equal to one

$$\underset{n=1}{\overset{\mathrm{}}{}}\underset{\omega _0\omega _1\mathrm{}\omega _{n1}}{}\mathrm{}_{\omega _0\omega _1\mathrm{}\omega _{n1}}=1$$
(37)
the constant $`K`$ is given by

$$K=\frac{1\mathrm{\Lambda }_0^1\mathrm{\Lambda }_1^1}{\mathrm{\Lambda }_0^1+\mathrm{\Lambda }_1^1}$$
(38)
and the previously introduced constant $`A`$ can now be fixed to the value:

$$A=\frac{1+K}{\mathrm{\Lambda }_0}.$$
(39)
A pair is uncertain if it falls in a subinterval which is smaller than $`\epsilon `$. Therefore, the fraction of uncertain pairs is equal to

$$f(\epsilon )=\underset{n=1}{\overset{\mathrm{}}{}}\underset{\omega _0\omega _1\mathrm{}\omega _{n1}}{}\mathrm{}_{\omega _0\omega _1\mathrm{}\omega _{n1}}|_{\mathrm{}_{\omega _0\omega _1\mathrm{}\omega _{n1}}<\epsilon }.$$
(40)
The sum can be reordered in terms of the numbers $`n_0`$ (resp. $`n_1`$) of symbol $`0`$ (resp. $`1`$) in the symbolic sequence $`\omega _0\omega _1\mathrm{}\omega _{n1}`$. We have that $`n_0+n_1=n`$. The fraction of uncertain pairs becomes

$$f(\epsilon )=\underset{n_0=0}{\overset{\mathrm{}}{}}\underset{n_1=0}{\overset{\mathrm{}}{}}\frac{(n_0+n_1)!}{n_0!n_1!}\frac{K}{\mathrm{\Lambda }_0^{n_0}\mathrm{\Lambda }_1^{n_1}}|_{n_0\mathrm{ln}\mathrm{\Lambda }_0+n_1\mathrm{ln}\mathrm{\Lambda }_1>\mathrm{ln}(K/\epsilon )}.$$
(41)
We use the steepest-descent method for this double sum. The domain of summation is the set of pairs of non-negative integers $`(n_0,n_1)`$ above the straight line $`n_0\mathrm{ln}\mathrm{\Lambda }_0+n_1\mathrm{ln}\mathrm{\Lambda }_1=\mathrm{ln}(K/\epsilon )`$ in the plane $`(n_0,n_1)`$. The function to be summed decreases exponentially to zero as the exponents $`n_0`$ and $`n_1`$ tend to infinity above the straight line. Therefore, the sum will be dominated by the maximum of the function on the straight line. If we introduce the alternative variables

$$x_0\frac{n_0}{\mathrm{ln}(K/\epsilon )}\text{and}x_1\frac{n_1}{\mathrm{ln}(K/\epsilon )}$$
(42)
we thus find

$$f(\epsilon )\epsilon \left(\frac{K}{\epsilon }\right)^D\epsilon ^{1D}$$
(43)
with

$$D=x_0\mathrm{ln}\frac{x_0+x_1}{x_0}+x_1\mathrm{ln}\frac{x_0+x_1}{x_1}$$
(44)
and the constraint

$$x_0\mathrm{ln}\mathrm{\Lambda }_0+x_1\mathrm{ln}\mathrm{\Lambda }_1=1.$$
(45)
The maximum of the function is given by the steepest-descent condition

$$\frac{dD}{dx_0}=\mathrm{ln}\frac{x_0+x_1}{x_0}+\frac{dx_1}{dx_0}\mathrm{ln}\frac{x_0+x_1}{x_1}=0.$$
(46)
With the constraint (45), we find that the steepest-descent condition reduces to

$$\mathrm{ln}\mathrm{\Lambda }_1\mathrm{ln}\frac{x_0+x_1}{x_0}=\mathrm{ln}\mathrm{\Lambda }_0\mathrm{ln}\frac{x_0+x_1}{x_1}.$$
(47)
We now show that the exponent $`D`$ is equal to the Hausdorff dimension. In the present piecewise-linear map, the Hausdorff dimension of the repeller must satisfy by definition the following condition

$$\frac{1}{\mathrm{\Lambda }_0^D}+\frac{1}{\mathrm{\Lambda }_1^D}=1.$$
(48)
Let us show that it is indeed the case. For this purpose, we first express the quantity $`D`$ in terms of $`x_0`$ and $`x_1`$ by Eq. (44). Then, we use the steepest-descent condition (47) and, thereafter, the constraint (45). The calculations are the following:

$`{\displaystyle \frac{1}{\mathrm{\Lambda }_0^D}}+{\displaystyle \frac{1}{\mathrm{\Lambda }_1^D}}`$ $`=`$ $`e^{x_0\mathrm{ln}\frac{x_0+x_1}{x_0}\mathrm{ln}\mathrm{\Lambda }_0x_1\mathrm{ln}\frac{x_0+x_1}{x_1}\mathrm{ln}\mathrm{\Lambda }_0}`$ (49)
$`+e^{x_0\mathrm{ln}\frac{x_0+x_1}{x_0}\mathrm{ln}\mathrm{\Lambda }_1x_1\mathrm{ln}\frac{x_0+x_1}{x_1}\mathrm{ln}\mathrm{\Lambda }_1}`$
$`=`$ $`e^{(x_0\mathrm{ln}\mathrm{\Lambda }_0+x_1\mathrm{ln}\mathrm{\Lambda }_1)\mathrm{ln}\frac{x_0+x_1}{x_0}}+e^{(x_0\mathrm{ln}\mathrm{\Lambda }_0+x_1\mathrm{ln}\mathrm{\Lambda }_1)\mathrm{ln}\frac{x_0+x_1}{x_1}}`$ (50)
$`=`$ $`e^{\mathrm{ln}\frac{x_0+x_1}{x_0}}+e^{\mathrm{ln}\frac{x_0+x_1}{x_1}}=1`$ (51)
Q.E.D.

The previous calculation provides a direct proof that the uncertainty exponent of the Maryland algorithm is the Hausdorff codimension in a piecewise-linear map with two branches. Consequently, the uncertainty exponent should be expected to differ from the information codimension on general multifractal repellers. Indeed, by using the example (34), we have numerically observed that the uncertainty exponent agrees with the Hausdorff codimension but differs from the information codimension if $`\mathrm{\Lambda }_0\mathrm{\Lambda }_1`$.

The direct proof presented here can be extended to piecewise-linear one-dimensional maps with more than two branches. In the case of $`m`$ branches, there would be $`m`$ variables $`n_i`$ or $`x_i`$ and $`(m1)`$ steepest-descent conditions such as (46) because of the constraint that the maximum point belongs to the hyperplane $`_{i=1}^mn_i\mathrm{ln}\mathrm{\Lambda }_i=\mathrm{ln}(K/\epsilon )`$. Otherwise, a Lagrange multiplier may be introduced, to deal with this constraint in the search of the maximum. In the general case, the argument of Ref.  provides the proof.


format: text

text:
# Quest for Black Holes and Superstring Excitations in Cosmic Ray Data

## I Introduction

A promising route towards reconciling the apparent mismatch of the fundamental scales of particle physics and gravity is to modify the short distance behavior of gravity at scales much larger than the Planck length. This can be accomplished in a straightforward manner Arkani-Hamed:1998rs ; Antoniadis:1998ig  if one assumes that the standard model (SM) fields are confined to a 4-dimensional world (corresponding to our apparent universe), while gravity lives in a higher dimensional space. One virtue of this assumption is that very large extra dimensions are allowed without conflicting with current experimental bounds Hewett:2002hv , leading to a fundamental Planck mass much lower than its effective 4-dimensional value, even as low as the electroweak scale. In particular, if the spacetime is taken as a direct product of a 4-dimensional spacetime and a flat spatial torus $`T^{D4}`$ (of common linear size $`2\pi r_c`$), one obtains a definite representation of this picture in which the effective 4-dimensional Planck scale, $`M_{\mathrm{Pl}}10^{19}`$ GeV, is related to the fundamental scale of gravity, $`M_D`$, according to $`M_{\mathrm{Pl}}^2=8\pi M_D^{D2}r_c^{D4}`$ Arkani-Hamed:1998rs .

One of the most startling predictions of scenarios of this sort is that microscopic black holes (BHs) could be observed in forthcoming colliders Banks:1999gd  and at existing cosmic ray detectors Feng:2001ib . In fact, the non-observation of cosmic neutrinos has already been used to set lower bounds on the fundamental scale of $`D`$-dimensional gravity Anchordoqui:2003jr . This work relied on classical models of BH production, which are valid only for sufficiently massive BHs, $`M_{\mathrm{BH}}M_D`$. For masses close to $`M_D`$, gravity becomes strong and the classical description can no longer be trusted. String theory provides the best hope for understanding the regime of strong quantum gravity, and in particular for computing cross sections at energies close to the Planck scale Dimopoulos:2001qe . In principle embedding TeV-scale gravity models in realistic string models might facilitate the calculation of cross sections for BHs (and string excitations) of masses close to $`M_D`$.

Our purpose here is to explore this possibility within a simple string theory. We begin by reviewing the standard classical and semiclassical pictures of BH production and decay, and the string $``$ BH correspondence principle. We then calculate the $`\nu N`$ cross section for production of excited string states within superstring theory, and find it to be larger than the BH production cross section Cheung:2002aq . This makes it plausible that the BH cross section can be used as a lower bound down to $`M_{\mathrm{BH}}=M_D`$, which would significantly strengthen existing lower bounds on $`M_D`$.

To be specific, we will consider embedding of a 10-dimensional low-energy scale gravity scenario within the context of SO(32) Type I superstring theory, where gauge and charged SM fields can be identified with open strings localized on a 3-brane and the gravitational sector consists of closed strings that propagate freely into the internal dimensions of the universe Antoniadis:1998ig . After compactification on $`T^6`$ down to four dimensions, $`M_{\mathrm{Pl}}`$ is related to the string scale, $`M_s`$, and the string coupling constant, $`g_s`$, by $`M_{\mathrm{Pl}}^2=(2\pi r_c)^6M_s^8/g_s^2.`$ Hereafter, $`D=10`$. Generalization to arbitrary number of dimensions is straightforward.

## II Black Holes and String Balls with Low-Scale Gravity

Analytic and numerical studies have revealed that gravitational collapse takes place at sufficiently high energies and small impact parameters, as conjectured years ago by Thorne Thorne:ji . In the case of 4-dimensional head-on collisions D'Eath:hb , as well as those with non-zero impact parameter Eardley:2002re , a horizon forms when and only when a mass is compacted into a hoop whose circumference in every direction is less than $`2\pi `$ times its Schwarzschild radius up to a factor of order 1. In the 10-dimensional scenario the Schwarzschild radius still characterizes the maximum impact parameter for horizon formation Yoshino:2001ik . In the course of collapse, a certain amount of energy is radiated in gravitational waves by the multipole moments of the incoming shock waves D'Eath:hb , leaving a fraction $`yM_{\text{BH}}/\sqrt{\widehat{s}}`$ available to Hawking evaporate Hawking:1975sw . Here, $`M_{\text{BH}}`$ is a lower bound on the final mass of the BH and $`\sqrt{\widehat{s}}`$ is the center-of-mass energy of the colliding particles, taken as partons. This ratio depends on the impact parameter of the collision, as well as on the dimensionality of space-time Yoshino:2002br .

Subsequent to formation, the BH proceeds to decay Chamblin:2003wg . The emission rate per degree of particle freedom $`i`$ of particles of spin $`s`$ with initial total energy between $`(Q,Q+dQ)`$ is found to be Han:2002yy

$$\frac{d\dot{N}_i}{dQ}=\frac{\sigma _s(Q,r_s)\mathrm{\Omega }_{d3}}{(d2)(2\pi )^{d1}}Q^{d2}\left[\mathrm{exp}\left(\frac{Q}{T_{\mathrm{BH}}}\right)(1)^{2s}\right]^1,$$
(1)
where $`T_{\mathrm{BH}}=7/(4\pi r_s)`$ is the BH temperature,

$$r_s(M_{\mathrm{BH}})=\frac{1}{M_{10}}\left[\frac{M_{\text{BH}}}{M_{10}}\mathrm{\hspace{0.17em}\hspace{0.17em}8}\pi ^{3/2}\mathrm{\Gamma }\left(9/2\right)\right]^{1/7}$$
(2)
is the Schwarzschild radius Myers:un ,

$$\mathrm{\Omega }_{d3}=\frac{2\pi ^{(d2)/2}}{\mathrm{\Gamma }[(d2)/2]}$$
(3)
is the volume of a unit $`(d3)`$-sphere, and $`\sigma _s(Q,r_s)`$ is the absorption coefficient (a.k.a. the greybody factor). Recall that SM fields live on a 3-brane ($`d=4`$), while gravitons inhabit the entire spacetime ($`d=10`$). The prevalent energies of the decay quanta are $`T_{\mathrm{BH}}1/r_s`$, resulting in $`s`$-wave dominance of the final state. Indeed, as the total angular momentum number of the emitted field increases, $`\sigma _s(Q,r_s)`$ rapidly gets suppressed Kanti:2002nr . In the low energy limit, $`Qr_s1,`$ higher order terms get suppressed by a factor of $`3(Qr_s)^2`$ for fermions and $`25(Qr_s)^2`$ for gauge bosons. For an average particle energy, $`Qr_s^1,`$ higher partial waves also get suppressed, though the suppression is not as large. This strongly suggests that the BH is only sensitive to the radial coordinate and does not make use of the extra angular modes available in the internal space Emparan:2000rs . A recent numerical study Harris:2003eg  has explicitly shown that the emission of scalar modes into the bulk is largely suppressed with respect to the brane emission. In order to contravene the argument of Emparan–Horowitz–Myers Emparan:2000rs , the bulk emission of gravitons would need to show the opposite behavior – a substantial enhancement into bulk modes. There is no a priori reason to suspect this qualitative difference between $`s=0`$ and $`s=2`$, and hence no reason to support arguments Cavaglia:2003hg  favoring deviation from the dominance of visible decay. With this in mind, we assume the evaporation process to be dominated by the large number of SM brane modes.

The total number of particles emitted is roughly equal to the BH entropy,

$$S_{\mathrm{BH}}=\frac{\pi }{2}M_{\text{BH}}r_s.$$
(4)
At a given time, the rate of decrease in the BH mass is just the total power radiated

$$\frac{d\dot{M}_{\mathrm{BH}}}{dQ}=\underset{i}{}c_i\frac{\sigma _s(Q,r_s)}{8\pi ^2}Q^3\left[\mathrm{exp}\left(\frac{Q}{T_{\mathrm{BH}}}\right)(1)^{2s}\right]^1,$$
(5)
where $`c_i`$ is the number of internal degrees of freedom of particle species $`i`$. Integration of Eq. (5) leads to

$$\dot{M}_{\mathrm{BH}}=\underset{i}{}c_if\frac{\mathrm{\Gamma }_s}{8\pi ^2}\mathrm{\Gamma }(4)\zeta (4)T_{\mathrm{BH}}^4A_4,$$
(6)
where $`f=1`$ ($`f=7/8`$) for bosons (fermions), and the greybody factor was conveniently written as a dimensionless constant, $`\mathrm{\Gamma }_s=\sigma _s(Q,r_s)/A_4`$, normalized to the BH surface area Emparan:2000rs

$$A_4=\frac{36}{7}\pi \left(\frac{9}{2}\right)^{2/7}r_s^2$$
(7)
seen by the SM fields ($`\mathrm{\Gamma }_{s=1/2}0.33`$ and $`\mathrm{\Gamma }_{s=1}0.34`$ greybody ). Now, since the ratio of degrees of freedom for gauge bosons, quarks and leptons is 29:72:18 (the Higgs boson is not included), from Eq. (6) one obtains a rough estimate of the mean lifetime,

$$\tau _{_{\mathrm{BH}}}1.67\times 10^{27}\left(\frac{M_{\text{BH}}}{M_{10}}\right)^{9/7}\left(\frac{\mathrm{TeV}}{M_{10}}\right)\mathrm{s},$$
(8)
which indicates that BHs evaporate instantaneously into visible quanta. This is within a factor of 2 of the heuristic estimate made in Ref. Argyres:1998qn ,

$$\tau _{_{\mathrm{BH}}}6.58\times 10^{28}\left(\frac{M_{\text{BH}}}{M_{10}}\right)^{9/7}\left(\frac{\mathrm{TeV}}{M_{10}}\right)\mathrm{s},$$
(9)
which is generally used in the literature.

The semiclassical description outlined above is only reliable when the energy of the emitted particle is small compared to the BH mass, i.e.,

$$T_{\mathrm{BH}}M_{\mathrm{BH}},\mathrm{or}\mathrm{equivalently},M_{\text{BH}}M_{10},$$
(10)
because it is only under this condition that both the gravitational field of the brane and the back reaction of the metric during the emission process can be safely neglected Preskill:1991tb . For BH with initial masses well above $`M_{10}`$, most of the decay process can be well described within the semiclassical approximation. However, the condition stated in Eq. (10) inevitably breaks down during the last stages of evaporation. At this point it becomes necessary to introduce quantum considerations. To this end we turn to a quantum statistical description of highly excited strings.

It is well-known that the density of string states with mass between $`M`$ and $`M+dM`$ cannot increase any faster than $`\rho (M)=e^{\beta _HM}/M,`$ because the partition function,

$$Z(\beta )=_0^{\mathrm{}}𝑑M\rho (M)e^{M\beta },$$
(11)
would fail to converge Hagedorn:st . Indeed, the partition function converges only if the temperature is less than the Hagedorn temperature, $`\beta _H^1`$, which is expected to be $`M_s`$. As $`\beta `$ decreases to the transition point $`\beta _H`$, the heat capacity rises to infinity because the energy goes into the many new available modes rather than into raising the kinetic energy of the existing particles Frautschi:1971ij . In the limit, the total probability diverges, indicating that the canonical ensemble is inadequate for the treatment of the systems. However, one can still employ a microcanonical ensemble of a large number of similar insulated system each with a given fixed energy $`E`$. With the center-of-mass at rest, $`E=M`$. This means that the density of states is just $`\rho (M)`$ and the entropy $`S=\mathrm{ln}\rho (M)`$. In this picture the equilibrium among systems is governed by the equality of the temperatures, defined for each system as

$$T\left(\frac{S}{M}\right)^1=\frac{M}{\beta _HM1}.$$
(12)
Equilibrium is achieved at maximum entropy when the total system heat capacity, $`C`$, is positive. Ordinary systems (on which our intuition is founded) have $`C>0`$. However, for a gas of massive superstring excitations the heat capacity,

$$C\frac{1}{T^2}\left(\frac{^2S}{M^2}\right)^1=\left(\frac{M}{T}\right)^2,$$
(13)
is negative, as is the case for BHs Hawking:de . The positivity requirement on the total specific heat implies that strings and BHs cannot coexist in thermal equilibrium, because any subsystem of this system has negative specific heat, and thus the system as a whole is thermodynamically unstable. This observation suggests that BHs may end their Hawking evaporation process by making a transition to an excited string state with higher entropy, avoiding the singular zero-mass limit Bowick:1985af . The suggestion of a string $``$ BH transition is further strengthened by three other facts: (i) in string theory, the fundamental string length should set a minimal radius for the Schwarzschild radius of any BH Veneziano:1986zf ; (ii) $`T_{\mathrm{BH}}\beta _H^1`$ for $`r_sM_s^1`$ Susskind:ws . (iii) There is a seeming correlation between the greybody factors in BH decay and the level structure of excited strings Das:1996wn . The string $``$ BH “correspondence principle” Horowitz:1996nw  unifies these concepts: When the size of the BH horizon drops below the size of the fundamental string length $`\mathrm{}_s\mathrm{}_{10},`$ where $`\mathrm{}_{10}`$ is the fundamental Planck length, an adiabatic transition occurs to an excited string state. Subsequently, the string will slowly lose mass by radiating massless particles with a nearly thermal spectrum at the unchanging Hagedorn temperature Amati:1999fv .(The probability that a BH will radiate a large string, or else that a large string would undergo fluctuation to become a BH is very small Horowitz:1997jc .)

The continuity of the cross section at the correspondence point, at least parametrically in energy and string coupling, provides independent supportive argument for this picture Dimopoulos:2001qe . Specifically, (in the perturbative regime) the Virasoro-Shapiro amplitude leads to a “string ball” (SB) production cross section $`g_s^2\widehat{s}/M_s^4`$. This cross section saturates the unitarity bounds at $`g_s^2\widehat{s}/M_s^21`$ Amati:1987wq , so before matching the geometric BH cross section $`r_s^2`$ there is a transition region at which $`\widehat{\sigma }M_s^2`$. All in all, the rise with energy of the parton-parton $``$ SB/BH cross section can be parametrized as Dimopoulos:2001qe

$`\widehat{\sigma }(\sqrt{\widehat{s}})\{\begin{array}{cc}{\displaystyle \frac{g_s^2\widehat{s}}{M_s^4}}\hfill & M_s\sqrt{\widehat{s}}M_s/g_s,\hfill \\ {\displaystyle \frac{1}{M_s^2}}\hfill & M_s/g_s<\sqrt{\widehat{s}}M_s/g_s^2,\hfill \\ {\displaystyle \frac{1}{M_{10}^2}}\left[{\displaystyle \frac{\sqrt{\widehat{s}}}{M_{10}}}\right]^{2/7}\hfill & M_s/g_s^2<\sqrt{\widehat{s}},\hfill \end{array}`$ (17)
where $`M_{10}=(8\pi ^5)^{1/8}M_s/g_s^{1/4}.`$

Colliders have not yet attained the energies required to probe TeV scale gravitational collapse; nevertheless cosmic rays may already provide some clues. To date a handful of cosmic rays have been observed with energy in excess of $`10^{11}`$ GeV Anchordoqui:2002hs . When these particles impinge on a stationary nucleon in the upper atmosphere they probe center-of-mass energies as high as $`\sqrt{s}400`$ TeV. In what follows we estimate the sensitivity of cosmic ray observatories to SB/BH production.

## III Probes of large extra dimensions with cosmic rays

The SB/BH production cross section, $`𝒪(M_{\mathrm{EW}}^1)`$, is about 5 orders of magnitude smaller than QCD cross sections, $`𝒪(\mathrm{\Lambda }_{\mathrm{QCD}}^1),`$ thus making it futile to hunt for SB/BHs in hadronic cosmic rays. On the other hand, SM neutrinos interact so weakly that any quantum gravitational enhancement of the cross section can be experimentally distinguished from background. A substantial neutrino component may accompany the observed cosmic ray flux because of neutrino production in cosmic beam dumps (such as gamma ray burst fireballs Razzaque:2002kb , $`X`$-ray binaries Anchordoqui:2002xu , or active galactic nuclei Szabo:qx ). In addition, the inelastic collision of ultra-high energy ($`>10^{10.6}`$ GeV) nucleons with the relic photons permeating the universe should deplete the hadronic cosmic ray intensity at the end of the spectrum and breed neutrinos from decay products of charged pions Stecker:1979ah . The reaction generating these cosmogenic neutrinos is well known physics and the flux depends only on the existence of ultra-high energy nucleon sources more distant than $`8`$ Mpc. In this case, the expected $`\nu `$-flux can be accurately predicted by fitting the observed spectrum with a homogeneous population of nucleon sources Protheroe:1996ft . The existence of the cosmogenic neutrinos is a safe bet, implying a minimum guaranteed flux for testing any new physics that promotes the neutrino cross section to a sub-hadronic size.111In the simple $`T^6`$ compactification, there are some unitarization procedures for the multi Kaluza-Klein graviton exchange amplitudes which can yield hadronic cross sections at $`E_\nu 10^{10}`$ GeV Jain:2000pu . However, these cross sections have quadratic energy dependence, and so will generate moderately penetrating showers with definite profiles at lower energies Anchordoqui:2000uh . These have not been reported to date, and their absence places serious constraints on the enhancement of the neutrino-nucleon cross section to sub-hadronic size.

The inclusive production of BHs proceeds through different final states for different classical impact parameters $`b`$ Yoshino:2002br . These final states are characterized by the fraction $`y(z)`$ of the initial parton center-of-mass energy, $`\sqrt{\widehat{s}}=\sqrt{xs}`$, which is trapped within the horizon. Here, $`z=b/b_{\mathrm{max}},`$ and $`b_{\mathrm{max}}=1.3r_s(\sqrt{\widehat{s}})`$ Yoshino:2002br . With a lower cutoff $`M_{\mathrm{BH},\mathrm{min}}`$ on the BH mass required for the validity of the semi-classical description, this implies a joint constraint

$$y(z)\sqrt{xs}M_{\mathrm{BH},\mathrm{min}}$$
(18)
on the parameters $`x`$ and $`z`$. Because of the monotonically decreasing nature of $`y(z)`$, Eq. (18) sets an upper bound $`\overline{z}(x)`$ on the impact parameter for fixed $`x.`$ The corresponding parton-parton BH cross section is $`\widehat{\sigma }_{_{\mathrm{BH}}}(x)=\pi \overline{b}^2(x),`$ where $`\overline{b}=\overline{z}b_{\mathrm{max}}.`$ The total BH production cross section is then Anchordoqui:2003jr

$$\sigma _{_{\mathrm{BH}}}(E_\nu ,M_{\mathrm{BH},\mathrm{min}},M_{10})_{\frac{M_{\mathrm{BH},\mathrm{min}}^2}{y^2(0)s}}^1𝑑x\underset{i}{}f_i(x,Q)\widehat{\sigma }_{_{\mathrm{BH}}}(x),$$
(19)
where $`i`$ labels parton species and the $`f_i(x,Q)`$ are parton distribution functions (pdfs) Pumplin:2002vw . The momentum scale $`Q`$ is taken as $`r_s^1,`$ which is a typical momentum transfer during the gravitational collapse. A useful criterion for a BH description is that the entropy is sizable, so that Eq. (10) is satisfied. For $`M_{\mathrm{BH},\mathrm{min}}=3M_{10}`$, $`S_{\mathrm{BH}}131`$.

In the perturbative string regime, i.e., $`M_{\mathrm{SB},\mathrm{min}}<\sqrt{\widehat{s}}M_s/g_s`$, the SB production cross section is taken as

$$\sigma _{_{\mathrm{SB}}}(E_\nu ,M_{\mathrm{SB},\mathrm{min}},M_{10})=_{\frac{M_{\mathrm{SB},\mathrm{min}}^2}{s}}^1𝑑x\underset{i}{}f_i(x,Q)\widehat{\sigma }_{_{\mathrm{SB}}}(\widehat{s}),$$
(20)
where $`\widehat{\sigma }_{_{\mathrm{SB}}}(\widehat{s})`$ contains the Chan Paton factors which control the projection of the initial state onto the string spectrum. In general, this projection is not uniquely determined by the low-lying particle spectrum, yielding one or more arbitrary constants. The analysis in the $`\nu q\nu q`$ channel illustrates this point Cornet:2001gy . The $`\nu g`$ scattering, relevant for $`\nu N`$ interactions at ultra-high energies, introduces additional ambiguities. For simplicity, in the spirit of Dimopoulos:2001qe , we take $`\widehat{\sigma }_{_{\mathrm{SB}}}=(\pi ^2\text{g}_s^2\widehat{s})/(8M_s^4),`$ derived from the Virasoro-Shapiro amplitude. The momentum scale in this case varies from $`Q\sqrt{\widehat{s}}`$ at parton center-of-mass energies corresponding to low-lying string excitations, to the BH value $`Qr_s^1`$ for $`\sqrt{\widehat{s}}M_s/g_s^2,`$ the transition point energy. This corresponds to very high string excitations energies. In practice, this may be approximated as a two stage designation:

$`Q(\sqrt{\widehat{s}})\{\begin{array}{cc}\sqrt{\widehat{s}}\hfill & \mathrm{if}\sqrt{\widehat{s}}M_s/g_s^2,\hfill \\ & \\ r_s^1\hfill & \mathrm{if}M_s/g_s^2<\sqrt{\widehat{s}},\hfill \end{array}`$ (24)
with small ($`10\%`$) uncertainty Anchordoqui:2002vb .

In Fig. 1 we show a comparison between the total BH production cross section and the particular model we are considering for SBs. As can be seen in Fig. 1, for $`E_\nu 10^7`$ GeV, $`\sigma _{_{\mathrm{SB}}}>\sigma _{_{\mathrm{BH}}}>\sigma _{_{\mathrm{SM}}}`$. This makes plausible that in the quantum range, $`M_{10}<M_{\mathrm{BH},\mathrm{min}}<3M_{10},`$ the BH cross can be thought of as a lower bound on the anomalous $`\nu N`$ cross section. Since the BH production cross section is suppressed by $`M_{10}^2,`$ bounds on the $`\nu N`$ cross section can be translated into lower limits on $`M_{10}.`$ If the inequality above holds true in more realistic string models, the bounds on $`M_{10}`$ can be significantly strengthened by taking $`M_{\mathrm{BH},\mathrm{min}}=M_{10}.`$

Several techniques have been used to search for the cosmic neutrino component Halzen:2002pg . Showers in the atmosphere can be detected either by observing the fluorescence or Čerenkov light induced by the ionizing tracks in the air, or by directly detecting the charged particles in the shower tail with scintillation counters scattered on the Earth’s surface. Neutrino showers mediated by SB/BHs may be distinguished from hadronic primary particles by looking for very inclined showers. This is because for a horizontal path (90 zenith angle) the neutrino interaction length for SB/BH production is larger than the thickness of the atmosphere (36000 g/cm2) and thus, unlike hadrons that shower high in the atmosphere, neutrinos could travel through most of this matter before gravitational collapse is triggered. An alternative technique exploits naturally occurring large volume Čerenkov radiators such as the Antarctic ice-cap. Several radio antennae monitor cold ice for radio frequency Čerenkov radiation resulting from neutrino in-ice cascades. Potential signal events are distinguished from background using vertex location and signal shape.

The experimental situation is that, in spite of extensive searching, only one neutrino-like event has been reported, with an expected background from hadronic cosmic rays of 1.72 events. This implies, at 95% C.L., a maximum of 3.5 events, which combined with the sum of the exposures of the various experiments, and convoluted with the guaranteed cosmogenic neutrino flux, allows the extraction of a lower bound on the fundamental Planck scale, $`M_{10}>1.01.4`$ TeV, for $`M_{10}<M_{\mathrm{BH},\mathrm{min}}<3M_{10}`$ Anchordoqui:2003jr . This bound is competitive with that obtained from Tevatron data Abbott:2000zb , and these represent the most stringent lower limits to date for the size of 6 large extra dimensions.

## IV Concluding remarks

We have reviewed the possibility of searching for black holes and superstring excitations in cosmic rays using neutrino interactions in the Earth’s atmosphere and the Antarctic ice. The fact that neutrino showers have not yet been observed allows us to put a lower bound on the fundamental Planck scale $`M_{10}>1.01.4`$ TeV. As discussed previously, the string $``$ BH correspondence principle provides support for the higher bound. The next generation of cosmic ray experiments Anchordoqui:2001cg ; Ringwald:2001vk  and dedicated neutrino telescopes Kowalski:2002gb  should probe values of $`M_{10}`$ up to about 5 TeV.

Conversely, if neutrino showers are observed above predicted SM rate, they can be ascribed to either new physics or a larger than expected cosmogenic flux. The most powerful technique to distinguish between these two possibilities exploits the absorption of SB/BH secondaries by the Earth, by separately binning Earth skimming events Feng:2001ue  which arrive at very small angles to the horizontal. An enhanced flux will increase both quasi-horizontal and Earth-skimming event rates, whereas a large SB/BH cross section suppresses the latter, because the hadronic decay products of SB/BH evaporation do not escape the Earth’s crust Anchordoqui:2001cg . If such a suppression is observed then string balls and black holes will provide a compelling explanation.

###### Acknowledgements.

LAA would like to thank Sacha Davidson and the organizers of COSMO-03 for a most stimulating meeting. The work of LAA and HG has been supported in part by the US National Science Foundation (NSF), under grants No. PHY–0140407 and No. PHY–0073034, respectively. The work of JLF is supported in part by NSF CAREER grant No. PHY-0239817. The work of ADS is supported in part by Department of Energy Grant No. DE–FG01–00ER45832 and NSF Grants PHY–0071312 and PHY–0245214.


text:
# Infrared Features of the Lattice Landau Gauge QCD

## 1 The gluon propagator and the ghost propagator

We measured $`SU(3)`$ gluon propagator of $`\beta =6.0`$ on lattice $`16^4`$, $`24^4`$, $`32^4`$ and $`\beta =6.4`$ on $`32^4`$, $`48^4`$. Shown in Fig.1 is the gluon dressing function $`Z(q^2)`$ which is defined by the gluon propagator $`D_{\mu \nu }(q)`$ through

$$D_{\mu \nu }(q)=(\delta _{\mu \nu }\frac{q_\mu q_\nu }{q^2})D_A(q^2)$$
(1)
as $`Z(q^2)=q^2D_A(q^2)`$. Here, we adopted the $`\mathrm{log}U`$ type definition of gauge field $`A`$. Although the size of the sample is not large, the data scale with each other and consistent with the result of. The fitted line is obtained from the model based on the principle of minimal sensitivity(PMS). The PMS allows one to obtain a solution that bridges infrared region and ultraviolet region, but its uniqueness is not guaranteed and in the infrared region, there is discrepancy from other scheme that satisfies the commensurate scale relation. We leave the continuation to the infrared region to the future study and fix the parameter $`y`$ by the solution at the factorization scale i.e. $`q=1.97GeV`$ and call it the $`\stackrel{~}{MOM}`$ scheme.

The ghost propagator is defined as the Fourier transform of an expectation value of the inverse Faddeev-Popov operator $``$

$$D_G^{ab}(x,y)=\mathrm{tr}\lambda ^ax|([U])^1|\lambda ^by$$
(2)
The lattice data scales and the PMS in $`\stackrel{~}{MOM}`$ scheme fits in the region $`q>0.5GeV`$.

## 2 QCD running coupling and Kugo-Ojima parameter

In the contour-improved perturbation theory, physical quantities $``$ are expressed in a series

$$(q^2)=_1(q^2)+\underset{n=1}{\overset{\mathrm{}}{}}A_n_{n+1}(q^2)$$
(3)
where $`A_n`$ is known by perturbation theory in loop expansion, and $`_n`$ are expressed by the Lambert W function. The result of running coupling using $`\mathrm{\Lambda }_{\overline{MS}}=237MeV`$ fits the data of $`\beta =6.4,48^4`$ data, except the peak around $`1GeV`$.

In the $`\beta =6.4,48^4`$ lattice in $`\mathrm{log}U`$ version, we obtained the Kugo-Ojima parameter $`c=0.793(61)`$, trace $`e/d=0.982(1)`$ and the horizon function deviation $`h=0.19`$.

## 3 Gauge fixing in use of PT and Gribov copy problem in $`SU(2)`$

Fundamental modular gauge (FMG: global maximum of the optimizing function $`F_U(g)`$) is the unique gauge known as the legitimate Landau gauge without Gribov copies, where $`F_U(g)=_{\mathrm{}}\mathrm{Re}\mathrm{tr}U_{\mathrm{}}^g`$. Simulated annealing methods are now widely used to cope with Gribov copies in the maximal abelian gauge. In the present work, we combine parallel tempering (PT) method as a new trial with gauge fixing for FMG, and test it in $`SU(2)`$ gauge theory on $`16^4`$ lattice with $`\beta =2.2`$ of the Wilson action. The PT was invented as one of suitable methods for simulating correctly low temperature equilibrium of systems possessing complex energy profile like spin glass. We do not need the exact equilibrium of Boltzmann distribution $`\mathrm{Prob}(g)exp\{\beta F_U(g)\}`$, but only the lowest energy state $`g`$ of $`Min_gE(g)`$ efficiently where $`E(g)=F_U(g)`$. In the PT method, we prepare $`N_r`$ replicas of site-spin (gauge transformation) system with gauge link interaction $`U`$. Each replica $`g_j`$ is given different $`\beta _{i_j}`$, $`\beta _{min}\beta _i\beta _{max}`$. After $`N_{mc}`$ Monte Carlo steps at the given $`\beta `$, the $`\beta `$’s are exchanged between adjacent values according to Metropolis acceptance check. During $`N_{ex}`$ iterations of the $`\beta `$-exchange step, the lowest energy state among all replicas is monitored and stored. These accumulated $`N_c`$ gauge spins are all fixed to the Landau gauge by the standard over-relaxation (OR) method ($`\mathrm{Max}_x(A)_x<10^4`$). We pick up the best copy among those, and name it as PT copy. In the PT method, the distribution of $`\beta `$’s in the interval $`[\beta _{min},\beta _{\mathrm{max}}]`$ is technically the most important point so that the exchange rates are almost equal everywhere between replicas. We found the good choice for adjustment of the new $`\beta ^{}`$’s from the old $`\beta `$’s as $`\mathrm{}\beta ^{}\sqrt{\overline{\mathrm{\Delta }}/(\mathrm{}E/\mathrm{}\beta )}`$ where $`\mathrm{\Delta }=\mathrm{}\beta \mathrm{}E`$ and $`\overline{\mathrm{\Delta }}`$ is the replica average of $`\mathrm{\Delta }`$. Due to the actual situation of our choice, $`\beta _{min}=0.4`$ and $`\beta _{max}=30`$, and the system size of $`16^4`$, and $`N_r=24`$, direct application of the PT is practically meaningless since the Metropolis acceptance rate is very small. We dare to shrink the energy scale by a factor $`\alpha =0.0002`$ only within the $`\beta `$ exchange procedure to have the Metropolice acceptance rate around $`0.5`$. Other parameters of ours are as follows, $`N_{mc}=120`$, $`N_{ex}=12`$, $`N_c=1632`$. Total number of PT gauge fixed samples of $`\beta =2.2`$ of $`SU(2)`$ Wilson action is 67, and we compared them with gauge copies that are directly gauge fixed by OR method and are named as first copies. Among 67 samples, 63 PT copies give larger $`F_U(1)`$ with $`N_c=16`$ than first copies and the remaining 4 PT copies with $`N_c=2032`$. Using these two groups of copies, Gribov noise is investigated in Kugo-Ojima parameter, gluon propagator and ghost propagator. The averaged Kugo-Ojima parameter is suppressed about 4% in PT in comparison with first copy. Effects of noise about 6% are seen in Fig.4 for the ghost propagator, but no significant effects in the gluon propagator. For a few samples among 67, comparison is made for each sample between two groups of $`16(=N_c)`$ PT generated copies, and 32 copies generated(CRGT) by direct OR gauge fixing after random gauge transformation, in a scatter plot of Kugo-Ojima parameter vs $`F_U(1)`$. It is observed in this plot that a group of PT copies scatter almost in a line on top-end of $`F_U(1)`$, with all CRGT scattering in lower $`F_U(1)`$ area, and this fact raises questions if this ”flat-end” feature is real one of the FMG, or if we have not reached the real FMG yet? This is a rather unexpected feature with respect to FMG.

## 4 Discussion and outlook

The momentum dependences of the gluon propagator and the ghost propagator in the infrared are analyzed in use of PMS method. Deviation of $`u(0)=c`$ from -1 (Kugo-Ojima confinement criterion) does not necessarily imply violation of colour symmetry. It could imply a question on the realization of BRST quartet mechanism on the lattice and/or even further validity of the continuum Lagrangian in the infrared region. So far as our data in off-diagonal elements of ghost propagator show, there appear no definite signals of colour global symmetry breaking.

From the scatter plot the Kugo-Ojima vs $`F_U(1)`$, it seems that our PT gauge fixing works well, but more careful parameter tuning and comparison with simulated annealing method deserve a future study. Gribov noise in the ghost propagator is qualitatively similar to Cucchieri’s noise at smaller $`\beta =1.6`$. The infrared singularity is weakly suppressed in FMG.

Some revisions in numerical data of ghost and Kugo-Ojima parameter have been made in this paper due to use of conjugate gradient method.

This work was supported by the KEK supercomputing project No.03-94.


text:
# Structure and Evolution of Starburst and Normal Galaxies

## 1 INTRODUCTION

A coherent picture of the formation and evolution of galaxies can only emerge through an understanding of starbursts over a wide range of cosmic lookback time. Starbursts drive galactic superwinds (e.g., Lehnert & Heckman 1996) which play a key role in the metal enrichment of the intergalactic medium and the establishment of the strong mass-metallicity relation in ellipticals and bulges (Lynden-Bell 1992). Moreover, the evolution of starbursts and active galactic nuclei (AGN) may be intimately linked through a common central gas supply and evolutionary effects (e.g., Norman & Scoville 1988; Kauffmann & Haehnelt 2000).

To date, detailed studies of star-formation have primarily targeted galaxies at high redshifts ($`z2.55`$) using optical (e.g., Steidel et al. 1999) and submillimeter (e.g., Hughes et al. 1998; Barger et al. 2001) emission. Conversely, in the local Universe ($`z<1`$), studies of starbursts have mainly focused on dusty ultra-luminous infrared galaxies (ULIRGs) which show signs of strong interaction or double nuclei (Joseph & Wright 1985; Armus, Heckman & Miley 1987; Scoville et al 2000) and large gas densities (Sanders, Scoville, Soifer 1991). These are likely counterparts of high redshift sub-mm sources (Chapman et al 2003).

This paper focuses on a population of starburst galaxies which is different from those studied before. The sample of starbursts studied here is optically selected, complete to $`R_{AB}=24`$ mag, and extends to $`z1`$. Therefore, the sample preferentially consists of relatively dust-free blue starbursts which differ in nature or in evolutionary stage from those selected at infrared wavelengths. Using multi-waveband high resolution images from the Advanced Camera for Surveys (ACS), obtained as a part of the Great Observatories Origins Deep Survey (GOODS), we compare the optically-selected starbursts at $`z1`$ to a control sample consisting of normal late and early type systems, and discuss their role in the evolution of these galaxies.

## 2 Observations and Sample Selection

### 2.1 Observations

The GOODS project consists primarily of space-borne ($`HST`$/ACS) and ground-based multi-waveband imaging of two fields, the HDF-N and CDF-S. In this study we consider only the Southern GOODS field (CDF-S) where most of the initial data have been obtained and processed.

The GOODS-ACS observations are performed in four bands ($`BViz`$), covering an area of $`10\mathrm{}\times 16\mathrm{}`$. The ground-based imaging observations cover the entire CDF-S area in $`U`$$`UBVRIJHK`$ to the depths of $`R_{\mathrm{AB}}`$=25 mag and $`K_{\mathrm{s}_{\mathrm{AB}}}`$=23 mag. A more detailed discussion of the GOODS observations is given in Giavalisco et al (2003).

### 2.2 Photometric Redshifts and Spectral types

The ground-based CDF-S data provide photometry in 14 passbands consisting of WFI ($`U`$$`UBVRI`$), FORS ($`RI`$), SOFI, and ISAAC ($`JHK`$). The observations cover the entire field except for the FORS and ISAAC observations which are deeper and only cover a sub-area of the field. Magnitudes measured over an aperture with a $`3\mathrm{}`$ diameter were used to construct the observed SEDs for all galaxies in the $`R`$-band selected survey. Photometric redshifts and spectral types were estimated for individual galaxies by fitting the observed and template SEDs using the Bayesian photometric redshift technique (Benitez 2000). The template SEDs consist of E, Sbc, Scd, and Im (Coleman, Wu & Weedman 1980), and starburst (Kinney et al 1996). A sub-sample of galaxies in the CDF-S have spectroscopic redshifts from the K20 survey (Cimatti et al 2002). Comparison between the estimated photometric and spectroscopic redshifts for these galaxies gives $`\frac{z_{phot}z_{spec}}{(1+z_{spec})}0.10`$, which is taken as the error in photometric redshifts here. Additionally, the derived spectral types of the galaxies were found to be consistent with visually classified elliptical and spiral galaxies on the ACS images. A detailed study of the photometric redshift measurements for the GOODS field will be presented in Mobasher et al (2003).

### 2.3 Sample Selection

The starburst sample in this paper includes galaxies from the $`R`$-band selected ground-based survey which satisfy the following criteria: (a) They have starburst spectral types and $`R_{\mathrm{AB}}<24`$ mag. The magnitude limit ensures that photometric redshifts and spectral types are measured more accurately ($`\frac{\mathrm{\Delta }(z)}{(1+z_{spec})}<0.1`$). (b) They have redshifts in the range $`0.2<z<1.3`$. Lower redshifts are excluded to minimize relative uncertainties in photometric redshifts and resolution dependent effects (see $`\mathrm{\S }`$ 3). (c) They lie in the CDF-S area covered by the ACS. A control sample was also constructed (to the same magnitude limit) consisting of normal galaxies (i.e. non-starbursts), with early (E/Sa) and late (Sb–Sd) spectral types. Im galaxies are not included.

Given the selection criteria, our sample of starburst galaxies is biased towards relatively dust-free blue starbursts and against dusty star-forming galaxies. This study thus complements previous investigations of dusty starbursts and ULIRGs which were selected based on their infrared properties (e.g., Scoville et al 2000). It is to be noted that the SED fitting method used here does not fully guarantee that the selected objects are indeed starbursts. For example, a very metal poor galaxy with only modest amount of star formation could have a very blue color and hence, be classified as a starburst here. However, the number of such galaxies is expected to be small with the sample mostly dominated by galaxies experiencing bursts of star formation.

Table 1 presents the number of starburst and control galaxies in their observed frames, so that in the redshift intervals indicated, they correspond to the rest-frame $`B`$-band. The median absolute magnitudes for the three different types of galaxies in the present sample ($`R_{\mathrm{AB}}<24`$) are $`M_B=19.6`$ (starburst), $`19.9`$ (Sb–Sd) and $`20.6`$ (E/Sa). The median $`UV`$ ($`VK`$) colors for the starburst, late and early-type galaxies here are respectively estimated as $`0.70`$ (1.80), $`0.20`$ (2.40) and 0.35 (2.8). In order to check the sensitivity of our results to the adopted magnitude limit, the analysis in this letter ($`\mathrm{\S }`$ 3) was repeated to $`R_{\mathrm{AB}}<`$ 25 mag, with no difference found.

## 3  Structure and Morphology

We quantify the concentration and asymmetry of the sample galaxies using the CAS concentration ($`C`$) (Bershady et al. 2000) and asymmetry ($`A`$) indices (Conselice et al. 2000). The concentration index is proportional to the logarithm of the ratio of the 80 % to 20 % curve of growth radii

$$C=5\times log(r_{80\%}/r_{20\%}).$$
(1)
The asymmetry index is obtained by rotating a galaxy image by $`\pi `$, subtracting it from its pre-rotated image, summing the intensities of the absolute value residuals, and normalizing the sum to the original galaxy flux. The concentration and asymmetry parameters are measured for the starburst sample in the rest-frame $`B`$-band ($`C_\mathrm{B}`$, $`A_\mathrm{B}`$) to $`z1`$, using different observed ACS bands in the appropriate redshift ranges, as shown in Table 1. Corresponding quantities were also derived in the rest-frame $`R`$-band ($`C_\mathrm{R}`$, $`A_\mathrm{R}`$). The rest-frame $`B`$ and $`R`$ bands are respectively sensitive to young stars and the underlying older stellar population.

In order to minimize resolution dependent effects and relative uncertainties in photometric redshifts, we constrain the redshift range to $`z`$ 0.2–1.3, where the change in spatial resolution ($`0.05\mathrm{}`$) is only a factor $``$ 2, ranging from 160 pc to 400 pc, assuming a flat cosmology with $`\mathrm{\Omega }_M=1\mathrm{\Omega }_\mathrm{\Lambda }=0.3`$ and a Hubble constant $`H_0`$ = 70 km s-1 Mpc-1. The structural comparisons in $`\mathrm{\S }`$ 3.1–3.2 probe galaxy evolution over a time period where the age of the Universe varied from 6 to 11 Gyr. It is important to note that the starburst and control samples are defined purely based on SED fits to multi-band ground-based images, without using any morphological information. Thus, the derived structural parameters provide additional independent probes of the nature and evolutionary stages of these systems.

### 3.1 Concentration and AGN fraction

Figure 1 shows the distribution of rest-frame $`B`$ concentration indices ($`C_\mathrm{B}`$) for the starburst galaxies and the control sample, divided into two redshift bins. The distribution of $`C_\mathrm{B}`$ in early-type systems is significantly skewed towards higher values compared to starburst galaxies. Results from both the $`T`$-test and the Kolmogorov-Smirnov test, listed in Table 2, confirm that the distribution of concentration indices for starburst and early-type galaxies are statistically different, both in the mean and the overall distribution. About 72 % of the early-type galaxies have $`C_\mathrm{B}`$ $`>`$ 3.0 compared to 12 % of the starburst and 18 % of the late-type galaxies (Table 2).

The data suggest that by $`z`$ 1, early-type systems have already developed large central light concentrations, which demarcate them from optically selected starbursts and late-type galaxies (Figs. 1a and 1b). These large concentrations are likely associated with bulges and spheroidal components. This was demonstrated by simulation of galaxies with de Vaucouleurs (r1/4) law, which were found to have large $`C_\mathrm{B}`$ ($`>`$ 3.4) values (Jogee et al. 2003, in prep.). Furthermore, comparison of our samples with the CDF-S X-ray catalog (Alexander et al 2003) shows that a much larger fraction ($`>`$ 25 %) of the early-type systems host AGNs compared to only 2% of the starbursts111Assuming $`1.5arcsec`$ as the maximum distance between the ACS $`B`$-band coordinates and the X-ray identified sources (Koekemoer et al. 2003). Extending the search to $`6.0\mathrm{}`$ increases the number of starbursts to 7%.. The larger AGN fraction and high $`C_\mathrm{B}`$ we find in early-types are in qualitative agreement with the high concentrations in AGN hosts reported by other studies (e.g., Grogin at al. 2003). One general scenario consistent with our results is that the early-type systems with large $`C_\mathrm{B}`$ developed their central structures (bulges and spheroids) in much earlier ($`z`$1) violent gas-rich dissipative episodes such as major mergers (e.g., Naab & Burkert 2001). A concurrent growth of bulges/spheroids and massive black holes may occur in some of these mergers (Kauffmann & Haehnelt 2000), consistent with the large AGN fraction in the early types.

What is the nature of the optically selected starbursts at $`z1`$ and what will they evolve into? Most of these galaxies have concentrations in the range $`C_\mathrm{B}`$=3.3–4.5, which are well below those shown by the majority ($`>`$ 50 %) of early-type galaxies, but comparable to those of late-type spirals in the same redshift range. Low values in $`C_\mathrm{B}`$ can result from disk-like systems as revealed by fits to artificial galaxies with exponential disks (Jogee et al. 2003, in prep.). However, $`C_\mathrm{B}`$ in starbursts is more strongly influenced by very young stars than in early-type systems so that low $`C_\mathrm{B}`$ may also result from localized, off-centered star-forming regions which often exist in early stages of tidal interactions or mergers. Visual inspection of the starbursts confirms that they include disk-like systems as well as a large number of morphologically disturbed objects which show tails, double nuclei, and nearby companions with a range of luminosity ratios (1:10 to 1:1). In summary, the lower $`C_\mathrm{B}`$ and AGN fraction in the majority of optically-selected starbursts, compared to early-type systems, suggest two possibilities. Either the starbursts and early types are entirely different in nature (being respectively disk and bulge dominated), or/and are in different evolutionary phases (that some of the starbursts in major mergers evolve into early-types as violent relaxation sets in and starburst spectral signature fades out).

### 3.2 Asymmetry

Figure 2a shows the distributions of rest-frame asymmetry ($`A_\mathrm{B}`$) indices in the rest-frame $`B`$-band for the starburst and control sample. A larger fraction (50 %) of the starburst galaxies have a high $`A_\mathrm{B}`$ $`>`$0.3, compared to early-types (13 %) and late-types (27 %) (Table 2). It is likely that large values of $`A_\mathrm{B}`$ are driven by irregular patterns of star formation which may either be externally triggered in tidally interacting/merging systems or spontaneously triggered in isolated galaxies. Asymmetries induced by dust lanes, caused by high inclinations, are not significant as most objects in our sample are moderately inclined, based on their estimated ellipticities. Fig. 2b shows the corresponding rest-frame $`A_\mathrm{R}`$ index, derived out to $`z0.58`$ (Table 1). $`A_\mathrm{R}`$ is less sensitive than the rest-frame $`B`$-band to massive young stars. Thus, tidally induced features in mergers (e.g., Jogee, Kenney, & Smith 1999), can lead to large values of $`A_\mathrm{R}`$. The presence of both large $`A_\mathrm{B}`$ and $`A_\mathrm{R}`$ in a large fraction of the starbursts suggests that a significant part of their SF is tidally triggered. High asymmetry values have also been previously noted in a local sample of infrared luminous starburst galaxies (Conselice et al 2000). Visual inspection of our images confirms this claim as the majority of the systems with large $`A_\mathrm{B}`$ ($`>`$0.4) and $`A_\mathrm{R}`$ ($`>`$0.35) indeed appear to be optically disturbed or/and to have nearby companions.

## 4 Summary

From a comparative study of optically-selected starburst galaxies out to $`z1`$ and a control sample of normal (early-type E/Sa and late-type) galaxies, we find:

(1) The early-type systems have rest-frame $`B`$ concentration indices and AGN fraction ($`>`$25 %) which are significantly higher than those of late-type spirals and optically-selected starbursts at $`z1`$ . The results are generally consistent with the idea that early-epoch ($`z1`$) gas-rich dissipative processes (e.g., major mergers) have likely played an important role in developing large central concentrations in early-type E/Sa galaxies. The larger AGN fraction qualitatively supports scenarios leading to a concurrent growth of central black holes and bulges in some early merger events.

(2) The lower $`C_\mathrm{B}`$ and AGN fraction in the majority of $`z1`$ optically-selected starbursts, compared to early-type systems, suggest that either the starbursts and early types are different in nature (being respectively disk and bulge dominated), or/and are in different evolutionary phases (such that some of the starbursts in major mergers evolve into early-types).

(3) The $`z1`$ optically-selected starbursts, on average, have larger asymmetry indices than normal galaxies, both in rest-frame $`B`$ and $`R`$-band. This suggests that a significant fraction of the starbursts may be tidally triggered.

Support for this work was provided by NASA through grant GO09583.01-96A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, under NASA contract NAS5-26555.