text:
# Black Hole fragmentation and holography

## I Introduction

One of the most striking novelties in the research on gravitational physics is the possibility that gravity in $`d`$ dimensions could be described by a local field theory in $`d1`$ dimensions Maldacena:1997re ; Witten:1998qj ; Gubser:1998bc . The theoretical evidence for such holographic description of gravity is mounting. Indications that holography could be a fundamental feature of the gravitational interaction come from different directions: string theory, black hole physics, cosmology 'tHooft:gx ; Susskind:1994vu ; Verlinde:2000wg  (For a recent review see Bousso:2002ju ). A particularly interesting output of these investigations has been the formulation of stringent holographic bounds for the entropy of a system occupying a given region of space Susskind:1994vu ; Bousso:1999xy .

An explicit realization of the holographic principle has been found only in particular cases, essentially for anti- de Sitter (AdS) (and de Sitter) gravity, the so-called anti de Sitter/conformal field theory (AdS/CFT) correspondence Maldacena:1997re ; Witten:1998qj ; Gubser:1998bc . A general way to explicit realize the holographic principle for generic gravitational systems, in particular for gravity in asymptotically flat spacetimes, is still lacking. In particular, it is not clear if the realizations of the holographic principle take always the form a correspondence between $`d`$-dimensional gravity and a field theory in $`d1`$ dimensions, or if there could be some alternative, still unknown, realization of it. If the holographic principle has to be considered a genuine feature of every quantum theory of gravity, one could explain our lacking in understanding the holographic principle as a lacking in understanding quantum gravity. However, there are strong indications that holography is a feature of gravity that appears and therefore should be explained, already at the semiclassical level. The Bekenstein-Hawking area law for the black hole entropy is the most striking example of holographic behavior of a gravitational system that has to be explained already at the semiclassical level.

An alternative strategy one can use in this context, is to explore the similarities and the differences between gravity and local field theories in order to check at a fundamental level the possibility of finding correspondences between the two classes of theories. This approach can be very powerful. One nice example, discussed in almost every introducing paper on the holographic principle, is the scaling behavior of the entropy as a function of the volume of the system for a local field theory compared to that of a black hole. For a local field theory the entropy is an extensive quantity, it scales as the volume of the space. On the other hand the entropy of a black hole scales as the area of the horizon. This simple fact enables one to conclude that the correspondence between gravity and field theory, if it exists, must be holographic.

In this paper we will focus on an other aspect of the relationship between gravity and local field theory, namely on the dependence for composite systems of the entropy from the energy. Working in the microcanonical ensemble the entropy /energy relation for a free field theory in $`d`$-dimensions is given by

$$SE^{(d1)/d},$$
(1)
which considering two arbitrary excitations with energy $`E_1,E_2`$ satisfies the inequality

$$S(E_1+E_2)<S(E_1)+S(E_2),$$
(2)
Conversely asymptotically flat black holes of masses $`M_1,M_2`$ in $`d4`$ dimensions satisfy

$$S(M_1+M_2)>S(M_1)+S(M_2),$$
(3)
i.e fragmentation of a black hole of mass $`M_1+M_2`$ into two smaller black holes of masses $`M_1`$ and $`M_2`$ is entropically not preferred. Assuming the existence of a gravity/field theory correspondence, one has to identify the black hole masses $`M_1`$, $`M_2`$ as excitations $`E_1`$, $`E_2`$ of the field theory. It follows that Eq. (2) contradicts Eq. (3) and consequently that the assumed correspondence gravity/field theory cannot be true. Also skipping this problem, i.e assuming the existence of black hole solutions satisfying inequality (2) rather then (3), one is faced by an other problem. At first sight Eq. (2) seems incompatible with every entropy bound because a black hole could always increase its entropy by fragmentating into smaller black holes.

In this paper we will analyze in detail for various classes of black holes in various dimensions the validity of the relation (3). We will show that Eq. (3) holds true only for asymptotically flat black holes. Black holes with different asymptotic behavior (for instance AdS black holes) may satisfy Eq. (2). For the two-dimensional (2D) case we will be able to give a complete characterization of the entropy for composite black holes in terms of the asymptotic behavior of the solution. We will also show the existence of black hole solutions with mixed behavior, i.e solutions satisfying either Eq. (2) or Eq. (3) when the masses $`M_1`$ and $`M_2`$ are above or below some critical value $`\stackrel{~}{M}_0`$. We will argue that this mixed behavior as a natural counterpart in a field theory with a mass gap. Finally, we also discuss the compatibility of Eq. (2) with entropy bounds.

## II Asymptotically flat black holes

Let us first consider the $`d`$-dimensional Schwarzschild black hole ($`d4`$)

$$ds^2=\left(1\frac{k_dGM}{r^{d3}}\right)dt^2+\left(1\frac{k_dGM}{r^{d3}}\right)^1dr^2+r^2d\mathrm{\Omega }_{d2}^2,$$
(4)
where $`M`$ is mass and $`k_d=16\pi /(d2)\mathrm{\Omega }_{d2}`$, $`\mathrm{\Omega }_{d2}`$ being the volume of the unit $`S^{d2}`$ transverse sphere. The black hole entropy is given by the area law

$$S=\frac{A}{4G}=\frac{\mathrm{\Omega }_{d2}}{4}(k_dL_pM)^{(d2)/(d3)},$$
(5)
where $`L_p=G^{1/(d2)}`$ is the Planck length.

In order to discuss the inequality (3) on a physical ground, we will make here and throughout the all paper the following assumptions: 1) The theory admits, at least in some approximation, multi black hole solutions and 2) The gravitational potential energy of the multi black hole configuration can be neglected with respect to the black hole masses. Assumption 1) is necessary in order to give to a multi black hole configuration a precise meaning whereas assumption 2) assures us that the multi black hole configuration can be treated as a composite system with zero binding energy, whose mass and entropy is simply the sum of that of the elementary constituents.

Using Eq. (5) one easily find that Eq. (3) is satisfied. As expected, for Schwarzschild black holes the entropy is maximized by the single black hole configuration with mass $`M_1+M_2`$. The same holds true for asymptotically flat charged black holes. Considering for instance the four-dimensional Reissner-Nordstrom solution one finds for the entropy

$$S=\left(ML_p+\sqrt{M^2L_p^2Q^2}\right)^2,$$
(6)
where $`Q`$ is the electric charge of the black hole. Again, from Eq. (6) it follows $`S(M_1+M_2)>S(M_1)+S(M_2)`$. There is one simple argument that can be used to argue that inequality (3) holds in general for asymptotically flat black holes at least when $`ML_p`$ is much bigger then the (eventually present) black hole charges. The area law gives $`S(r_+/L_p)^{d2}`$, where $`r_+`$ is the horizon radius. For asymptotically flat black holes and when $`ML_p>>Q_i`$, where $`Q_i`$ are the charges associated with the black hole, the gravitational potential is dominated by the Newtonian term so that we have $`r_+L_p^{(d2)/(d3)}M^{1/(d3)}`$. It follows $`S(L_pM)^{(d2)/(d3)}`$,

## III  AdS black holes

The argument presented at the end of the previous section does not apply to non-asymptotically flat black holes. For 2D black holes a general discussion we will presented in the next section. Here we will discuss only the most interesting case, namely asymptotically AdS black holes. The entropy-mass relation for the $`d`$-dimensional Schwarzschild-anti de Sitter black hole ($`d4`$),

$$ds^2=\left(1+\lambda ^2r^2\frac{k_dGM}{r^{d3}}\right)dt^2+\left(1+\lambda ^2r^2\frac{k_dGM}{r^{d3}}\right)^1dr^2+r^2d\mathrm{\Omega }_{d2}^2,$$
(7)
is rather complicated. However a simple formula can be found for black holes with $`ML_p^{d2}\lambda ^{d3}>>1`$. In this latter case we have

$$S\left(\frac{M}{L_p\lambda ^2}\right)^{(d2)/(d1)}.$$
(8)
Differently from the asymptotically flat case, we see that now inequality (2) is satisfied. For AdS black holes ( at least for those with large enough mass) fragmentation is entropically preferred. This fact is perfectly consistent with the existence of a AdS/CFT correspondence between $`d`$-dimensional AdS gravity and $`d1`$-dimensional conformal field theory. The $`d`$-dimensional black hole entropy (8) reproduces correctly, after the $`E=M`$ identification, the entropy (1) for a free field theory in $`d1`$-dimensions.

Let us now discuss in some detail the three-dimensional (3D) case, the Banados-Teitelboim-Zanelli (BTZ) black hole Banados:wn . This case is very instructive not only because it is possible to solve exactly the inequality (2) but also because it clarifies the role played by the black hole ground state in our considerations. The BTZ black hole solution with zero angular momentum is

$$ds^2=(\lambda ^2r^28GM)dt^2+(\lambda ^2r^28GM)^1dr^2+r^2d\varphi ^2,$$
(9)
where the black hole mass $`M`$ is defined with reference to the $`M=0`$ black hole ground state, $`ds^2=\lambda ^2r^2dt^2+\lambda ^2r^2dr^2+r^2d\varphi ^2`$. The entropy is given by $`S=2\pi r_+/4G=(\pi /\lambda )\sqrt{2M/G}`$ so that Eq. (2) is satisfied for every value of the mass $`M`$. On the other hand, if the black hole mass is defined with reference to the full, geodetically complete AdS spacetime $`ds^2=(1+\lambda ^2r^2)dt^2+(1+\lambda ^2r^2)^1dr^2+r^2d\varphi ^2`$, the mass spectrum looks rather different. The $`M=0`$ ground state is separated from the continuos part of the spectrum with $`M1/8G`$ by a mass gap. For $`M1/8G`$ the entropy is given by $`S=(\pi /2\lambda G)\sqrt{8GM1}`$. Taking for simplicity $`M_1=M_2=M`$ we see that Eq. (2) is satisfied only for $`M3/16G`$, whereas for $`1/8GM3/16G`$ we have $`S(M_1+M_2)S(M_1)+S(M_2)`$. We will come back to this point in Sect. V , where will argue that this behavior is typical for a spectrum with a mass gap.

## IV Two-dimensional black holes

In $`d3`$ spacetime dimensions it is very difficult to formulate general criteria that enable us to decide if a black hole satisfies either Eq. (2) or Eq. (3). These criteria can be found for 2D black holes. The 2D case is interesting for several reasons. Two-dimensional gravity supports a realization of the AdS/CFT correspondence Cadoni:1998sg . In two spacetime dimensions we can formulate entropy bounds Mignemi:2003uz . Moreover, 2D black holes can be used to describe black holes in higher dimensions. Two dimensional black holes arise as effective description of the near-horizon, near-extremal behavior of $`d`$-dimensional charged black hole and branes. Every $`d`$-dimensional spherically symmetric solution can be described, after dimensional reduction, by a 2D dilaton gravity model.

The generic 2D dilaton gravity model (For a recent review see Grumiller:2002nm ) can be completely characterized by a dilaton potential $`V(\mathrm{\Phi })`$, the action for the model being given by $`A=(1/2)d^2x\sqrt{g}(\mathrm{\Phi }R+\lambda ^2V(\mathrm{\Phi }))`$. The general 2D black hole solution takes the form

$$ds^2=\left(J(\mathrm{\Phi })\frac{2M}{\lambda }\right)dt^2+\left(J(\mathrm{\Phi })\frac{2M}{\lambda }\right)^1dr^2,\mathrm{\Phi }=\lambda r,$$
(10)
where $`J(\mathrm{\Phi })=V(\mathrm{\Phi })`$ and $`M`$ is the black hole mass. To be sure that Eq. (10) describes a black hole we will take $`M0`$, $`\mathrm{\Phi }>0`$ and $`J(\mathrm{\Phi })>0`$ a strictly increasing function of $`\mathrm{\Phi }`$, (we will therefore have $`J(\mathrm{})=\mathrm{}`$). The temperature and entropy associated with the black hole are given by

$$S=2\pi \mathrm{\Phi }_h,T=\frac{\lambda }{4\pi }V(\mathrm{\Phi }_h),$$
(11)
where $`\mathrm{\Phi }_h=J^1(M)`$ is the value of the dilaton at the black hole horizon.

Using Eq. (10) one easily finds

$$S=2\pi J^1(M).$$
(12)
One can now find sufficient conditions to be imposed on the function $`J(\mathrm{\Phi })`$ such that Eq. (2) is satisfied $`M_1,M_2`$. Let us first show that if both conditions (the prime denotes derivation with respect to $`M`$)

$`(J^1)^{\prime \prime }(M)<0,M,`$ (13)
$`J^1(M=0)=0,`$ (14)
are satisfied then

$$J^1(M_1+M_2)<J^1(M_1)+J^1(M_2),$$
(15)
which, owing to Eq. (12) implies that the inequality (2) for the entropy is also satisfied. The convexity condition (13) implies for every positive $`\widehat{M},M_1,`$ $`(J^1)^{}(\widehat{M}+M_1)<(J^1)^{}(\widehat{M})`$, from which it follows $`J^1(\widehat{M}+M_1+M_2)+J^1(\widehat{M})<J^1(\widehat{M}+M_1)+J^1(\widehat{M}+M_2)`$. Evaluating the previous inequality at $`\widehat{M}=0`$ and using condition (14) one easily recovers Eq. (15).

Analogously, we can show that if conditions

$`(J^1)^{\prime \prime }(M)>0,M,`$ (16)
$`J^1(M=0)=0,`$ (17)
hold then the inequality (3) is satisfied for every $`M_1,M_2`$. If $`(J^1)^{\prime \prime }=0`$, identically, then the black hole solution of the model will satisfy $`S(M_1+M_2)=S(M_1)+S(M_2)`$. If $`(J^1)^{\prime \prime }`$ changes sign, we cannot make any definite statement about inequalities (2) or (3).

Any function $`J^1(M)`$ which is everywhere convex (concave) and strictly growing must diverge for $`M\mathrm{}`$ lesser (faster) then $`J^1=M`$. It follows that $`J(\mathrm{\Phi })`$ must diverge for $`\mathrm{\Phi }\mathrm{}`$ faster (lesser) then $`J=\mathrm{\Phi }`$. We have reached an important result: for models satisfying the required criteria about the derivatives of $`J^1`$ and $`J(0)=0`$, fragmentation of 2D black holes can be or not be entropically preferred depending on the asymptotic behavior of the function $`J`$. The black solutions of the 2D dilaton gravity model, characterized by a function $`J`$, which asymptotically diverges faster (lesser) then $`\mathrm{\Phi }`$, will satisfy Eq. (2) (Eq. (3)).

When $`J^1(0)0`$ but condition (13) still holds we will have a mixed behavior. The black hole will satisfies Eq. (2) for $`M>\stackrel{~}{M}_0`$ and Eq. (3) for $`M<\stackrel{~}{M}_0`$, where $`\stackrel{~}{M}_0`$ is some threshold mass. Let us first notice that only the case $`J^1(0)=\mathrm{\Phi }_0<0`$ has physical relevance. Being $`J^1`$ a strictly growing function, $`\mathrm{\Phi }_0>0`$ implies $`J^1(M_0)=0`$ with $`M_0`$ negative, i.e the presence of negative masses in the spectrum. Conversely, $`\mathrm{\Phi }_0<0`$ implies

$$J^1(M_0)=0,$$
(18)
with $`M_0`$ positive, i.e the black hole spectrum is limited from below by the extremal, non vanishing value of the mass $`M_0`$. In general the extremal mass $`M_0`$ is simply related with the threshold mass $`\stackrel{~}{M}_0`$. Using Eqs. (12) and (18) one easily realizes that the entropy of the extremal state is zero, $`S(M_0)=0`$. Depending on the behavior of the dilaton potential at $`M_0`$ we can have two cases. 1) $`J^1(M_0)=0`$ and $`V(M_0)=0`$. The extremal state has zero entropy and from Eq. (11) also zero temperature. 2) $`J^1(M_0)=0`$ and $`V(M_0)0`$. The extremal state has zero entropy but nonvanishing temperature. In the next section we shall see that in both cases this behavior can be explained with the presence of a mass gap.

Let us now give some examples to illustrate our general results. The simplest model satisfying both conditions (13),(14) is the Jackiw-teitelboim model, $`V=2\mathrm{\Phi }`$. We have $`J^1(M)=\sqrt{2M/\lambda }`$, which satisfies both $`(J^1)^{\prime \prime }<0`$ and $`J^1(0)=0`$. For the entropy we get $`S\sqrt{M}`$, from which Eq. (2) follows for every value of the mass. A more general model is given by Mignemi:1994wg

$$V=(h+1)\mathrm{\Phi }^h,h>1.$$
(19)
For $`h>0`$ the model satisfies the condition (13) and (14). The black hole entropy is given by

$$S(M/\lambda )^{1/(h+1)},$$
(20)
which as expected satisfies Eq. (2). This class of models contains, as particular cases, 2D gravity models arising as near-horizon limit of dilatonic zero-branes Cadoni:2001ew , black 3-branes Cadoni:2003vi  and heterotic string black holes Cadoni:1999gh . It is interesting to note that the entropy-energy relation (20) becomes that of a free field theory given by Eq. (1) identifying $`h=1/(d1)`$. This fact gives an hint about the possibility of finding a correspondence between these 2D dilaton gravity models and a free field theory.

For $`1<h<0`$ the model satisfies the conditions (16) and (17). The black hole entropy now satisfies Eq. (3). An important particular case is given by $`h=1/2`$, which describes the spherical dimensional reduction of the Schwarzschild black hole. For $`h=0`$ we have the CGHS model Callan:rs  . Eq. (19) gives $`(J^1)^{\prime \prime }=0`$ identically. The entropy depends linearly on the mass, so that, as expected, $`S(M_1+M_2)=S(M_1)+S(M_2)`$.

As an examples of a model fulfilling condition (13) but not (14) let us first consider the exponential potential $`V=\beta \mathrm{exp}(\beta \mathrm{\Phi })`$, with $`\beta >0`$. The black hole horizon is located at $`\lambda r=\mathrm{\Phi }=J^1(M)=(1/\beta )\mathrm{ln}(2M/\lambda )`$. $`J^1`$ diverges asymptotically lesser then $`M`$ but we have $`J^1(0)0`$. We therefore expect the presence of a threshold mass $`\stackrel{~}{M}_0`$ separating the two regions of the spectrum where Eq. (2), respectively Eq. (3) hold.

The black hole mass has an extremal value, $`MM_0=\lambda /2`$. However, $`J^1(M_0)=0`$ and $`V(M_0)0`$. The entropy and temperature of the black hole are

$$S=\frac{2\pi }{\beta }\mathrm{ln}\frac{2M}{\lambda },T=\frac{\beta }{2\pi }M.$$
(21)
The ground state has zero entropy but finite temperature $`T(M_0)=(\beta /4\pi )\lambda `$. Using Eq. (21) an taking for simplicity $`M_1=M_2=M`$, one can easily find that the inequality (2) is satisfied only for $`M>\stackrel{~}{M}_0=\lambda `$, whereas for $`\lambda /2<M<\lambda `$ Eq. (3) holds. Black hole fragmentation becomes entropically preferred for $`M`$ bigger then the threshold mass $`\stackrel{~}{M}_0`$. Notice that the threshold mass $`\stackrel{~}{M}_0`$ although of the same order do not coincide with the extremal mass $`M_0`$, we have $`\stackrel{~}{M}_0=2M_0`$.

As a second example of models with $`J^1(0)0`$ let us consider $`J(\mathrm{\Phi })=\mathrm{\Phi }^2+1`$. This is the 2D analogue of the 3D BTZ black hole discussed at the end of the previous section. Also in this case the black hole mass has an extremal value $`M_0=\lambda /2`$. However, now $`ST\sqrt{(2M/\lambda )1}`$, so that the extremal state has zero entropy and temperature, $`S(M_0)=T(M_0)=0`$. The threshold mass at which black hole fragmentation becomes entropically preferred is given by $`\stackrel{~}{M}_0=(3/2)M_0=(3/4)\lambda `$.

Until now we have considered a system in which the total energy is constant, i.e fragmentation of a black hole of mass $`M`$ into two smaller black holes of masses $`M_1+M_2=M`$. Let us now consider systems at constant volume. We will show that the energy is minimized when the conditions which maximize the entropy are satisfied.

We have to consider the black hole mass $`M`$ as a function of its radius $`R`$ and to solve the inequality

$$M(R_1+R_2)>M(R_1)+M(R_2).$$
(22)
Using Eq. (10) one finds $`M(R)=(\lambda /2)J(\mathrm{\Phi }_h)=(\lambda /2)J(R)`$. Eq. (22) becomes $`J(R_1+R_2)>J(R_1)+J(R_2)`$. Because $`(J^1)^{\prime \prime }<0`$ and $`J^1(0)=0`$ imply, respectively, $`J^{\prime \prime }>0`$ and $`J(0)=0`$ it follows that conditions (13), (14) for the function $`(J^1)`$ are equivalent to conditions (16), (17) for the function $`J`$. As a consequence, whenever conditions (13) and (14) are satisfied we have not only $`S(M_1+M_2)<S(M_1)+S(M_2)`$ but also $`M(R_1+R_2)>M(R_1)+M(R_2)`$. The process of fragmentation of a black hole maximizes the entropy if the total mass is constant and minimizes the mass if the total volume is constant. Conversely, if conditions (16) and (17) hold we have $`S(M_1+M_2)>S(M_1)+S(M_2)`$ but also $`M(R_1+R_2)<M(R_1)+M(R_2)`$. The configuration with maximal entropy and minimal mass is now given by the single black hole.

## V  Mass gap

It is evident that for what concerns the entropy of composite solutions 2D dilaton gravity models satisfying the conditions (13) and (14) of the previous section are very similar to free field theories. The mass of the gravitational solutions are naturally identified with the energy of the excitation in the field theory. The same is true for generic $`d`$-dimensional AdS solutions when the mass of the gravitational solution is much bigger then $`L_p^{d2}\lambda ^{d3}`$. We can ask ourselves if the 2D dilaton gravity models satisfying condition (13) but not (14) have also a field theory counterpart. More in general, one would like to find field theoretical counterparts of $`d`$-dimensional black hole solutions whose spectrum exhibits extremal $`M0`$ solutions. The most natural candidates are field theories with mass gaps. This is rather obvious for the 2D dilaton gravity models discussed at the end of Sect. IV (for instance the model with exponential dilaton potential), which are characterized by a ground state of mass $`M_0`$ with zero entropy and nonvanishing temperature. The only way to have a nondegenerate ground state at finite temperature is the presence of a mass gap in the spectrum. The energy gap $`M_0`$ must be of the order of the temperature. In fact for the model with exponential potential discussed in Sect. IV we have $`E_{gap}=M_0T(M_0)\lambda `$.

The relationship between black hole solutions and field theories with mass gaps seems to be more general. Crucial for the existence of such correspondence are both the asymptotic behavior of the metric and the presence of an extremal solution with $`M_00`$. To illustrate this correspondence let us consider a generic field theory whose spectrum has a mass gap of energy $`E_0`$ separating the $`E=0`$ state from the continuos part of the spectrum $`EE_0`$. We will also assume that for $`E>>E_0`$ the spectrum is that of a generic free field theory and that the entropy of the $`E=E_0`$ state vanishes. For $`E>>E_0`$ the entropy/energy relation will be given by Eq. (1) which satisfies $`S(E_1+E_2)<S(E_1)+S(E_2)`$. On the other hand for $`EE_0`$ we have $`S=S(EE_0)`$ with $`S(E_0)=0`$. Thus, the inequality $`S(E_1+E_2E_0)>S(E_1E_0)+S(E_2E_0)`$ will be satisfied at least for $`E_1=E_2=E_0`$. Because for $`E>>E_0`$ the opposite inequality holds, this implies the existence of a threshold energy $`\stackrel{~}{E}_0`$ separating the two regimes in complete analogy with what happens for black hole solutions. This behavior is rather intuitive. For small excitations near $`E_0`$ the single state of energy $`E=E_1+E_2`$ has more degeneracy then the two states of energies $`E_1`$ and $`E_2`$ because the states below the gap do not contribute. For excitations of energy $`E>>E_0`$ the contribution of the gap is irrelevant and the entropy is dominated by the contribution coming from the continuos part of the spectrum.

## VI  Entropy bounds

The holographic principle in its usual formulation puts an upper bound to the entropy, hence to the amount of information, which can be stored in given region of space, $`SA/4G`$, where $`A`$ is the area of the surface enclosing the region. At fixed volume the entropy cannot increase by splitting a system into parts, i.e $`S(V_1+V_2)S(V_1)+S(V_2)`$. This inequality is satisfied not only for extensive systems, when the entropy scales as the volume of the system , but also when the entropy scales as the area of the surface enclosing the system. An holographic bound cannot be violated by fragmentation. Hence the holographic principle is compatible with both Eqs. (2) and (3). This is true for the holographic bound but it is no necessarily true for entropy bounds of different type, which can make reference not only to the volume the system (or to the area enclosing it) but also to its energy. For instance the Bekenstein bound $`S2\pi ER`$ Bekenstein:jp  refers not only to the linear size $`R`$ but also to the energy $`E`$ of the system.

In principle one could also try to formulate entropy bounds, which make reference to nothing but the energy of system. This could be done following the same line of reasoning one uses to formulate the Bekenstein-Hawking bound. In this latter case one takes a system in a given region of space and finds that its entropy is bounded by the area of the surface enclosing the region. More in detail, when one pumps energy in the system keeping its spatial extension constant, the entropy increases until the energy equals the mass pertaining to a black hole fitting in the region. A black hole forms and the bound $`SA/4G`$ is saturated.

Instead, one could take a system with a constant amount of energy $`M`$ and look for a bound on the associated entropy. Spreading the system over regions of space with ever decreasing volumes the entropy will increase until the corresponding Schwarzschild radius is reached. A black hole forms and a entropy bound $`S\pi GM^2`$ is found. The problem is that this bound cannot be universal. Although, the Schwarzschild black hole cannot increase its entropy by fragmentating, we have seen that there exist many black hole solution that can increase their entropy trough fragmentation. For this simple reason an entropy bound which refers only to the energy of the system cannot be universal.

## VII conclusions

In this paper we have discussed the holographic principle and entropy bounds for composite gravitational systems. We have shown that with respect to fragmentation black hole solutions exhibit three different types of behavior. The entropy may decrease, increase or have a mixed behavior. In the latter case there will exist a threshold mass separating the regime where fragmentation is entropically preferred from the regime where it is not. Moreover, for 2D black holes we have been able to find a complete characterization of the entropy behavior, in the form of sufficient conditions imposed on the function $`J`$, which defines the 2D gravitational model.

Our results have a strong impact on the way one can realize the holographic principle. In fact only for gravity theories satisfying relation (2) one can hope to find a realization of the holographic principle in terms of a correspondence gravity/field theory. For gravity theories satisfying (3) this correspondence cannot be realized as a gravity /field theory duality and must necessarily have an alternative, still unknown, form. Because asymptotically flat black holes, and in particular the Schwarzschild black hole, belongs to this latter case we can exclude the existence of a correspondence between asymptotically flat black holes and a field theory.

On the other hand we have seen that the class of gravity theories satisfying Eq. (2) is, at least in the 2D case, rather broad. It contains not only AdS gravity, for which the gravity/ field theory duality is well established, but also other models, for instance those with a exponential dilaton potential. The possibility of finding a gravity/ field theory duality for these theories is an open question, which deserves further investigations.

An other point of interest is the existence of black holes with mixed behavior, i.e exhibiting a transition from a regime where Eq. (2) holds to a regime where instead Eq. (3) is satisfied. We have argued that the most natural candidates for holographic duals of this black holes are field theories with a mass gap. Finally, our results seem to exclude the possibility of formulating an entropy bound only in terms of the energy of a system: the existence of black holes satisfying Eq. (2) will always allow to violate the bound by black hole fragmentation.

###### Acknowledgements.

We thank S. Mignemi for discussions and valuable comments.


format: text

text:
# Aging in short-ranged attractive colloids: A numerical study

## I Introduction

In the last years, the physics of systems characterized by short-ranged attractive inter-particle interactions — short in comparison to the particle size — has been the focus of several investigations. A good example of this type of systems are colloidal solutions, since the intensity and the range of the potential can be finely controlled by changing the solvent properties or the chemistry of the dispersed particles. The short range of the attractive interaction produces thermodynamic features which are not observed in atomic or molecular system (where the interaction range is always long ranged)Pusey (1991); A.P.Gast et al. (1983); Meijer and Frenkel (1991); Anderson and Lekkerkerker (2002). Recently, theoretical, numerical and experimental studies have focused on the dynamical properties of these systems. One of the most astonishing discovery is that, at high density, the metastable liquid is characterized by a non monotonic temperature dependence of the characteristic structural times: the dynamics slows down not only upon cooling (as commonly observed in molecular systems), but also upon heating. The slowing down upon heating can be so intense that a novel mechanism of arrest takes place high $`T`$ Sciortino (2002). In the high density part of the phase diagram, for sufficiently short-ranged potentials, a re-entrant liquid-glass line is observed and a liquid phase emerges between two glasses. In other words, moving along a constant density path, it is possible to pass from a glass phase to a liquid phase and then again in a new glass phase just by progressively lowering the temperature. The two glasses are named “attractive glass” and “repulsive glass” due to the different mechanisms, respectively attraction and excluded volume, which control the dynamical slowing down NOT . In the re-entrant liquid region, between the attractive and repulsive glasses, dynamics is very unusual. Correlation functions show a subtle logarithmic decay, while the particle mean squared displacement follows a power-law in time.

The complex dynamics in short-ranged attractive systems has been interpreted within Mode Coupling Theory (MCT) Götze (1991). Theoretical predictions do not depend on the detailed shape of the attractive potential and have been confirmed for many different modelization of the attractive tail. Fabbian et al. (1999); Bergenholtz and Fuchs (1999); Dawson et al. (2001); Foffi et al. (2002a); Götze and Sperl (2003). A set of simulationsPuertas et al. (2002); Foffi et al. (2002b); Zaccarelli et al. (2002a); Puertas et al. (2003) and experiments Mallamace et al. (2000); Pham et al. (2002); Eckert and Bartsch (2002); Chen et al. (2003a, b); Eckert and Bartsch (2003) have confirmed the existence of such complex dynamics.
In numerical work based on the Square Well System (SWS), one of the simplest models studied, it has been shown that a diffusivity maximum in temperature appears when the width of the square well is short enough compared to the hard core diameter Zaccarelli et al. (2002a). The calculation of the iso-diffusivity lines in equilibrium Foffi et al. (2002b) confirms the reentrant behavior and the novel dynamics which take place in the reentrant region. Along an iso-diffusivity path, the decay of the density-density correlation functions, which at high temperature can be well described by a stretched exponential, becomes more and more uncommon upon cooling, showing at low T signs of logarithmic decayPuertas et al. (2002); Zaccarelli et al. (2002a); Sciortino et al. (2003).
Previous numerical work has focused on equilibrium properties and on the dynamics near the predicted glass-glass transition Puertas et al. (2002); Foffi et al. (2002b); Zaccarelli et al. (2002a, 2003a). Here we report a study of the non-equilibrium dynamics. Starting from an equilibrium state in the re-entrant fluid region, the system is quenched both at low and high temperature. We are then able to explore — within the same model — the aging dynamics for the attractive and repulsive glass, to find out if the differences which characterize the dynamics in the glassy state show up also in the out-of-equilibrium evolution. We discover that, while the aging dynamics of the repulsive glass is similar to the one observed in atomic and molecular liquidsKob and Barrat (1997); Sciortino and Tartaglia (2001a), the aging dynamics of the attractive glass shows novel unexpected features. We also report comparison with a recent experimental study of aging in short ranged colloidal systems Pham et al. (2003)

## II Simulation details

We investigate a system that has been extensively studied earlier in the equilibrium regime, a binary square well mixtureZaccarelli et al. (2002a). The binary system is a $`50\%`$-$`50\%`$ mixture of $`N=700`$ particles. The two species (labeled $`A`$ and $`B`$) are characterized by a diameter ratio $`\delta =\sigma _A/\sigma _B=1.2`$. Masses are chosen to be equal and unitary. The attraction is modeled by a square well interaction defined according to:

$$V_{ij}(r)=\{\begin{array}{cc}\mathrm{}\hfill & \text{ }r_{ij}<\sigma _{ij}\hfill \\ u_0\hfill & \sigma _{ij}<r_{ij}<\sigma _{ij}+\mathrm{\Delta }_{ij}\hfill \\ 0\hfill & r_{ij}>\sigma _{ij}+\mathrm{\Delta }_{ij}\hfill \end{array}$$
(1)
where $`\sigma _{ij}=(\sigma _i+\sigma _j)/2`$, $`i,j=A,B`$ and $`\mathrm{\Delta }_{ij}`$ has been chosen in such a way that $`ϵ_{ij}=\mathrm{\Delta }_{ij}/(d_{ij}+\mathrm{\Delta }_{ij})=0.03`$. We choose $`k_B=1`$ and the depth of the potential $`u_0=1`$. Hence $`T=1`$ corresponds to a thermal energy equal to the attractive well depth. The diameter of the small specie is chosen as unity of length, i.e. $`\sigma _B=1`$. Density is parametrized in terms of packing fraction $`\varphi =(\rho _Ad_A^3+\rho _Bd_B^3)\pi /6`$, where $`\rho _i=N_i/L^3`$, $`L`$ being the box size and $`N_i`$ the number of particles for each species. Time is measured in units of $`\sigma _B(m/u_0)^{1/2}`$. The unit has been chosen consistently with Ref. Zaccarelli et al. (2002a), where the same system has been studied in equilibrium. A standard MD algorithm has been implemented for particles interacting with SWS potentials Rapaport (1997). Between collisions, particles move along straight lines with constant velocities. When the distance between the particles becomes equal to the distance where $`V(r)`$ has a discontinuity, the velocities of the interacting particles instantaneously change. The algorithm calculates the shortest collision time in the system and propagate the trajectory from one collision to the next one. Calculations of the next collision time are optimized by dividing the system in small subsystems, so that collision times are computed only between particles in the neighboring subsystems.
With the present choice of parameters, the formation of crystallites during the simulation run is not observed. As compared to the monodisperse case, the slight asymmetry in the diameters allows us to follow the development of the slow dynamics over several orders of magnitudeFoffi et al. (2002b). We choose to work at the packing fraction $`\varphi =0.608`$, being this density in the middle of the re-entrant region of the phase diagram, as shown in previous calculations Zaccarelli et al. (2002a). At this density we prepared 60 independent configurations, equilibrated at the initial temperature $`T_i=0.6`$. Subsequently we quenched these independent configurations to the desired final temperatures, and then we followed the evolution in time at constant temperature. The characteristic time of the thermostat as been chosen to be much smaller than the structural relaxation time and such that the system may equilibrate within one time unit. Each of the 60 independent configurations has been quenched to $`T_f=1.2`$, i.e. in the repulsive glass, and to $`T_f=0.3`$, i.e. in the attractive glass and run up to $`t_f=210^5`$. The numerical protocol is illustrated in Fig. 1. Iso-diffusivity curves from Ref. Zaccarelli et al. (2002a) and the extrapolated ideal glass line from Ref. Sciortino et al. (2003) are plotted together with the quench path selected in this work. Since the system, after the quench, is out of equilibrium, all averages — represented in the text by the symbol $``$ — have been performed over the ensemble of the 60 initial independent configurations.

## III Results and Discussion

Numerical investigations of aging in glassy materials are usually accomplished by rapidly quenching the system from an equilibrium state to a new state, where the characteristic time scale of the structural relaxation is much longer than the observation time. When the time scale of the final state exceeds the characteristic time of the simulation, equilibrium cannot be reached within the simulation time. Out of equilibrium time translation invariance does not hold and the state of the system does not only depend from the external parameters, but also from the time elapsed from the quench, the so-called waiting time $`t_w`$. Indeed in the aging regime other notable effects emerge, i.e. violation of the fluctuation dissipation theorem Cugliandolo and Kurchan (1993); Parisi (1997); Latz (2000); Sciortino and Tartaglia (2001b); De Gregorio et al. (2002) and strong dependence on the story of the system Kob and Barrat (1997). Thus, in out-of-equilibrium conditions correlation functions will depend both on the time of observation $`t`$ and on the waiting time $`t_w`$, i.e. $`C_A(t_w,t)=A^{}(t_w)A(t_w+t)`$. Similarly, the one time properties (e.g. the static structure factor or the energy), will depend on the waiting time $`t_w`$.
The first quantity that we investigate is the average potential energy of the system $`U`$. At $`t_w=0^+`$, i.e. instantaneously after the quenches to respectively $`T_f=1.2`$ and $`T_f=0.3`$, the two systems possess the same potential energy. Then, with $`t_w`$, they start evolving toward different energy values. The situation is represented in Fig. 2. While the system quenched to $`T_f=1.2`$, i.e. in the repulsive regime, shows an increase in the potential energy (i.e. a progressive breaking of the interparticle bonding), the one quenched to $`T_f=0.3`$ shows the opposite behavior. In both cases, the time dependence of the potential energy shows a fast short-time ($`t_w<1`$) evolution (which reflects the equilibration of the fast degrees of freedom) followed by a much slower evolution, proper of the aging dynamics. The aging dynamics has detectable different trends at high and low temperatures. While at high temperatures a clear $`log(t_w)`$ dependence of $`U`$ is observed, at low $`T`$ the time dependence in $`log`$ scale is not linear and can be well fitted by $`UU_{eq}t_w^a`$, with $`a0.14`$. A similar form was successfully adopted by Parisi for a binary system of soft spheres in the aging regime Parisi (1997).

The time evolution of the aging structure can be followed by numerical calculation of the partial static structure factors, $`S_{\alpha \beta }(q,t_w)`$ defined as

$$S_{\alpha \beta }(q,t_w)=\varrho _\alpha ^{}(𝐪,t_w)\varrho _\beta (𝐪,t_w)/\sqrt{n_\alpha n_\beta }$$
(2)
where the partial density variables are defined as $`\varrho _\alpha (\stackrel{}{q})=_k\mathrm{exp}[i𝐪𝐫(t)_k^{(\alpha )}]/\sqrt{N}`$, and $`n_i=N_i/N`$ are the partial concentrations. In equilibrium, time translation invariance implies $`S_{\alpha \beta }(q,t)=S_{\alpha \beta }(q,0)`$, i.e. the static structure factor is a time independent function.
The inset of Fig. 3 shows $`S_{AA}(q,t_w)`$ for the quench to $`T_f=1.2`$, for waiting times spanning $`4`$ decades. Similarly to binary Lennard-Jones systems Kob and Barrat (1997), the variation of the static structure factor is particularly small, for all $`\alpha ,\beta `$ components. With increasing $`t_w`$ a very weak increase at all peaks of the structure factors is observed. This is shown in the main figure. Similarly for the $`T_f=0.3`$ case, the $`S_{\alpha \beta }(q,t_w)`$ changes are very small. Interestingly enough, for the first peak (and only there) the trend is reverted with respect to the previous case, i.e. the intensity of $`S_{AA}(q,t_w)`$ around the peak decreases on aging. This observation agrees with equilibrium studies. Indeed, for the square well system, an increase in temperature produces an increase in the first peak toward the hard sphere limit whereas lowering the temperature only promotes the structure factor oscillations at large wave-vectorZaccarelli et al. (2003b). It is perhaps interesting to recall that within MCT, the increase of $`S_{\alpha \beta }`$ at the first peak is responsible for the formation of the repulsive glass, whereas the long lasting oscillations in the structure factor generates the kinetic arrest in the attractive glass Dawson et al. (2001).

Next we study the intermediate scattering function, i.e. dynamic density-density correlators. For a binary mixture this quantity is defined by:

$$\varphi _{\alpha \beta }(q,t,t^{})=\varrho _\alpha (𝐪,t)^{}\varrho _\beta (𝐪,t^{})/S_{\alpha \beta }(q,t^{})$$
(3)
where the density variables $`\varrho _\alpha (𝐪,t)`$ have been defined before. In equilibrium $`\varphi _{\alpha \beta }(q,t,t^{})`$ will depend only on the difference $`tt^{}`$. When the system ages the density-density correlators is a function of $`t_w`$ and $`t`$. In this paper, we show $`\varphi _{\alpha \beta }(q,t,t^{})`$ as a function of $`t_w`$ and the time difference $`(tt_w)`$, i.e. $`\varphi _{\alpha \beta }(q,t_w,tt_w)`$. The density correlation is normalized by $`S_{\alpha \beta }(q,t_w)`$. The upper panel of Fig. 5 shows the $`t_w`$ dependence of $`\varphi _{AA}(q,t_w,tt_w)`$, with $`q\sigma _B=8.72`$, which lies between the first and the second peak in the static structure factor. On increasing $`t_w`$, the decay of the correlation functions becomes slower and slower. The time needed to lose memory of the starting configuration grows. This is a typical feature of aging systems and it has been encountered both in experiments Pham et al. (2003) and simulations Kob and Barrat (1997). The slowing down of the dynamics on increasing $`t_w`$ resembles the way the structural arrest emerges in supercooled liquids on approaching the glass transition. Indeed, in this case, upon cooling two characteristic time scales emerge, a fast decay around a plateau value and a slow decay to zero. The time duration of the plateau becomes longer, the closer the system is to the glass transition. On the other hand, in the aging regime the control parameter is the waiting time $`t_w`$, which plays a role analogous to that played by $`T`$ in equilibrium, i.e. the longer $`t_w`$ the longer is the time requested to lose memory. In the same way as temperature measures the distance from the glass transition in equilibrium liquids, $`t_w`$ measures the progressive thermalization of the structural degrees of freedom. In summary, with increasing $`t_w`$ the system develops a clear plateau, whose height is consistent with the value expected for the repulsive glass.
Similar results are shown in the lower panel of Fig. 5 for $`q\sigma _B=6.37`$, which is close to the maximum of the static structure factor. Two main differences emerge between the correlators at the the two wavevectors: (i) the relaxation time is longer for $`q\sigma _B=6.37`$, (ii) the plateau is higher. Both this issues follow, as expected, the trend of the equilibrium results, where at the maximum of the structure factor a maximum of the relaxation time and of the non-ergodicity parameter is encountered Foffi et al. (2003).
The upper panel of Fig. 6 shows the density-density correlation functions for the quench at $`T_f=0.3`$ for $`q\sigma _B=8.72`$. Similarly to what we found in the previous case, the time of decay of the correlators grows with the waiting time $`t_w`$. However the shape of the functions is drastically different, confirming that even the aging dynamics reflects the differences noted in the previous equilibrium studies. Two notable differences with the $`T_f=1.2`$ case are: (i) the strength of the $`\alpha `$-relaxation(i.e. the height of the plateau) in aging is much more intense, being now close to one ; (ii) the plateau does not significantly stretch in time (as shown in the inset) beyond $`tt_w1`$. Similar trends are found for a $`q`$-vector closer to the first peak of the structure factor, as shown in the lower panel of Fig. 6.

A recent study of the glass-glass transition has suggested the possibility that the attractive glass is significantly destabilized by hopping processes Zaccarelli et al. (2003a), which in the present context can be associated with the breaking of the interparticle bonds. The data shown in the inset of Fig.6 support such interpretation, since no clear plateau can be observed for $`t_w>1`$.

A very interesting feature of the correlation functions in the aging regime is the property of scaling with waiting time. Theoretical studies suggest that correlators at different waiting time can be rescaled by the equation:

$$\varphi _{\alpha \beta }(q,t_w,tt_w)=\varphi _{\alpha \beta }^{st}(q,t)+\varphi _{\alpha \beta }(q,\frac{h(tt_w)}{h(t_w)})$$
(4)
where $`\varphi _{\alpha \beta }^{st}(q,t)`$ accounts for the short time Cugliandolo . The $`h(t)=t`$ case is the so-called simple aging case. For a binary mixture of Lennard-Jones spheres the system follows (4) with $`h(t)=t^\alpha `$ and $`\alpha 0.88`$ Kob and Barrat (1997).
We tried to apply this scaling relation to our system. We focus on the low temperature quench, since the high temperature one shows a phenomenology, which is similar to the Lennard-Jones case. In the case of aging of $`T_f=0.3`$ no clear scaling of the correlation functions can be performed. In Fig 7a we show the crude tentative of data scaling. Correlators have been shifted in time to superimpose the correlators in the region where $`0.7<\varphi _q<0.9`$. Fig. 7b shows the same data, this time shifted to maximize superposition for $`0.35<\varphi _q<0.5`$. In Fig. 8 the scaling coefficients $`t_r`$ used to obtain the scaling in Fig. 7 are shown as function of the waiting time $`t_w`$. It can be seen that even if they could be fitted with power laws (with an exponent $`\alpha 0.38\pm 0.1`$ different from the exponent found for the Lenard-Jones case), the collapse of the different $`t_w`$ curves is not good.

It is worth to make a comparison of our numerical results with recent data measured by Pham et.al Pham et al. (2003), for a system of colloids with depletion interaction. In their experiment, they measure the density-density correlation function in the aging regime for two system at different polymer concentration, corresponding respectively to conditions where the repulsive and the attractive glass are expected. To better compare with Pham et. al data, we calculated the density dynamical structure factor (irrespective of the particle type) for the same $`q`$-vector of the cited experiment. Results are shown in Fig. 9 and Fig. 10 where the $`T_f=1.2`$ and $`T_f=0.3`$ quenches are respectively plotted together with the experimental results from Ref. Pham et al. (2003) (lower panel). Slight differences in the plateau values can be ascribed to differences in the actual composition of the simulated (binary mixture) and experimental (slightly polydisperse system). The agreement between the numerical and the experimental data is impressive.

Another important two-time observable is the mean squared displacement (MSD). In aging, it is defined as $`|𝐫(tt_w)𝐫(t_w)|^2`$, where $`𝐫(t)`$ is the position of the particle at time $`t`$ and the average $``$ is performed over all the particles and the independent configurations. Here we focus on the total MSD, i.e. evaluated with no distinction between particles of the two species.

Fig. 11a shows the MSD for $`T_f=1.2`$. In analogy with the case of the density correlators, at short times there is no significant $`t_w`$-dependence. At larger times the diffusive process is visible. When the waiting time becomes large, a plateau starts to emerge at intermediate times, between the ballistic (short time) and the diffusive (long time) regime. The duration of the plateau grows as $`t_w`$ increases. As above, this behavior is similar to the one observed in previous simulations of a binary mixture of Lenard-Jones particles Kob and Barrat (1997).

This three time-region behavior (ballistic, plateau, diffusive) is generally explained in term of the so called cage effect. When repulsive interactions start to be dominant the system slows down since each particle starts to get trapped in a shell of first neighbor particles. The length of the plateau indicates the typical time that a particle needs to escape this cage. The value of the plateau stands at the typical size of a cage, i.e. around $`10\%`$ of the diameter in agreement with Lindemann’s melting criterion Ashcroft and Mermin (1976). For $`T_f=0.3`$, the MSD shows a completely different behaviour, (Fig.11b). Its evolution can again be divided in three time-regions, but now the short time (ballistic) and the long time (diffusive) regions bracket a long intermediate region where the MSD appears to grow as a power-law in time, with an effective exponent which decreases with $`t_w`$. Differently from the previous case, no evident sign of a plateau is present, except for a hint of a plateau at MSD$`=0.001`$ which starts to develop for large waiting times for $`tt_w0.5`$. This value is consistent with a displacement of the order of the width of the potential well, and can be interpreted as a result of particles bonding in the attractive well, i.e. as some sort of attractive cages. However, activated bond-breaking processes appear to destroy such a confinement on time scales longer than $`tt_w1`$.

## IV Conclusions

In systems with short-ranged attractive potentials, an efficient competition between attraction and excluded volume generates a highly non trivial equilibrium dynamics and two kinetically distinct glasses. In the “attractive glass”, the particle localization is controlled by the bonding distance, while in the“repulsive glass” localization is controlled by the neighbor location. This offers the unique possibility of studying two different aging dynamics within the same model. In this article, we have realized this numerically for a model system which shows both types of disordered arrested structures, the square well potential. We have found that the aging dynamics is clearly different for the two glasses, both in the waiting time evolution of the one-time quantities, and in the time evolution of the two-times quantities.

The results reported in this article refer to a simple model system which capture the essence of the physics introduced by the short-range attraction. The square well potential is indeed a satisfactory representation of the coarse-grained particle-particle potential in colloidal system, for example in the presence of depletion interactions. Of course, dynamical features of the solvent constituents and of the solvent-particle interactions may significantly affect the particle dynamics. In this respect, the experimental validation of the equilibrium results is very valuable.

In particular, we have shown that the partial static structure factor evolves with $`t_w`$, showing opposite trends in the two cases, i.e. the first peak grows with $`t_w`$ for the repulsive glass while decreases for the attractive one. We have also shown that the decay of correlation functions becomes slower and slower on increasing $`t_w`$ in both cases, but with noteworthy differences in the strength of the relaxation (much more intense for the attractive glass case) and in the shape of the relaxation itself. The analysis of the MSD brings clear evidence that, while in the repulsive case there is the emergence of a plateau, whose time duration increases significantly with $`t_w`$, in the attractive glass case this does not happen. Indeed, in short-ranged attractive colloids, activated processes can be associated with thermal fluctuations of order $`u_0`$, which allow particles to escape from the bonds. These processes generate a finite bond lifetime and, at the same time, may destabilize the attractive glass. It is important to find out how the difference in the aging dynamics arising from the different localization mechanisms are affected by the presence of activated processes.

Finally, we have shown that a good agreement is found between experimental and numerical data, also in the aging dynamic.

###### Acknowledgements.

We acknowledge support from MIUR COFIN 2002 and FIRB and from INFM Iniziativa Calcolo Parallelo. We thank W.C.K. Poon and K.N. Pham for useful discussion and for brining their results to our attention.


text:
# Special Lagrangian submanifolds in the complex sphere

## Introduction

Special Lagrangian submanifolds may be defined as those submanifolds which are both Lagrangian (an order 1 condition) and minimal (an order 2 condition). Alternatively, they are characterised as those submanifolds which are calibrated by a certain $`n`$-form (cf \[HL\]), so they have the remarkable property of being area minimizing. Their study have received many attention recently since connections with string theory have been discovered. More particularly, understanding fibrations of special Lagrangian (possibly with singularities) in Calabi-Yau manifolds of (complex) dimension 3 is crucial for mirror symmetry (cf \[SYZ\],\[Jo\]). Since the pioneering work of Harvey and Lawson \[HL\], where this notion were introduced, several authors (\[Ha\],\[Jo\],\[CU2\]) have discovered many families of special Lagrangian submanifolds in the complex Euclidean space, however very few examples (cf \[Br\]) of such submanifolds are known in other Calabi-Yau manifolds, where this notion is naturally extended. Maybe the main reason is that such manifolds are somewhat rare, in particular the existence of compact, Calabi-Yau manifolds is a hard result of S.-T. Yau \[Y\] involving non explicit solutions to some PDE. However there exists some intermediary examples of non-flat, non-compact Calabi-Yau manifolds, maybe the simplest of them being the complex sphere.

In this paper we describe a family of Lagrangian submanifolds of a complex variety of $`^{n+1}`$ with a $`SO(n)`$-invariance. In the case of the complex sphere, we obtain among them new examples of special Lagrangian submanifolds, which are a kind of generalization of the Lagrangian catenoid discovered in \[HaLa\] (cf also \[CU1\]).

In \[CMU\], minimal Lagrangian submanifolds have been studied in the complex hyperbolic and complex projective spaces respectively. This is somewhat close to our situation for the following reason: these spaces are not Calabi-Yau manifolds, and so these submanifolds are not minimizing a priori ; however Y.-G. Oh has proved in \[Oh\] that a minimal Lagrangian submanifold of the complex hyperbolic space is stable and without nullity, a property necessary for being a minimizer.

The paper is organized as follows:

In Section 1 we give a definition of Calabi-Yau manifolds, we state some basic facts on them and describe the Calabi-Yau structure yielded on the complex sphere by the Stenzel metric. In Section 2 we describe a general class of Lagrangian submanifolds whith some $`SO(n)`$-invariance immersed in some complex submanifold of $`^{n+1}.`$ In the particular case of the complex sphere, we then check that they are also Lagrangian for the Calabi-Yau structure. In the last section we find in this class the special Lagrangian and describe them.

The nice argument about $`SO(n+1)`$-invariance at the end of next section is due to Dominic Joyce and Mu-Tao Wang suggested to me the Lagrangian Ansatz; it is my pleasure to thank them. The author, as a member of EDGE, Research Training Network HPRN-CT-2000-00101, has been supported by the European Human Potential Programme.

## 1 The Calabi-Yau structure on the complex sphere

We start with the following definitions — which are not the standard ones — proposed in \[Jo\]:

###### Definition 1

Let $`n2.`$ An almost Calabi-Yau manifold is a quadruple $`(X,J,\omega ,\mathrm{\Omega })`$ such that $`(X,J)`$ is a $`n`$-dimensional complex manifold, $`\omega `$ the symplectic form of a Kähler metric $`g`$ on $`X`$ and $`\mathrm{\Omega }`$ a non-vanishing holomorphic $`(n,0)`$-form.

$`(X,J,\omega ,\mathrm{\Omega })`$ is called a Calabi-Yau manifold if in addition

$$\frac{\omega ^n}{n!}=(1)^{n(n1)/2}\mathrm{\Omega }\overline{\mathrm{\Omega }},$$
$`(1)`$
where $`\omega `$ is the symplectic $`2`$-form corresponding to $`g.`$

An important property of a Calabi-Yau manifold is that it is Ricci-flat, i.e. the Ricci curvature vanishes. Conversely, on a simply connected open subset a Ricci-flat Kähler metric yields a Calabi-Yau structure. The most simple example of Calabi-Yau manifold is $`^n`$ equipped with its standard structures, that is

$$\omega _0=i\underset{j}{}dz_jd\overline{z}_j$$

and

$$\mathrm{\Omega }_0=dz_1\mathrm{}dz_n.$$

###### Definition 2

A submanifold $`L`$ of real dimension $`n`$ of $`X`$ is said to be special Lagrangian if $`\omega |_L\mathrm{Im}\mathrm{\Omega }|_L0.`$

The first condition says that $`L`$ is Lagrangian and the second one that $`L`$ is calibrated with respect to $`\mathrm{Re}\mathrm{\Omega }.`$ This implies that a special Lagrangian is always necessarily minimizing (cf \[HL\]).

We now describe a class of almost Calabi-Yau manifolds. Let $`P`$ be a complex polynomial. Then the set

$$Q:=\{(z_0,\mathrm{},z_n),P(z_0)=\underset{j=1}{\overset{n}{}}z_j^2\}$$

is a complex submanifold of $`^{n+1}`$ which is smooth except if $`P`$ admits double poles. For example, if $`P(z)=z^2,`$ $`Q`$ is made of two copies of $`^n`$ singularly connected at $`0.`$ $`Q`$ inherits from $`^{n+1}`$ a Kähler structure.

On $`Q`$ we define the following holomorphic $`(n,0)`$-form $`\mathrm{\Omega },`$ defined by

$$\mathrm{\Omega }d(P(z_0)+z_1^2+\mathrm{}+z_n^2)=dz_0\mathrm{}dz_n$$

We may give an explicit form on the set $`\{P(z_0)0\}:`$

$$\mathrm{\Omega }=\frac{1}{\dot{P}(z_0)}dz_1\mathrm{}dz_n$$

In the remainder of the section we shall see that in the particular case of the complex sphere, that is when $`P(z)=1z^2,`$ $`Q`$ is a Calabi-Yau manifold. This is a consequence of the following result due to Stenzel (cf \[St\]):

###### Proposition 1

On the complex sphere $`\{_{j=0}^nz_j^2=1\},`$ there exists a Ricci-flat metric, whose corresponding symplectic form is

$$\omega _{St}=i\overline{}u(r^2)=i\underset{j,k}{}\frac{^2}{z_j\overline{z}_k}u(r^2)dz_jd\overline{z}_k$$

where $`r^2=_{j=0}^nz_j\overline{z}_j`$ and $`u`$ is some smooth real function.

In dimension 2, we have an explicit form: $`u(r^2)=\sqrt{1+r^2},`$ and in higher dimension, $`u`$ is defined as the solution of some differential equation. As the complex sphere is simply connected, this result shows that it is Calabi-Yau.

It remains to show that the holomorphic $`n`$-form $`\mathrm{\Omega }`$ that we have introduced is the one of the Calabi-Yau structure, up to a multiplicative constant. So see this, let $`f`$ be the ratio of these two holomorphic forms. We observe that $`\mathrm{\Omega }`$ is invariant under the natural action of $`SO(n+1)`$ defined by $`z=x+iyAz=Ax+iAy,ASO(n+1).`$ But this is also the case for the Stenzel metric, and thus for the corresponding holomorphic $`n`$-form. This implies that $`f`$ is constant on the orbits of the action of $`SO(n+1).`$ But generic orbits have codimension $`1`$ in the complex sphere, so the image of $`f`$ in $``$ has at most real dimension 1 locally, which forces $`f`$ to be constant, as it is holomorphic. Thus the two forms are proportional.

## 2 A Lagrangian Ansatz

Let $`\gamma `$ be a curve into $`.`$ Then

$$\begin{array}{cccc}X:\hfill & I\times 𝕊^{n1}& & Q\\ & (s,x)& & (\gamma (s),\sqrt{P(\gamma (s))}x_1,\mathrm{},\sqrt{P(\gamma (s))}x_n),\end{array}$$

where $`(x_1,\mathrm{},x_n)=x,`$ is a an Lagrangian immersion for $`\omega _0.`$ In spite of the indetermination of the complex square root, the immersion is well defined because of the invariance of $`𝕊^{n1}`$ by $`xx.`$ However it becomes singular when $`\gamma `$ is a zero of $`P.`$

From now on we consider only the case $`P(\gamma )=1\gamma ^2,`$ that is the complex sphere equipped with the Stenzel form

Due to the larger symmetry of the complex sphere, we may slightly generalise the previous ansatz: we embed $`𝕊^{n1}`$ into $`𝕊^n`$ in the following way: $`𝕊^{n1}=𝕊^n\{x_0=0\}.`$ For some $`p𝕊^n,`$ we shall note $`A_p`$ any element of $`SO(n+1)`$ such that $`A_p(\{x_0=0\})=p^{}.`$ Let $`\gamma `$ be a curve into $`^{}.`$ Then we define

$$\begin{array}{cccc}X_{p,\gamma }:\hfill & I\times 𝕊^{n1}& & Q\\ & (s,x)& & \gamma (s)p+\sqrt{1\gamma (s)^2}A_px,\end{array}$$

and it is clear that the image does not depend on the choice of $`A_p.`$

###### Lemma 1

$`X`$ is also Lagrangian for the Stenzel metric.

Proof. Let $`(z_0,\mathrm{},z_n)`$ be coordinates on $`^n`$ such that $`p=(1,0,\mathrm{},0).`$ As long as $`\gamma `$ does not vanish we can use $`(z_1,\mathrm{},z_n)`$ as coordinates on a neighbourhood of the image of $`X`$ in $`Q.`$ In particular, we have

$$\frac{z_0}{z_j}=\frac{z_j}{z_0}.$$

We compute that

$$\frac{}{z_j}r^2=\overline{z_j}\frac{\overline{z}_0}{z_0}z_j$$

$$\frac{^2}{z_j\overline{z}_k}r^2=\delta _{jk}+\frac{z_j\overline{z}_k}{|z_0|^2}$$

So we deduce that

$$\omega _{St}=i\underset{j,k}{\overset{n}{}}a_{jk}dz_jd\overline{z}_k$$

with

$$a_{jk}=\left(\delta _{jk}+\frac{z_j\overline{z}_k}{|z_0|^2}\right)u^{}+2\mathrm{R}\mathrm{e}\left(\overline{z}_jz_k\frac{\overline{z}_0}{z_0}z_jz_k\right)u^{\prime \prime }$$

Hence we can decompose $`\omega _{St}`$ as

$$\omega _{St}=u^{}\omega _0+\omega _1,$$

where

$$\omega _1:=i\underset{j,k}{\overset{n}{}}\left(\frac{z_j\overline{z}_k}{|z_0|^2}u^{}+2\mathrm{R}\mathrm{e}\left(\overline{z}_jz_k\frac{\overline{z}_0}{z_0}z_jz_k\right)u^{\prime \prime }\right)dz_jd\overline{z}_k$$

As $`X`$ is Lagrangian for $`\omega _0,`$ it remains to show that it is also Lagrangian for $`\omega _1.`$ This will be a consequence of the fact that the two following $`2`$-forms vanish on $`X:`$

$$\underset{j,k}{\overset{n}{}}z_j\overline{z}_kdz_jd\overline{z}_k$$

$$\underset{j,k}{\overset{n}{}}z_jz_kdz_jd\overline{z}_k$$

We shall note $`T_\alpha =(X_\alpha ^1,\mathrm{},X_\alpha ^n),`$ $`1\alpha n1,`$ a basis of tangent vectors to $`𝕊^{n1}`$ at $`x=(x_1,\mathrm{},x_n).`$

This yields a basis of the tangent space at $`X(s,x)`$:

$$X_s=(\dot{\gamma },\frac{\dot{\gamma }\dot{P}(\gamma )}{2\sqrt{P(\gamma )}}x_1,\mathrm{},\frac{\dot{\gamma }\dot{P}(\gamma )}{2\sqrt{P(\gamma )}}x_n)$$

$$X_\alpha =(0,\sqrt{P(\gamma )}X_\alpha ^1,\mathrm{},\sqrt{P(\gamma )}X_\alpha ^n).$$

Now we have:

$$\left(\underset{j,k}{}\overline{z}_jz_kdz_jd\overline{z}_k\right)(X_s,X_\alpha )=\underset{j,k}{}|P(\gamma )|x_jx_k\left|\begin{array}{cc}\frac{\dot{\gamma }\dot{P}(\gamma )}{2\sqrt{P(\gamma )}}x_j& \sqrt{P(\gamma )}X_j^\alpha \\ & \\ \overline{\frac{\dot{\gamma }\dot{P}(\gamma )}{2\sqrt{P(\gamma )}}}x_k& \overline{\sqrt{P(\gamma )}}X_k^\alpha \end{array}\right|$$

$$=\underset{j,k}{}|P(\gamma )|x_jx_k\left(\frac{\dot{\gamma }\dot{P}(\gamma )}{2\sqrt{P(\gamma )}}\overline{\sqrt{P(\gamma )}}x_jX_\alpha ^k\overline{\frac{\dot{\gamma }\dot{P}(\gamma )}{2\sqrt{P(\gamma )}}}\sqrt{P(\gamma )}x_kX_\alpha ^j\right)$$

$$=|P(\gamma )|\frac{\dot{\gamma }\dot{P}(\gamma )}{2\sqrt{P(\gamma )}}\overline{\sqrt{P(\gamma )}}\underset{j}{}x_j^2\underset{k}{}x_kX_\alpha ^k|P(\gamma )|\overline{\frac{\dot{\gamma }\dot{P}(\gamma )}{2\sqrt{P(\gamma )}}}\sqrt{P(\gamma )}\underset{k}{}x_k^2\underset{j}{}x_jX_\alpha ^j$$

The latter vanish because tangent vectors to the real sphere at some point are orthogonal to this point.

$$\left(\underset{j,k}{}z_j\overline{z}_kdz_jd\overline{z}_k\right)(X_\alpha ,X_\beta )=\underset{j,k}{}|P(\gamma )|x_jx_k\left|\begin{array}{cc}\sqrt{P(\gamma )}X_j^\alpha & \sqrt{P(\gamma )}X_j^\beta \\ & \\ \overline{\sqrt{P(\gamma )}}X_k^\alpha & \overline{\sqrt{P(\gamma )}}X_k^\beta \end{array}\right|$$

$$=\underset{j,k}{}|P(\gamma )|x_jx_k\left(P(\gamma )X_j^\alpha X_k^\beta P(\gamma )X_k^\alpha X_j^\beta \right)=0$$

The computations for showing that also $`_{j,k}z_j\overline{z}_kdz_jdz_k`$ vanishes on $`X`$ are analogous.

## 3 Special Lagrangian submanifolds

Into the class of Lagrangian defined in the previous section, and in the two cases where we have a Calabi-Yau structure, we shall look for those calibrated by $`\mathrm{Re}\mathrm{\Omega },`$ that is those on which $`\mathrm{Im}\mathrm{\Omega }`$ vanishes.

### 3.1 Flat case of $`P(z)=z^2`$

The holomorphic $`(n,0)`$-form we use is the standard one

$$\mathrm{\Omega }_0=dz_1\mathrm{}dz_n$$

The equation is $`\mathrm{Im}(\dot{\gamma }\gamma ^{n1})=0,`$ which is easily integrated as $`\mathrm{Im}(\gamma ^n)=C,`$ where $`C`$ is some real constant. It turns that we have two types of solutions:

* If $`C=0,`$ $`\gamma `$ is the union of $`n`$ lines passing trough the origin, with angles $`\pi /n[2\pi ]`$ and the corresponding special Lagrangian are just Lagrangian $`n`$-spaces,
* If $`C0,`$ the curve is asymptotic to two of the lines described above. All the curves (and so the corresponding Lagrangian as well) are congruent. We recognize the Lagrangian catenoid which was first identified in \[HL\] and characterised in \[CU1\] as the only special Lagrangian of $`^n`$ which is foliated by $`(n1)`$-dimensional round spheres.

### 3.2 Case of the complex sphere $`P(z)=1z^2`$

Before starting the computations, we observe that we already know a special Lagrangian in the complex sphere: if a Calabi-Yau has a real structure which is in some sense compatible (cf \[Au\], pp. 62–64), then the set of real points, if non empty, is a special Lagrangian submanifold. In our case, real points constitute an real sphere embedded in $`Q,`$ so this will be a trivial solution to our problem.

We now compute the holomorphic form on $`X:`$

$$\mathrm{\Omega }(X_s,X_1,\mathrm{},X_{n1})=\frac{1}{\dot{P}(\gamma )}\left|\begin{array}{cccc}\frac{\dot{\gamma }\dot{P}(\gamma )}{2\sqrt{P(\gamma )}}x_1& \sqrt{P}X_1^1& \mathrm{}& \sqrt{P}X_{n1}^1\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \frac{\dot{\gamma }\dot{P}(\gamma )}{2\sqrt{P(\gamma )}}x_n& \sqrt{P}X_1^n& \mathrm{}& \sqrt{P}X_{n1}^n\end{array}\right|$$

$$=\frac{\dot{\gamma }}{2\sqrt{P(\gamma )}}\sqrt{P(\gamma )}^{n1}\left|\begin{array}{cccc}x_1& X_1^1& \mathrm{}& X_{n1}^1\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ x_n& X_1^n& \mathrm{}& X_{n1}^n\end{array}\right|$$

$$=\dot{\gamma }\sqrt{P(\gamma )}^{n2}C_{},$$

where $`C_{}`$ is some real constant.

Thus the differential equation for $`\gamma `$ is

$$\mathrm{Im}\left(\dot{\gamma }\sqrt{1\gamma ^2}^{n2}\right)=0.$$

This equation has two singular points $`\pm 1`$ and is regular elsewhere. Moreover, there is always a simple solution to this equation, the real segment $`[1,1].`$ This corresponds to the standard embedding of the real sphere. We shall analyse separately the even and odd cases.

#### 3.2.1 Even case

In this case, the equation is polynomial and it is easy to find a first integral: we integrate $`\mathrm{Im}\left(\dot{\gamma }(1\gamma ^2)^{n/21}\right)=0`$ and we get $`\mathrm{Im}(Q(\gamma ))=C,`$ where $`Q`$ is some polynomial such that $`Q^{}(z)=(1z^2)^{n/21}.`$ The degree of $`Q`$ is $`n1`$ and the the integral curves are algebraic curves of degree $`n1`$ as well.

In order to have a better description of the solutions, we perform an asymptotic analysis of the equation when $`|\gamma |\mathrm{}`$ and when $`\gamma \pm 1.`$

If $`|\gamma |>>1,`$ we have

$$\mathrm{Im}\left(\dot{\gamma }(1\gamma ^2)^{n/21}\right)\mathrm{Im}(\dot{\gamma }\gamma ^{n2}),$$

so asymptotically, the phase portrait of the equation looks like the one of the flat case (however not of the same dimension). In particular all integral curves are asymptotic to the half lines $`\{\mathrm{arg}(z)=k\pi /(n1)\}`$ , $`1k2n2.`$

Next we write $`\gamma =1+y`$ and we find that for $`y`$ close to $`0,`$

$$\mathrm{Im}\left(\dot{\gamma }(1\gamma ^2)^{n/21}\right)\mathrm{Im}\left(\dot{y}(2y)^{n/21}\right)$$

So close to the singular point $`1,`$ the integral lines look like the level sets of

$$\mathrm{Im}(1/n(2y)^{n/2})=C.$$

In particular, the level set of level $`0`$ is the union of $`n`$ branches asymptotic in $`1`$ to the half-lines $`\{\mathrm{arg}(z)=2k\pi /n+\pi \}`$ , $`1kn,`$ one of them being of course the real segment $`[1,1].`$

At the point $`1`$ the situation is exactly symmetric with respect to the reflexion $`z\overline{z}.`$

#### 3.2.2 Odd case

When $`n`$ is odd, there is also a first integral which takes the following form

$$\mathrm{Im}\left(\gamma \sqrt{1\gamma ^2}R(\gamma )+A\mathrm{arcsin}(\gamma )\right),$$

where $`R`$ is a polynomial of degree $`n3`$ and $`A`$ a real constant. This quantity is not well-defined on the whole complex plane, however the level set lines can be defined locally.

If $`|\gamma |>>1,`$ we have

$$\mathrm{Im}\left(\dot{\gamma }\sqrt{1\gamma ^2}^{n2}\right)=\mathrm{Im}\left(\dot{\gamma }\sqrt{1+(i\gamma )^2}^{n2}\right)\mathrm{Im}(\dot{\gamma }(i\gamma )^{n2})=\pm \mathrm{Im}(i\dot{\gamma }\gamma ^{n2})$$

Again, the phase portrait looks asymptotically as in the flat case, but this time up to a rotation of angle $`\pm \pi /2(n1).`$ In particular, all curves are asymptotic to the half lines $`\{\mathrm{arg}(z)=k\pi /(n1)+\pi /2(n1)\}`$ , $`1k2n2.`$

The asymptotic analysis next to the singular points $`\pm 1`$ is the same than in the even case.

#### 3.2.3 Conclusion

From these remarks, we can describe the general picture of the phase portrait :

* There is a singular curve made of the real segment $`[1,1]`$ and of $`2n2`$ branches, one half of them starting from $`1`$ and the other half from $`1,`$ making between them an angle of $`2\pi /n`$ and going to infinity where they are asymptotic to the half lines $`\{\mathrm{arg}(z)=k\pi /(n1)\}`$ , $`1k2n2`$ when $`n`$ is even and to $`\{\mathrm{arg}(z)=k\pi /(n1)+\pi /2(n1)\}`$ , $`1k2n2`$ when $`n`$ is odd.
* the other curves are smooth and have two ends asymptotic to two successive branches described above

For the corresponding special Lagrangian in $`Q,`$ we deduce the following:

* The standard embedding of the real sphere,
* $`2n2`$ special Lagrangian, touching the real sphere at the north pole $`(1,0,\mathrm{},0)`$ for one half of them and at the south pole $`(1,0\mathrm{},0)`$ for the other half.
* A one-parameter smooth families of special Lagrangian with two ends asymptotic to two members of the above family.

###### Remark 1

In the case of dimension $`2,`$ the integral lines are simply the horizontal lines.

Henri Anciaux,

IMPA, Estrada Dona Castorina, 110

22460–320 Rio de Janeiro, Brasil


text:
# Cellular automaton for bacterial towers

## 1 Introduction

Bacteria can either grow freely or adhere to surfaces. In the latter case one obtains structures ranging from compact colonies with self-affine surfaces, over diffusion-limited self-similar aggregates, to towering-pillar or mushroom shaped biofilms . Bacterial biofilms are commonly present in daily life; they typically flourish on substrates in aqueous environments, from sewage pipes to human teeth. In addition to performing physical measurements, it is necessary to introduce theoretical models for biofilm development. These models can be used in the prediction and/or understanding of the system’s evolution in experiment or simulation and they can also improve the understanding of intrinsic features. As an example, the efficiency of an antibiotic can depend on the spatial distribution and organization of bacteria.

Two different approaches have been developed over time: they use either microscopic or macroscopic rules. While macroscopic models can give more easily a global image of the biofilm, the microscopic ones describe better the behaviour at the level of the smallest biofilm unit - the “cell”. Ideally suited as elemental units of a cellular automaton approach, these cells contain typically a very small number of bacteria, or lots of nutrient particles, or simply consist of water-filled compartments. Biophysical properties are not imposed a priori but emerge as a result of local microscopic rules .

One of the best known theoretical biology models for population growth, which has since long been applied in various other fields, is the Eden Model . Various features have been added to this model, in order to describe different types of situations that appear in specific cases. Other approaches are diffusion limited aggregation (DLA) and, more generally, the cellular automata methods . In this article we shall present a method based on this last technique, which amounts to nutrient-limited aggregation in which chemotaxis is the driving mechanism and diffusion plays a minor role.

## 2 Model description

The general framework we employ is the cellular automata (CA) method. CA are discrete dynamical systems whose behavior is completely specified in terms of a local relation (rules). CA can be thought of as a stylized universe. Space is represented by a uniform grid, with each site or “cell” containing a few bits of data. Time advances in discrete steps, and the laws of the universe are expressed by a simple recipe - say, a small look-up table - through which at each step each cell computes its new state from that of its close neighbors. The system’s laws are local and uniform.

For our purpose the environment is a three-dimensional simple cubic lattice of size ($`l\times w\times h=81\times 81\times 27`$). Each lattice site is in one of a finite number of distinct states. Usually the number of states is small but, in principle, any finite CA model over a finite alphabet can be defined. The state of a cell (site) changes according to some rules. There are two distinct classes of automata: synchronous, which allows for simultaneous update of all cells’ states, and asynchronous. The updating order can be deterministic or random in the asynchronous class. The local rules of a certain CA are essentially determined by the considered neighborhood, consisting of the 6 nearest neighbors, 12 next-nearest and 8 next-next-nearest neighbors of a given site in the cubic lattice, 26 in total. We have used asynchronous updating.

We start with a homogeneous substrate layer of bacteria of size $`81\times 81\times 10`$ submerged in water. The water initially occupies a $`81\times 81\times 17`$ volume. We deposit at random large nutrient chunks on the substrate, conveniently approximated in shape by compact rectangular pillars of cells, of size $`27\times 27\times 10`$. The pillars are deposited with probability $`P`$ over a grid of 9 macroscopic square substrate plaquettes of area $`27\times 27`$, each containing many substrate cells. This process creates a strongly non-uniform nutrient supply partly covering the substrate. Bacterial cells may multiply provided they eat neighboring nutrient cells (nutrient is then replaced with bacteria). This process is characterized by a growth probability $`G`$ and represents binary division coupled to elementary chemotaxis, as the bacteria grow and spread preferentially in the direction of highest nutrient concentration. Clearly, the greater the number of neighboring cells occupied by nutrient, the greater is the probability that a bacterial cell will divide. For the time being, we disallow diffusion, and leave an examination of nutrient diffusion effects to a later section.

Technically, we choose at random a pair of neighboring cells and if the pair consists of nutrient and bacteria we replace the nutrient by bacteria with probability $`G`$ and then pass on to another randomly chosen pair. We take measurements at time steps of ($`l\times w\times h`$) updates, so that each cell has on average been visited once in every time step, and record bacterial population and nutrient configuration.

## 3 Results

### 3.1 Nutrient-limited growth without diffusion

We assume a probability $`P`$ for putting at random nutrient pillars over the coarse-grained substrate grid and a probability $`G`$ for bacterial cells to consume adjacent nutrient cells and subsequently divide, thereby occupying also the nutrient cell. As can be seen in the vertical section shown in Figure 1, the bacteria develop towers in the places where the nutrient was deposited at random. Note the developing roughness of the originally flat bacteria/nutrient interface. Within a tower the growth process is Eden-like (in three dimensions), since all perimeter cells carry nutrient. Eventually, however, the bacteria fill up the towers and perfectly flat walls and rooftops result.

Let us denote by $`N(t)`$ the bacterial population, in units of cells, in excess of the initially present substrate population, so that $`N(0)=0`$. Plots of $`N(t)`$, measured every time step, reveal a characteristic sigmoid shape, common for various biological growth processes approaching a finite carrying capacity . Since the deposition of nutrient is stochastic, the carrying capacity $`K`$ depends on the sample, and consequently a $`K`$-distribution arises, each $`K`$ being simply equivalent to the total mass of deposited nutrient pillars. The speed of saturation of the population towards $`K`$ is, of course, faster when $`G`$ is bigger and for each population curve we can discern three distinct regimes as illustrated in Figure 2, for which we took $`P=0.33`$.

i) The initial growth is not exponential - it would be exponential in free binary division - but is limited by the amount of nutrient in the first horizontal layer adjacent to and in contact with the bacterial substrate. The population increase is just proportional to the contact area $`A`$ and to $`G`$, leading to a linear initial growth. This stage can be characterized by a linear saturation time $`\tau _{lin,sat}`$, proportional to $`1/G`$. This is the average time needed to build a bacterial tower in the absence of interface fluctuations (cf. “mean-field” approximation), and is found by simple linear extrapolation of the initial growth. ii) The growth soon accelerates as the area of the bacteria/nutrient contact boundary increases by interface roughening and the appearance of nutrient “holes” (see Fig.1). The roughness, caused by height fluctuations, is limited by the finite horizontal width of the nutrient pillar(s), resulting in a saturating contact area, again implying a linear growth regime. iii) This regime in turn crosses over to an exponentially fast convergence of the growth towards $`K`$, which is the asymptotic behaviour for long times. Indeed, the bacteria reaching the ceiling of the nutrient pillar fall short of nutrient and eventually all nutrient cells are consumed. The sigmoid shape arises from the succession of these three regimes i)-iii). The dependence of the bacterial population on time can be understood most clearly from the time evolution of the area $`A(t)`$ of the bacteria/nutrient contact interface, shown in Figure 3.

Actually, we have

$$\frac{dN(t)}{dt}GA(t)$$
(1)
Asymptotically for long times, the decrease of the nutrient mass, $`M(t)`$, which equals the increase of the bacterial population $`N(t)`$, is simply proportional to the number of nutrient cells, since these become isolated, whence independent, and surrounded by neighboring bacterial cells. Thus we have

$$dMMGdt$$
(2)
It follows that the asymptotic behavior of $`N`$ for long times can be written as

$$N(t)K(1e^{(tt_0)/\tau }),$$
(3)
where the saturation time $`\tau `$ is proportional to $`1/G`$ and depends only weakly on the probability $`P`$ of depositing nutrient pillars, because these can to a first approximation (especially for small $`P`$) be considered as “non-interacting”. Numerically, $`1/\tau 0.85\pm 0.1`$ for $`G=0.8`$, implying a very short time scale for the actual exponential saturation in Fig.2 (right) compared to the linear saturation time. The time $`t_0`$ denotes the moment for which a zero excess population, $`N(t_0)0`$, would result by extrapolating backwards the asymptotic behavior. This time equals about $`t_0=16`$ for all curves with $`G=0.2`$ and about $`t_0=5`$ for all curves with $`G=0.8`$ (Fig.2). It is remarkable that the three characteristic times $`\tau _{lin,sat}`$, $`\tau `$ and $`t_0`$ are, within the accuracy of our computations, independent of the carrying capacity $`K`$.

### 3.2 Nutrient-limited growth with nutrient diffusion

We have performed simulations similar to the ones described above, but in which we have included stochastic nutrient diffusion. A probability $`I`$ is introduced for the interchange of adjacent cells, one occupied with nutrient and the other with water. This model is closer to physical, or biological, reality since the nutrient in general is soluble in water. In our simulations diffusion turns out to be necessary for generating the observed porous morphologies with fractal-like surface roughness and possible overhangs leading to mushroom-shaped towers, also seen experimentally . Figures 4 and 5 show relevant vertical and horizontal sections through the developing biofilm structure.

In the presence of nutrient diffusion, bacterial reproduction generates a two-step biofilm growth. Firstly, Eden-like growth driven by chemotaxis takes place in the nutrient pillars (as for zero diffusion). Secondly, DLA-type of growth roughens and ramifies the bacterial towers as the wandering nutrient recombines by aggregation. Clearly, the long time DLA process determines the biofilm morphotype. A remarkable property of the fractal character of the emerging structure becomes detectable in Fig.5 (horizontal slices through the biofilm). Sections through the top layers reveal a cascade or hierarchy of length scales. Water “lakes” in bacterial “land” develop into bacterial “islands” surrounded by water, over a height difference of just a few layers. This is strikingly similar to the structure that is inherent in the hierarchical population model, where it is imposed ad hoc. In contrast, here in the simulations it occurs spontaneously.

In Figure 6 we observe once more the basic sigmoid shape of the $`N(t)`$-curves describing the time evolution of the bacterial population, for 10 samples. Thus the shape of the population curves is robust with respect to the new perturbation; it is not much influenced by diffusion, at least not for short and medium times. The new feature brought about by diffusion is the existence of a long tail in $`N(t)`$, beginning after the formation of bacterial towers and ending when all detached nutrient has diffused “back” towards, and is eaten by, the bacterial colonies. This tail is still exponential, but with a decay time roughly 100 times longer than in the absence of diffusion. For example, for $`P=0.33`$ and $`G=0.8`$, we find $`1/\tau =0.009\pm 0.001`$. Interestingly, this saturation rate turns out to be independent of $`I`$ provided $`I>0`$, at least to the accuracy of our computations so far (averaging over samples, and for $`I=0.2,0.4,0.6,0.8`$ and 1). One may speculate that this could have been anticipated qualitatively, since weak diffusion leads to minor detachment of nutrient and swift recombination. On the other hand, strong diffusion entails more detachment but also greater mobility, whence also swift aggregation. The DLA process’ efficiency seems largely independent of the diffusion probability.

The exponential character of the saturation can be understood as follows. For long times, long after the bacterial towers are formed, the remaining water-dispersed nutrient particles can be considered as independent walkers - an assumption verified by keeping track of the asymptotic monomer character of nutrient matter in the simulations - and thus the asymptotics is again governed by $`dMMdt`$, implying exponential decay.

Just like in the previous simulations without diffusion, a carrying capacity distribution is found, which is controlled by the initial stochastic deposition of nutrient pillars. For a given number of nutrient pillars, the carrying capacity $`K`$ is fixed. To illustrate how this $`K`$-distribution emerges we present a histogram (Fig.7) giving the frequency of occurrence of $`N(10)`$, the population after ten time steps, for 1000 samples with $`P=0.33`$, $`G=0.8`$ and $`I=0.8`$. Note the fine structure of the main peaks. To optimize the resolution we have excluded from the viewing window the trivial peak at $`N=0`$ (no nutrient was put) and a small but not negligible eighth peak (outside of view, to the right). Note that, as time progresses, these peaks sharpen and converge to a simple set of ten $`K`$-values corresponding to a population filling all nutrient cells originally contained in the pillars. This sharpening is particular to our simple model, and will no longer occur when additional realistic ingredients such as antibiotics, detachment of bacteria and rinsing water flow is allowed, causing fluctuations of $`N(\mathrm{})`$.

## 4 Comparison with the hierarchical population model

The random hierarchical population model  generates bacterial towers through a more macroscopic construction algorithm. A typical local population $`N(x,y;t)`$ grown over a substrate in the $`xy`$-plane, for $`t=4`$ (after four generations) and for $`P=0.33`$ is shown in Fig.8.

A vertical (left) and a horizontal (right) section through the density $`N(x,y)`$ of Fig.8 are shown in Fig.9. Obviously, bacterial colonies of the type depicted in Fig.4 (right) and in Fig.5 are fairly well represented by Figs.8-9. The differences reside in the details of the self-similar properties intrinsic to the asymptotic long-time behavior of the hierarchical model and not to that of the simulations.

Typical population curves $`N(t)`$ produced in the hierarchical model are displayed in Fig.10, for $`P=0.33`$ and length rescaling parameter $`\lambda =3`$. After a short transient, the populations saturate exponentially rapidly to their carrying capacity. For long times, self-similar carrying capacity distributions result , which bear a resemblance to that of Fig.7, except in the details of the fractal geometry. We now turn to the important question of whether we can provide a microscopic basis for the rescaling factor $`\lambda `$.

In the exactly soluble random hierarchical population model the average excess population of bacteria $`\overline{N}(t)`$ obeys the simple time evolution in discrete time $`t=n`$ ,

$$\overline{N}(n)=(1\lambda ^n)\overline{K},$$
(4)
where $`\overline{K}`$ denotes the quenched average carrying capacity. This quantity is proportional to the birth probability $`B`$ (in the present context the probability of death, $`D`$, equals zero, because no antibiotic or other bactericidal agent is considered). Therefore, $`\overline{K}B`$. The exponential saturation is characterized by a rescaling factor $`\lambda `$, which is introduced ad hoc in that model to incorporate various possible conditions leading to a convergent population (limited nutrient, temperature boost or quench, diffusion without convection, bacterial motility, etc.). The asymptotic behavior of $`\overline{N}`$ is conspicuous by rewriting (4) in the form

$$\overline{N}(t)=(1e^{(\mathrm{log}\lambda )t})\overline{K}$$
(5)
In our CA simulations without diffusion, the towers fill up completely and compactly with bacteria (eventually without surface roughening), and the average carrying capacity is proportional to $`P`$, the nutrient pillar deposition probability. Precisely the same asymptotic result can be achieved in the hierarchical population model, provided we set $`B=P`$ in the first generation, and $`B=1`$ in all subsequent generations, while restricting the growth to the plaquettes where towers were grown in generation 1 (i.e., setting $`B=0`$ for substrate areas where no nutrient was put). Clearly, the correspondence between $`\lambda `$ in that model and the variables of our present simulation is

$$\mathrm{log}\lambda =1/\tau G$$
(6)
This convergence rate applies not only to the average but also to each curve of Fig.10 individually, since the process is self-averaging for long times . From (6) we learn that a constant microscopic (“cellular”) growth probability $`G`$ combined with a finite nutrient supply, fixing $`K`$, leads to an exponential saturation, justifying the use of and providing a microscopic interpretation for the constant rescaling factor $`\lambda `$ in the hierarchical population model.

Including nutrient diffusion in the simulations affects the asymptotic long-time behavior in that a non-zero nutrient/water interchange probability $`I`$ strongly increases the saturation time. Further, rough biofilm surfaces result. The hierarchical model is ideally suited to generate similar roughness, with marginal fractal geometry, in the spatial distribution of the population density $`N(x,y;t)`$ (Fig.8). Note that this density does not provide full information about the three-dimensional biofilm structure - in contrast with the simulations which give the full 3D bacterial “cell” population $`N(x,y,z;t)`$ and nutrient “cell” mass $`M(x,y,z;t)`$. The hierarchical-model population density yields a compactification in the form of a single-valued density or height profile defined over a 2D support, the substrate. Thus, Fig.8 shows a relevant population landscape from the hierarchical model, conform with the present 3D simulations including diffusion. Once again, the length rescaling factor $`\lambda `$ can be interpreted using the more microscopic parameters $`G`$ and $`I`$ of the CA simulations. Our simulations so far suggest that (6) remains valid for $`I0`$.

## 5 Conclusions and perspectives

We have proposed a microscopic cellular automaton describing bacterial tower development on a planar substrate in a strongly non-uniform nutrient environment. The “cells” contain bacteria, nutrient or water. A growth probability $`G`$ governs the stochastic chemotactic consumption of nutrient cells by adjacent bacterial cells, causing the former to be replaced by bacteria. A probability $`I`$ controls the interchange of neighboring nutrient and water cells, provoking nutrient diffusion. Nutrient is deposited only once, at the beginning, in the form of drops or pillars, with probability $`P`$. The strongly non-uniform nutrient supply causes bacterial colonies to grow into towers, and the absence of nutrient refuelling induces exponential saturation of the bacterial population towards a carrying capacity $`K`$ which depends on $`P`$. In view of the quenched random deposition a $`K`$-distribution is found, rather than a single final population.

Typical sigmoid curves are obtained for the time-dependence of the bacterial population. The asymptotic behavior for long times can be mapped onto similar saturation behavior displayed in the hierarchical population model. The phenomenological length rescaling factor in that model, necessary to render a convergent population, can be justified at a more microscopic level through its correspondence to $`G`$ (and $`I`$) found in the present work. This finding lends support to the applicability of that model as a phenomenological tool for generating fractal towering-pillar biofilm structure.

The simulations can be improved easily by adding new features while keeping maximal simplicity. We consider especially the possibility to add some antibiotic, deposited with probability $`Q`$, which is capable of killing bacterial cells (with extinction probability $`E`$), and to look for its bactericidal efficiency. In this process an adjacent antibiotic/bacteria pair of cells gets replaced by a pair water/dead-cell. Further, bacterial motility can be included through an additional water/bacteria cell interchange probability, allowing bacterial cells to detach. We have implemented these ideas in a simulation in which on each plaquette nutrient is deposited with $`P=0.33`$, antibiotics is deposited with $`Q=0.33`$ and nothing is put with probability $`1(P+Q)`$. In a first approximation only dead cells are allowed to detach. All diffusive interchange probabilities (water/nutrient, water/dead-cell, and water/antibiotic) are set to $`I=0.6`$. Further, $`G=0.8`$ and $`E=1`$. Fig.11 shows a horizontal section through the resulting biofilm structure after about 20 time steps. Interestingly, a bacterial tower has arisen in a nutrient-rich region (right), but only a relatively small hole resulted in the neighbouring antibiotic-rich area (middle). The mechanism which prevents massive extinction is the formation of a layer of dead cells mixed with water, separating the lower-lying bacteria from the bactericidal agent. This feature has no counterpart in the hierarchical population model, which is symmetric with respect to interchange of “birth” and “death”, but may be instrumental in biofilm resistance to antibiotics.

In a next step, “rinsing” can be added, entailing removal of detached bacteria, antibiotics and/or nutrient by water flow. Bacterial division in low-nutrient surroundings can also be allowed for, decoupling reproduction from chemotaxis, etc. We stress that the avenues explored here, as well as in the hierarchical model, provide a non-uniform (random) environment for the growth of towering pillar and mushroom structures. We believe that these morphotypes are difficult to understand - and indeed do not emerge easily - in uniform environments assumed when using standard Eden or DLA type aggregation rules. Incidentally, we checked that our approach reproduces pure DLA patterns in the limit of infinite nutrient diffusion (or turbulent convection), randomizing the nutrient mass over space.

Acknowledgments.
This research was supported by the Flemish Programme FWO-G.0222.02 “Physical and interdisciplinary applications of novel fractal structures”. We thank Katarzyna Sznajd-Weron and Jan Żebrowski for useful discussions.