text:
# Large extinctions in an evolutionary model: The role of innovation and keystone species

## Abstract

The causes of major and rapid transitions observed in biological macroevolution as well as in the evolution of social systems are a subject of much debate. Here we identify the proximate causes of crashes and recoveries that arise dynamically in a model system in which populations of (molecular) species co-evolve with their network of chemical interactions. Crashes are events that involve the rapid extinction of many species and recoveries the assimilation of new ones. These are analyzed and classified in terms of the structural properties of the network. We find that in the absence of large external perturbation, ‘innovation’ is a major cause of large extinctions and the prime cause of recoveries. Another major cause of crashes is the extinction of a ‘keystone species’. Different classes of causes produce crashes of different characteristic sizes.

Sanjay Jain∗,†,‡,¶ and Sandeep Krishna∗,¶

Centre for Theoretical Studies, Indian Institute of Science, Bangalore 560 012, India

Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA

Jawaharlal Nehru Centre for Advanced Scientific Research, Bangalore 560 064, India

Emails: jain@cts.iisc.ernet.in, sandeep@physics.iisc.ernet.in

Major transitions in biological and social systems are attributed primarily to some combination of external perturbation, selection, novelty, and a complex internal dynamics and structure of the system . While empirical studies attempt to identify proximate causes for individual events (which requires stupendous effort because of the complexity of the systems and the difficulty of obtaining data) most modeling efforts have tended to focus on the statistics of the events taken as a whole (see, e.g.,  for reviews of models in macroevolution and finance, respectively). Here we discuss an evolutionary model that permits a precise dissection of the (often multiple) causes of individual events. While the model is abstract, highly simplified and motivated by chemical evolution, the structures and processes that arise in it seem to have the flavour of phenomena in biological and social evolution. These include the appearance of a ‘core’ and ‘periphery’ in the network structure of the system, a shifting balance between cooperative and competitive processes as structures evolve, ‘core-shifts’, history dependence, ‘keystone species’, ‘innovations’ that are ‘core-transforming’ or ‘dormant’, and others. Their precise mathematical formulation and analysis of their role in major transitions in this simplified context may help in constructing and analyzing more realistic models.

The Model
The system is an idealized prebiotic pond containing a set of $`s`$ chemical species. A given species $`j`$ can be a catalyst for the production of another species $`i`$ with some small a priori probability $`p`$. Then the presence of $`j`$ in the pond causes the population of $`i`$ to increase according to the rate equation for catalyzed chemical reactions. The catalytic relationship is represented graphically by an arrow from node $`j`$ to node $`i`$ in a directed graph representing the chemical network. The graph is completely specified by its adjacency matrix $`C=(c_{ij})`$, $`i,j\{1,\mathrm{},s\}`$. $`c_{ij}`$ is unity if there is a link from node $`j`$ to node $`i`$, i.e., if species $`j`$ catalyzes the production of species $`i`$, and zero otherwise. Each species has a population $`y_i`$ in the pond. $`x_iy_i/_jy_j`$ is its relative population; by definition $`0x_i1`$, $`_{i=1}^sx_i=1`$. On a certain timescale $`T`$ the relative populations of the species reach a steady state denoted $`𝐗(X_1,\mathrm{},X_s)`$ that depends on the catalytic network. We imagine the pond to be subject to periodic external perturbations on a timescale greater than $`T`$ in the form of tides, storms or floods. Such a perturbation can flush out existing molecular species from the pond (the ones with the least $`X_i`$ being the most likely to be eliminated) and bring in new chemical species from the environment whose catalytic relationships with existing species in the pond are quite different from the eliminated ones. Then, after the perturbation, the populations in the pond will evolve to a new steady state, be subject to another perturbation, and so on.

The precise rules are as follows: Initialization: Start with a random graph with $`s`$ nodes and ‘catalytic probability’ $`p`$. That is, for every ordered pair $`(i,j)`$ with $`ij`$, $`c_{ij}`$ is unity with probability $`p`$ and zero with probability $`1p`$. $`c_{ii}0`$ for all $`i`$ (to forbid self-replicating species). Fig. 1a is an example. Assign each $`x_i`$ a random number between 0 and 1 and uniformly rescale all $`x_i`$ so that $`_{i=1}^sx_i=1`$. Dynamics: First, keeping $`C`$ fixed, evolve $`𝐱`$ from its initial condition according to

$$\dot{x}_i=\underset{j=1}{\overset{s}{}}c_{ij}x_jx_i\underset{k,j=1}{\overset{s}{}}c_{kj}x_j$$
(1)
for a time $`T`$ large enough to reach its attractor. Denote $`x_i(T)X_i`$. Eq. 1 follows from the rate equation $`\dot{y}_i=_jc_{ij}y_j\varphi y_i`$ for the populations, which in turn is an idealization of rate equations in a well stirred chemical reactor.***The rate equation $`\dot{y}_i=k(1+\nu y_j)n_An_B\varphi y_i`$ follows from the reaction scheme $`A+B\stackrel{j}{}i`$, where $`A`$ and $`B`$ are reactants with populations $`n_A`$ and $`n_B`$, $`j`$ and $`i`$ are catalyst and product with populations $`y_j`$ and $`y_i`$ respectively, and $`\varphi `$ is a death rate or dilution flux in the reactor ($`k`$ is the rate constant for the spontaneous reaction and $`\nu `$ is the catalytic efficiency). If $`n_A,n_B`$ are large and fixed and the spontaneous reaction is much slower then the catalyzed reaction, this equation reduces to $`\dot{y}_i=cy_j\varphi y_i`$, where $`c`$ is a constant. A generalization of the latter equation is $`\dot{y}_i=_{j=1}^sc_{ij}y_j\varphi y_i`$ for the case where species $`i`$ has multiple catalysts. Eq. 1 follows from this by taking the time derivative of $`x_i=y_i/_{j=1}^sy_j`$. In the present model we make the idealization that all catalytic strengths are equal. The second (quadratic) term in Eq. 1 is needed to preserve the normalization of the $`x_i`$ under time evolution. Note that it follows automatically from the nonlinear relationship between $`x_i`$ and $`y_i`$ when the time derivative of $`x_i`$ is taken. Find the set of nodes with the least $`X_i`$. Second, pick a node (denoted $`k`$) from this set at random and remove this node from the graph along with all its links. Add a new node (also denoted $`k`$) to the graph which is connected randomly to the existing nodes according to the same catalytic probability $`p`$. Mathematically, this means that for every $`ik`$, $`c_{ki}`$ and $`c_{ik}`$ are reassigned to unity with probability $`p`$ and zero with probability $`1p`$, irrespective of the value they had earlier, and $`c_{kk}=0`$. Set $`x_k=x_0`$ (a small constant), perturb all other $`x_i`$ about their existing value $`X_i`$ by a small amount, and uniformly rescale all $`x_i`$ to preserve the normalization $`_{i=1}^sx_i=1`$. This procedure provides a new graph and a new initial condition for $`𝐱`$. Now return to the first step of the dynamics and iterate the procedure several times.The attractor configuration $`𝐗`$ is determined in this paper by its algebraic properties discussed later, not by numerically integrating Eq. 1. Hence we are effectively taking $`T=\mathrm{}`$ which allows us to prove certain results analytically. Preliminary simulations suggest that when Eq. 1 is numerically integrated for a fixed, but large, $`T`$, qualitatively similar results are obtained, although one can expect that some of the sharp transitions are now spread out over several graph update time steps.

This model, introduced in ref. , was inspired by the work in refs. . The removal of the least populated species implements selection , and its replacement by another randomly connected species implements the introduction of novelty into the system. The network character is shaped by the repeated application of these two incremental external perturbations along with the internal population dynamics.

Three regimes of behaviour and transitions between them. The system exhibits three regimes or phases of behaviour. This is illustrated in Fig. 2 which shows the number of populated species in the attractor, $`s_1`$, (those species for which $`X_i>0`$) vs. time ($`n`$, the number of graph updates) for a run with $`s=100`$ and $`p=0.0025`$. In the ‘random phase’ $`s_1`$ stays low with small fluctuations. In the ‘growth phase’, $`s_1`$ typically rises exponentially with occasional drops. Finally, in the ‘organized phase’ $`s_1`$ stays close to $`s`$, the maximum value it can take. The random and growth phases were discussed in . As is evident from Fig. 2, the organized and the growth phases exhibit occasional discontinuous transitions or ‘crashes’ in which a number of species suddenly go extinct (their $`X_i`$ become zero in a single time step). At the end of a crash the system is in the random or growth phase. This is followed by a recovery in which the system moves again towards the organized phase. In  it is shown that crashes are primarily ‘core-shifts’, a specific kind of change in the structure of the graphs (discussed below) and recoveries are due to the growth of ‘autocatalytic sets’. The main purpose of this paper is to elucidate the mechanisms which cause core-shifts. Crashes and recoveries also occur in a related model with negative links and variable link strengths , where they are more difficult to study analytically.

Definitions and notation
Autocatalytic set (ACS): An ACS is a subgraph each of whose nodes has at least one incoming link from a node of the same subgraph. (By a subgraph we mean a subset of nodes together with all their mutual links.) Thus an ACS contains a catalyst for each of its members . In each of Figs. 1b-k, the subgraphs formed by the set of all black nodes or all the black and grey nodes are ACSs. Figs. 1a and 1l do not have an ACS. For any subgraph $`A`$, let $`\lambda _1(A)`$ be the largest eigenvalue of the adjacency matrix of $`A`$. We denote $`\lambda _1(C)\lambda _1`$. It can be shown  that if the graph does not have an ACS, then $`\lambda _1=0`$, and if it does, then $`\lambda _11`$. $`\lambda _1`$ therefore represents a topological property of the network.

Dominant ACS: It can be shown that if $`\lambda _11`$ the subgraph comprising the populated species ($`X_i>0`$) must be an ACS , which will be referred to as the dominant ACS. The dominant ACS is uniquely determined by the graph and does not depend upon the initial condition for $`𝐱`$ (except for special initial conditions forming a set of measure zero, which we ignore). The subgraph formed by the set of all black and grey nodes in each of Figs. 1b-k is the dominant ACS for that graph. In addition to its topological significance, $`\lambda _1`$ also has a dynamical interpretation as being the ‘population growth rate’ of the dominant ACS. This follows from the fact that the attractor configuration $`𝐗`$ is always an eigenvector of $`C`$ with eigenvalue $`\lambda _1`$, i.e., $`_jc_{ij}X_j=\lambda _1X_i`$ . Thus, when $`\varphi =0`$, substituting $`y_iX_i`$ in the population dynamics equation $`\dot{𝐲}=C𝐲`$, one gets $`\dot{𝐲}=\lambda _1𝐲`$.

Core and periphery of the dominant ACS: The core of the dominant ACS of a graph $`C`$ (sometimes also referred to as the ‘core of $`C`$’) is the maximal subgraph, $`Q`$, from each of whose nodes all nodes of the dominant ACS can be reached along some directed path. The rest of the dominant ACS is its periphery. The subgraph of all black nodes in Figs. 1b-k constitutes the core of the dominant ACS in the graph.§§§Sometimes the dominant ACS consists of two or more disjoint subgraphs, as in Fig. 1i. Then the above definition applies to each component separately. There exist other ACS structures for which this definition is not adequate, e.g., two disjoint 2-cycles pointing to a single downstream node. Such structures arise rarely and can be treated by a more general definition of core and periphery without altering the main conclusions presented here. Every periphery node has an incoming path that originates from the core, but no outgoing path that leads to the core. The periphery can contain loops within itself (e.g., the 2-cycle between nodes $`36`$ and $`74`$ in Fig. 1e). One can prove that $`\lambda _1(Q)=\lambda _1(C)\lambda _1`$.

A subgraph with more than one node containing a directed path from each of its nodes to each of its other nodes is said to be irreducible . An irreducible subgraph is always an ACS but the converse is not true. However, the core (of each component of a dominant ACS) is an irreducible subgraph. Since $`\lambda _1(A)1`$ for any irreducible subgraph $`A`$, it follows that the latter is a ‘self-sustaining’ structure in the sense that if no other links were present in the graph, the nodes of $`A`$ would still have nonzero $`X_i`$ by virtue of their mutual links. $`\lambda _1`$ measures the ‘strength’ of the core in two ways: one, its intrinsic population growth rate, and two, its multiplicity of internal pathways. To see the latter, compare the increasing and decreasing pattern of $`\lambda _1`$ between $`n=2854`$ and $`5042`$ in the inset of Fig. 2 with the sequence of Figs. 1b-f. When the core (of every disjoint component of the dominant ACS) has exactly one cycle (Figs. 1b,f-i,k), then $`\lambda _1=1`$, and vice versa. Such a core is fragile due to the absence of any redundancy in its internal pathways; the removal of any link from such a core will cause the ACS property (of that component) to disappear.

Crash: a crash is a graph update event, $`n`$, in which a significant number (say $`>s/2`$) of the species go extinct, i.e., $`\mathrm{\Delta }s_1(n)s_1(n)s_1(n1)<s/2`$.

Core overlap: given any two graphs $`C`$ and $`C^{}`$ whose nodes are labeled, the core overlap between them, denoted $`Ov(C,C^{})`$, is the number of common links in the cores of $`C`$ and $`C^{}`$. If any one of them does not have an ACS, $`Ov(C,C^{})0`$.

Core-shift: A core-shift is a graph update event in which $`Ov(C_{n1},C_n)=0`$, i.e, there is no overlap between the cores of the dominant ACS before and after the event ($`C_n`$ is the graph at time step $`n`$).

Keystone species. One can consider the impact of the hypothetical removal of any species $`i`$ from the graph, irrespective of whether $`i`$ is the least populated or not. For example one can ask for the core of the graph $`Ci`$ that would result if species $`i`$ (along with all its links) were removed from $`C`$. A species $`i`$ will be referred to as a keystone species if $`Ov(C,Ci)=0`$. Thus a keystone species is one whose removal modifies the organizational structure of the graph (as represented by its core) drastically. In the ecological literature a keystone species is regarded as one whose elimination from the ecosystem would cause a significant fraction of species in the ecosystem to go extinct . We will see that likewise in the present model the removal of keystone species causes large crashes. (Indices of keystoneness based on $`Ov(C,Ci)`$ or on the change in $`s_1`$ caused by the removal of a species can also be defined.) Note that if $`\lambda _1=1`$ and the ACS has a single connected component (e.g., figs. 1b,f-h,k), every core node is necessarily a keystone species since its removal would destroy the cycle that constitutes the core.

Innovations. In the present model the new species $`k`$ at each time step together with its set of new links may be regarded as a ‘novelty’ introduced into the system. In the new attractor the new species may go extinct, i.e., $`X_k`$ may be zero, or it may survive, i.e., $`X_k`$ is non-zero. Let us define an innovation to be a novelty in which the relative population of the new species in the new attractor just after the novelty occurs is nonzero. This definition has the feature that an innovation always involves new connections. It does not use any exogenously defined notion of fitness. The only performance criterion it requires is that the new links should enable the new node to survive until the next graph update. Even this minimal requirement has nontrivial consequences. For instance when the new species receives an incoming link from an existing dominant ACS, it typically has a nonzero population in the new attractor. Each recovery process, which occurs due to the expansion of the dominant ACS during the growth phase , is an accumulation of just such innovations.

We will be interested in a special class of innovations in which the novelty creates a new populated irreducible subgraph. Such innovations create or add to a ‘self-sustaining structure’ in the graph, in the sense discussed earlier. For instance the appearance of the first ACS at $`n=2854`$ (see Fig. 1b) is such an innovation. There, species 90, which was a singleton before the event, went extinct and was replaced by a new species 90 that had an incoming and outgoing link to node 26. The two formed an irreducible subgraph whose $`\lambda _1`$ value was 1. This innovation in fact triggered the self organization of the network around this ACS. By $`n=3880`$ (see Fig. 1c) the core had grown to 18 nodes as a result of several events in which the new species was an addition to the core. Every event in which the core is strengthened by the addition of a new node is also an innovation in which a new populated irreducible subgraph is formed. When such an innovation arises, denote the new irreducible subgraph which includes the new species as $`N`$ (or $`N_n`$ at time step $`n`$.) $`N`$ or $`N_n`$ will stand for the maximal irreducible subgraph of which the new species $`k`$ is a member. It follows from the Perron-Frobenius theorem  that such an innovation necessarily increases the $`\lambda _1`$ value of some substructure in the graph, i.e., $`\lambda _1(N)>\lambda _1(Nk)`$. $`N_n`$ becomes the new core of the graph if it is ‘stronger’ than the old core. More precisely, one can show that $`N_n`$ will become the new core of the graph, replacing the old core $`Q_{n1}`$, whenever any one of the following two conditions hold:
(a) $`\lambda _1(N_n)>\lambda _1(Q_n^{})`$ or,
(b) $`\lambda _1(N_n)=\lambda _1(Q_n^{})`$ and $`N_n`$ is ‘downstream’ of $`Q_n^{}`$.We use the notation $`C_n^{}C_{n1}k`$ for the graph of $`s1`$ nodes just before the novelty at time step $`n`$ is brought in (and just after the least populated species $`k`$ is removed from $`C_{n1}`$). $`Q_n^{}`$ stands for the core of $`C_n^{}`$. A subgraph $`A`$ is ‘downstream’ of another subgraph $`B`$ if there exists a directed path from some node of $`B`$ to a node of $`A`$ but none from any node of $`A`$ to a node of $`B`$.
Such an innovation may be referred to as a core-transforming innovation. If $`Q_{n1}N_n`$, such an innovation enlarges the existing core. However, if $`Q_{n1}`$ and $`N_n`$ are disjoint, we get a core-shift.

In the present model, the generation of novelty does not depend upon existing structure since the links of a new node are chosen randomly from a fixed probability distribution: novelty is ‘noise’ in this model. (Variants of the model that depart from this are easily constructed.) However, whether the novelty constitutes an innovation is ‘context dependent’ (i.e., dependent upon the structure of the existing network). Also, the short and long term impact of an innovation depends upon the (historical evolution of the) ‘context’, as will be seen below.

Classification of core-shifts
In a set of runs with $`s=100,p=0.0025`$ totaling 1.55 million iterations we observed 701 crashes. 612 of these were core-shifts . Fig. 3 differentiates between the 612 core-shifts we observed. They fall into three categories: (i) complete crashes ($`136`$ events), (ii) takeovers by core-transforming innovations ($`241`$ events), and (iii) takeovers by dormant innovations ($`235`$ events).

Complete crashes. A complete crash is an event in which an ACS exists before but not after the graph update. Therefore for a complete crash at time step $`n`$, $`\lambda _1(C_{n1})>0`$ and $`\lambda _1(C_n)=0`$. These events take the system back to the random phase. For example at $`n=8232`$, in Fig. 1k, node $`54`$ is one of the least populated species and is hit at $`n=8233`$. It is replaced by a new species that has a single outgoing link to node $`50`$ and no incoming link, resulting in the complete disruption of the ACS. It is evident that complete crashes must always be caused by the elimination of a keystone species. Furthermore, Fig. 3 shows that $`\lambda _1(C_{n1})=1`$ for every complete crash observed in the runs. Hence the core of the ACS is a single cycle when the event occurs and the species removed is a member of that cycle.

Takeovers by core-transforming innovations. An example of a takeover by a core-transforming innovation is given in Figs. 1g,h. At $`n=6061`$ the core was a single loop comprising nodes $`36`$ and $`74`$. Node $`60`$ was replaced by a new species at $`n=6062`$. The new node $`60`$ created an innovation at $`n=6062`$, with $`N_{6062}`$ being the cycle comprising nodes $`60,21,41,19`$ and $`73`$, downstream from the old core. The graph at $`n=6062`$ has one cycle feeding into a second cycle that is downstream from it. (This is an example of condition (b) for a core-transforming innovation.) The population attractor of such a graph has the property that the $`X_i`$ of the nodes in the upstream cycle are all zero and only the nodes of the second cycle (as well as nodes further downstream from it) have nonzero $`X_i`$. Thus when the above innovation arises, the new cycle becomes the new core and all nodes that are not downstream from it get $`X_i=0`$, resulting in a large drop in $`s_1`$ from 89 to 32. For all such events in Fig. 3, $`\lambda _1(Q_n^{})=\lambda _1(C_{n1})`$ since $`k`$ happened not to be a core node of $`C_{n1}`$. Thus these core-shifts satisfy $`\lambda _1(C_n)=\lambda _1(N_n)\lambda _1(Q_n^{})=\lambda _1(C_{n1})1`$ in Fig. 3.

Takeovers by dormant innovations. Figs. 1e,f show an example of a takeover by a dormant innovation. At $`n=5041`$, the core has $`\lambda _1=1.24`$ and there is a cycle comprising nodes $`36`$ and $`74`$ in its periphery. Node $`85`$ is hit which results in a cycle ($`26`$ and $`90`$) feeding into another cycle ($`36`$ and $`74`$) at $`n=5042`$. Thus at $`n=5042`$, for the same reason as in the previous paragraph, $`36`$ and $`74`$ form the new core with only one other downstream node, $`11`$, being populated. All other nodes become depopulated resulting in a drop in $`s_1`$ by $`97`$. Such a core-shift is the result of an innovation that arose earlier (the cycle between $`36`$ and $`74`$ arose at $`n=4695`$) but lay dormant downstream of the existing core until one of the keystone species of the latter (node $`85`$) was hit and made it weak (i.e., reduced its $`\lambda _1`$ to a value less than or equal to the $`\lambda _1`$ of the downstream innovation). A dormant innovation can takeover as the new core only following a keystone extinction which weakens the old core. In such an event the new core necessarily has a lower (but nonzero) $`\lambda _1`$ than the old core, i.e., $`\lambda _1(C_{n1})>\lambda _1(C_n)1`$. Note that at $`n=5041`$ if the downstream cycle between $`36`$ and $`74`$ were absent, $`85`$ would not be a keystone species by our definition, since its removal would still leave part of the core intact (nodes $`26`$ and $`90`$). $`85`$ becomes keystone, and the core of which it is a part becomes fragile and susceptible to a core-shift, because a self-sustaining innovation has occurred in the distant periphery.

If we consider the set of all drops in $`s_1`$, large or small, where an ACS exists before the event (there are 126454 such events in the above mentioned runs) we find that the number of complete crashes remains the same, 136 (the mean size of the drop in such events is $`\overline{|\mathrm{\Delta }s_1|}=98.2`$ with a standard deviation $`\sigma =1.2`$), while core-shifts caused by dormant innovations go up to 359 (with $`\overline{|\mathrm{\Delta }s_1|}=62.2`$, $`\sigma =25.9`$) and those due to core-transforming innovations to 524 ($`\overline{|\mathrm{\Delta }s_1|}=48.2`$, $`\sigma =25.6`$). The rest of the events consist of 9851 ‘partial core-shifts’ (in which the core changes, but $`Ov(C_{n1},C_n)0`$; in this category $`\overline{|\mathrm{\Delta }s_1|}=2.18`$, $`\sigma =7.42`$), and 115584 events where there was no change in the core but the periphery is affected ($`\overline{|\mathrm{\Delta }s_1|}=1.05`$, $`\sigma =0.99`$). Thus different classes of proximate causes of drops arise dynamically with different frequencies and typically produce events in different size ranges. The ranges, however, overlap and some distributions have fat tails (e.g., of the 701 crashes with $`|\mathrm{\Delta }s_1|>50`$, there are 79 and 10 events, respectively, in the last two categories). Detailed distributions and their dependence on $`s,p`$ are open questions (but see ).

Discussion
The present model exhibits mechanisms by which innovations can play a major role in crashes and recoveries in a complex system.Analogues of innovations and core-shifts appear to be playing an important role in another related but quite different model  where rapid transitions are observed. It provides a mathematical example of ‘creative destruction’  at work in causing large upheavals. It distinguishes two processes involving innovations, both having analogues in the real world. One is exemplified by the appearance of the automobile which made the horse drawn carriage and its ancillary industries obsolete. This is like the example of the core-transforming innovation shown in Figs. 1g,h where the graph update produced a self-sustaining structure that was more vibrant than the existing core within the context of the present organization. This structure became the new core, rendering many nodes drawing sustenance from the old core dysfunctional. The subsequent development of other industries dependent on the automobile mirrors the growth of the ACS around the new core. The second process is exemplified by the emergence of the body plans of several phyla which are dominant today. It is believed that while these body plans originated in the Cambrian era more than 520 million years ago , the organisms with these body plans played no major role till about 250 million years ago. They started flourishing only when the Permian extinction depleted the other species that were dominant till that time . This is similar to the events shown in Figs. 1e,f where an earlier innovation had lain dormant for a while without disturbing the existing core, but when the latter became sufficiently weak, took over as the new core and flourished.

Recently there has been substantial progress in graph theoretic analyses of complex systems and in particular ‘small-world’  and ‘scale-free’  properties have been found for several real networks (for reviews and references see ). It may be interesting to study whether certain classes of real networks also have some kind of a ‘core-periphery’ structure, in view of its possible dynamical significance.

Acknowledgements. S.J. acknowledges the Associateship of Abdus Salam International Centre for Theoretical Physics, Trieste as well as the hospitality of the Max Planck Institute for Mathematics in the Sciences, Leipzig, and Jesus College, Cambridge. This work was supported in part by a grant from the Department of Science and Technology, Government of India.

Figure legends

Figure 1. The structure of the evolving graph at various time instants for a run with $`s=100,p=0.0025`$. Node numbers $`i`$ from $`1`$ to $`100`$ are shown in the circles representing the nodes. Nodes with zero relative population in the attractor configuration for the graph ($`X_i=0`$) are shown as white circles; the rest ($`s_1`$ in number) have nonzero $`X_i`$. In the graphs where an autocatalytic set (ACS) exists, black circles correspond to nodes in the ‘core’ of the ACS, and grey to the ‘periphery’, defined in the text. (a) $`n=1`$, the initial random graph, (b) $`n=2854`$, where the first ACS, a 2-cycle between nodes $`26`$ and $`90`$, appeared, (c) $`n=3880`$, the beginning of the organized phase when the ACS first spanned the entire graph, (d) $`n=4448`$, when the core reached a peak in the number of loops it contained, (e) $`n=5041`$, just before a ‘core-shift’, (f) $`n=5042`$, just after the core-shift caused by a ‘keystone’ extinction in the presence of a ‘dormant innovation’, (g) $`n=6061`$, just before another core-shift, (h) $`n=6062`$, just after the core-shift caused by a ‘core-transforming innovation’, (i) $`n=6070`$, when the old core stages a come back as a disconnected component after node $`32`$ becomes a singleton, (j) $`n=6212`$, when the new core strengthens itself and depopulates the recently resurgent old core, (k) $`n=8232`$, just before the first ‘complete crash’, (l) $`n=10000`$, between the first ‘complete crash’ and the subsequent ‘recovery’.

Figure 2. The number of populated species, $`s_1`$ (continuous line), and the largest eigenvalue of $`C`$ (whose significance is discussed later in the text), $`\lambda _1`$ (dotted line), versus time, $`n`$ for a run with $`s=100`$ and $`p=0.0025`$. The $`\lambda _1`$ values shown are $`100`$ times the actual $`\lambda _1`$ value. The first $`10000`$ time steps are enlarged in the inset. The run is the same for which the graph snapshots are shown in Fig. 1. The impact of the events described in Fig. 1 is clearly visible in this curve. E.g., at $`n=2854`$ $`\lambda _1`$ jumps from zero to one and $`s_1`$ exhibits the first sustained upward movement, at $`n=3880`$ $`s_1`$ hits its maximum value, $`100`$, and then fluctuates mainly between $`99`$ and $`100`$, and at $`n=4448`$ $`\lambda _1`$ reaches a local maximum. $`s_1`$ drops from $`100`$ to $`3`$ as a result of the ‘core shift’ at $`n=5042`$, and from $`89`$ to $`32`$ in the core-shift at 6062. At $`n=6070`$ a large recovery event is seen as the old core and the still intact part of its periphery get repopulated, only to be trounced again at $`n=6212`$ when the new core strengthens itself to a $`\lambda _1`$ value greater than $`1`$. At $`n=8233`$ $`s_1`$ crashes from $`100`$ to $`2`$ when the ACS is completely destroyed and $`\lambda _1`$ drops from one to zero.

Figure 3. Frequency, $`f`$, of the $`612`$ core-shifts observed in a set of runs with $`s=100`$ and $`p=0.0025`$ vs. the $`\lambda _1`$ values before, $`\lambda _1(C_{n1})`$, and after, $`\lambda _1(C_n)`$, the core-shift. Complete crashes (black; $`\lambda _1(C_{n1})=1`$, $`\lambda _1(C_n)=0`$), takeovers by core-transforming innovations (blue; $`\lambda _1(C_n)\lambda _1(C_{n1})1`$) and takeovers by dormant innovations (red; $`\lambda _1(C_{n1})>\lambda _1(C_n)1`$) are distinguished. Numbers alongside vertical lines represent the corresponding $`f`$ value.

Figure 3.


format: text

text:
# High Speed photometry of faint Cataclysmic Variables: I. V359 Cen, XZ Eri, HY Lup, V351 Pup, V630 Sgr, YY Tel, CQ Vel, CE-315

## 1 Introduction

Cataclysmic Variable stars (CVs) are a class of interacting close binaries which attract great attention. They enable a quantitative study of the many processes involved in mass transfer through accretions discs or through magnetically channeled accretion. A list of these physical processes, which occur also in many other astrophysical environments, is given in Warner (2000); a review of CVs in general is given in Warner (1995a).

Fewer than 40 CVs have been studied in any great detail (e.g., with doppler tomography, eclipse deconvolution). These, and the phenomena seen in other CVs, show that no two CVs are exactly alike. The different combinations of masses of primaries and secondaries, temperatures, magnetic field strengths, orbital inclinations and mass transfer rates ($`\dot{M}`$) suffice to produce great variety of behaviour. Much of this variety is already understood in terms of the component physics accompanying the gross differences. The effects of more subtle differences – e.g., original chemical composition of the secondary, effects on the secondary of past nova explosions, compositional and environmental (for instance, irradiation) dependence of disc viscosity – have still to be teased out of the rich observational phenomenology.

Much of the progress in the past two or three decades has come from the discovery of ever more variety of behaviour. There are often conspicuous effects in the longer term light curves; examples are the ER UMa stars (dwarf novae with outbursts every few days and superoutbursts every few weeks), ‘echo’ outbursts after a superoutburst, continuous small outbursts in the discs of nova-like variables. New phenomena continue to be found also in the ‘close-up’ details; examples are short period quasi-periodic oscillations (from $``$2 s in the X-ray flux of SS Cyg in outburst to $``$1000 s modulations in many high $`\dot{M}`$ discs), spiral waves in discs, tilted discs, desynchronisation of polars by nova explosions.

Although some of the new phenomena have accompanied the discovery of new CVs (e.g., the $`10^8`$ G magnetic fields in AR UMa and V884 Her), most have been discovered in long known CVs (e.g., TeV gamma rays from AE Aqr, ZZ Cet pulsations in GW Lib). It is clear, therefore, that while the great majority of CVs remain under-observed new phenomena may remain undiscovered. At a less exotic level, the determination of a significant number of new orbital periods, and in particular the discovery of eclipsing systems, provide objects for more detailed structural studies. With the advent of 8-m class telescopes such studies can be made even to 20th magnitude. The ongoing survey, of which this paper presents the first general results, was initiated with this in mind.

The need for such a survey, particularly in the southern hemisphere and capable of succeeding in the crowded fields in the general direction of the Galactic centre where the majority of nova remnants reside, was recognised in the early 1990s and as a result the CV group at the University of Cape Town (UCT) used several years’ equipment funds (in advance) for the purchase of a state-of-the-art CCD detector. The photometer and its software were designed and constructed by Dr. D. O’Donoghue, who was a member of the Astronomy Department at UCT at that time. The UCT CCD photometer, completed in 1993, has been extensively used on telescopes at the Sutherland site of the South African Astronomical Observatory. An outline of its functional properties can be found in O’Donoghue (1995) and Koen & O’Donoghue (1995). In the CV group at UCT it was first largely used for observations of well known CVs at higher time resolution or improved signal-to-noise. Its use for a general survey of faint CVs started in 1997, with the almost immediate discovery of the first ZZ Cet (i.e., non-radial) pulsations in the white dwarf primary of the SU UMa type dwarf nova GW Lib (Warner & van Zyl 1998) and independent discovery of the intermediate polar nature of the nova remnant HZ Pup (see Abbott & Shafter 1997 for detailed description). Follow-up observations of these led to neglect of the survey for some time, but we are now actively observing faint CVs again. We present here initial results.

The observations are given in Section 2 and some brief conclusions in Section 3.

## 2 Observations

All observations have been made with the UCT CCD photometer attached either to the 40-in or 74-in reflectors at Sutherland. In order to capture as many photons as possible ‘white light’ was used, with the result that extinction corrections and magnitude calibrations are only approximate.

The observing log for the stars included in this paper is given in Table 1.

### 2.1 V359 Centauri

V359 Cen was discovered bright in April 1929111Note that Duerbeck (1987) states 1939, and this has been propagated in later literature. and entered the lists as a probable nova, unconfirmed because no spectra were obtained. There are still no spectra available at minimum light, but in recent years further outbursts to a maximum of V $`13.8`$ have been observed at intervals $`1`$ y, showing the star to be of SU UMa type. The quiescent magnitude given in the Downes, Webbink & Shara (1997) catalogue is 21.0 from J plates, and Munari & Zwitter (1998) gave V $`20.5`$ in 1996, but we find V $``$ 18.7, which is a great deal easier to observe at moderate time resolution. Our May 2001 observation was obtained three weeks after a superoutburst, but the star had returned to its quiescent brightness.

Our light curves are displayed in Fig. 1 and show a recurrent hump with an amplitude $`0.40.5`$ mag. Our observations are not distributed adequately to remove aliases; the best determined period from the February 2001 run is 112 min (0.0779 days). The variations of hump profile and their occasional almost triangular shape resemble more what is seen in superhumps than orbital modulations, but it is rare to observe superhumps at quiescence – the only well documented examples are EG Cnc which sustained ‘late’ superhumps for at least 90 d after the 1996 superoutburst (at a time when ‘echo’ normal outbursts were occuring) (Patterson et al. 1998) and AL Com with a persistent superhump as well as an orbital hump at quiescence (Patterson et al. 1996). If what we are seeing is an orbital hump then the inclination of V359 Cen must be $`65^{}`$.

### 2.2 XZ Eridanus

XZ Eri has been known as a probable dwarf nova for nearly 70 years (Shapley & Hughes 1934). In recent years only three outbursts have been detected: Mar/Apr 1995, Jan 1998, Feb 1999. The first of these reached V $`14.1`$, the others reached V $`14.6`$. CCD photometry at quiescence gives V = 18.7 (Howell et al. 1991). The rarity and range of the outbursts suggests that they are superoutbursts. Any normal outbursts may have gone unnoticed. There has been no high speed photometry during an outburst.

A quiescence spectrum, obtained by Szkody & Howell (1992) shows the Balmer emission and continuum characteristic of a dwarf nova. There is a hint of doubling in the strongest Balmer lines, indicative of a high inclination system. Nearly 4 h of CCD photometry through a $`V`$ filter at quiescence, with 220 s integration on a 1.8-m telescope, showed a scatter of $`1.0`$ mag and, to quote the authors, no indications of periodic modulation (Howell et al. 1991).

Our observations (Fig. 2) show XZ Eri to be, in brief, a Z Cha lookalike, with the succession of immersions and emersions of the white dwarf and bright spot all quite clear. At V $`19`$ we cannot achieve sufficient time resolution to separate these components as clearly as in Z Cha (e.g., Wood et al. 1986). Using our best estimate for the time of mid-eclipse of the white dwarf we derive the following ephemeris:

$$\mathrm{HJD}_{\mathrm{min}}=2451905.4419+0.^\mathrm{d}0612\mathrm{E}.$$
(1)
The short $`P_{orb}`$ (1.47 h = 88.1 min) implies that XZ Eri will have superhumps when observed in outburst.

### 2.3 HY Lupi

HY Lup was Nova Lupi 1993, which reached V $`8`$ and was moderately fast ($`t_215`$ d). It is currently at V $`18.9`$ and has a resolved ejecta shell (Downes & Duerbeck 2000). No pre-outburst candidate brighter than V $`17`$ is detectable (McNaught & Garradd 1993) and it is probable that the remnant has not yet reached its quiescent brightness. No high speed photometry has been published for this nova remnant.

Our photometry (Fig. 3) shows little evidence for rapid flickering and only slow variations of small amplitude. There is no suggestion of any orbital modulation which would encourage further photometry. The $`P_{orb}`$ for HY Lup will need to be obtained spectroscopically.

### 2.4 V351 Puppis

Nova Puppis 1991 (V351 Puppis) was a moderately fast nova ($`t_2=16`$ d), probably a neon nova (Pachoulakis & Saizar 1995). The pre-eruption magnitude was variable over the range $`21>22`$ (McNaught 1992). The eruption light curve through 1996 is illustrated in Saizar et al. (1996); in 1998 the remnant was at V $`19.6`$ (Downes & Duerbeck 2000).

There are one or two additional aspects of this nova that are of interest. Sixteen months after optical maximum, hard X-rays were detected by ROSAT which, from the absence of a soft X-ray component, were deduced to be caused either by shocks in the ejected shell or by magnetically controlled accretion onto the white dwarf (Orio et al. 1996). Variability (time scale $`<`$ days) observed in the X-ray flux seemed to favour the latter interpretation. However, a later analysis explains the absence of soft X-rays as due to internal and external absorption, allowing the hard X-rays to be the product of surface nuclear burning, as in most novae (Vanlandingham et al. 2001). Multiwavelength observations of the nebular phase of the eruption made in 1992-1993 revealed a strong red excess attributed to the secondary star, implying that this component is either a giant or an irradiated dwarf (Saizar et al. 1996). In the former case an orbital period $`P_{orb}`$ days would be expected, in the latter $`P_{orb}4`$ h would be required in order for the secondary to be close enough to the hot primary to be strongly irradiated. The relatively small mass of the ejected shell of N Puppis 1991 implies a high mass for the white dwarf – near the Chandrasekhar limit (Vanlandingham et al. 2001).

We observed V351 Pup in December 2000 and found it to be modulated in brightness with a range of up to $``$1.2 mag and a period of 2.837 h. The light curves are given in Fig. 4 and show variable amplitude, low amplitude flickering, rounded maxima and a V-shaped minimum which may be filled in to form a flattened region. The ephemeris from these observations is

$$\mathrm{HJD}_{\mathrm{max}}=2451905.56746+0.^\mathrm{d}1182\mathrm{E}.$$
(2)
The additional observations we made in February 2001 are too far removed in time to enable us to improve on the precision of the period.

Our light curves for V351 Pup have a striking similarity to those of V1500 Cyg in the late 1980s (see especially Fig. 1 of Kaluzny & Chlebowski 1988). V1500 Cyg was Nova Cygni 1975 and in the late 1980s it was about 15.3 mag below its maximum brightness. Nova Pup 1991 was not discovered until after maximum, but there are reasons to believe it did not exceed V $`5.5`$ (Saizar et al. 1996) and it therefore at present is $``$14 mag below maximum. V1500 Cyg is a desynchronised polar (Stockman, Schmidt & Lamb 1988; Kaluzny & Chlebowski 1988) with an orbital period of 3.35 h, which suggests the possibility that V351 Pup may also be magnetic. With this in mind we asked Gary Schmidt to check for circular polarization and he reports “A sequence of spectropolarimetric observations was obtained by G. and P. Schmidt on 2 January 2001 using the 74-in reflector at Mt. Stromlo and covering the range 4220 – 7300 Å at $``$10 Å resolution. The summed spectrum, which spans a total of 3.7 h, reveals a weak continuum (V $`19.4`$) that is approximately flat in $`f_\lambda `$, plus broad ($`\pm 1700`$ km s-1) nebular emission lines of H, He II and O III. The polarimetric results fail to show significant circular polarization anywhere in the observational sequence, and the coadded sum is unpolarized to a 3$`\sigma `$ upper limit of 1%”.

Before dismissing the possibility that V351 Pup is a magnetic nova we note that (a) the level of circular polarization in V1500 Cyg is $`1.5`$% (Stockman et al. 1988), which, compared with the 10–30% in polars, shows the presence of much unpolarized light, probably that from irradiation of the secondary by the still hot primary; (b) V351 Pup is 10 y from its eruption and less far in its recovery than was V1500 Cyg in the late 1980s – and should therefore have even more unpolarised light; (c) GQ Mus (Nova Mus 1983) also has a polar-like orbital light curve and shows strong spectral signatures of magnetic accretion, but no circular polarization (Diaz & Steiner 1994); (d) similar remarks apply to V2214 Oph (Nova Oph 1988: Baptista et al. 1993).

### 2.5 V630 Saggitarii

V630 Sgr was the last of four novae detected in Saggitarius in 1936. It was discovered in October of that year, at a time when Harvard objective prism plates were being exposed, whereby it attained the status and alias HD 321353. The observed maximum was V $`4.5`$, but it probably reached V $`1.6`$ (Downes et al. 1997), and at minimum it is V $`19`$, giving it the unusually large eruption amplitude of $`17.5`$ mag, even for a $`t_2=4`$ d nova, suggesting a high inclination system (Warner 1986: this conclusion arises from the apparent intrinsic faintness of the nearly edge-on accretion disc which such an object has in quiescence).

There have been no modern photometric or spectroscopic observations of V630 Sgr, which is a faint object in a very crowded field. It is listed by Diaz & Steiner (1991), on the basis of its similarity to V1500 Cyg in being a very fast nova of large amplitude, as a candidate magnetic system.

Our light curves are shown in Fig. 5 and illustrate that V630 Sgr is indeed a moderately high inclination system, with eclipses 0.4 – 0.6 mag deep. The ephemeris derived from our light curves, measured from the fundamental Fourier component. The time of the first minimum of our first light curve of this set is given in the ephemeris, which is

$$\mathrm{HJD}_{\mathrm{min}}=2451781.2480+0.^\mathrm{d}1180\mathrm{E}.$$
(3)
The orbital period of 2.831 h is just at the top of the orbital period gap.

It is evident from Fig. 5 that there is a hump of complex and variable profile that has a period comparable to $`P_{orb}`$. The Fourier amplitude spectrum of the complete set of light curves is shown in Fig. 6. The main groups of power are associated with the fundamental and strong harmonics of the narrow eclipses. The window pattern of the data set can be judged from the first and second harmonics. Superimposed, but slightly offset from the $`P_{orb}`$ fundamental window pattern, is a second pattern. Prewhitening at the orbital frequency gives a clear pattern with an amplitude of 0.07 mag centred on $`P_sh=2.980`$ h, which is the period of the drifting humps. Therefore, V630 Sgr is a nova remnant that is also an eclipsing permanent superhumper. This is the first CV found to combine all of these features. The superhump excess of 5.3% and beat period of 2.4 d are typical for superhumps in superoutbursting dwarf novae and novalike variables (Warner 1995a). The beat period is evident in the night-to-night variation of mean brightness, as seen in Fig. 7. The Fourier transform of this delivers an approximate period of 2.1 d and an amplitude of 0.22 mag.

The eclipse depth of $`0.6`$ mag is relatively shallow, but at the $`P_{orb}`$ of V630 Sgr it implies a quite high orbital inclination. The width at half depth of the narrowest eclipses is $`\mathrm{\Delta }\psi 0.085P_{orb}`$. The secondary mass $`M(2)`$ at $`P_{orb}=2.8`$ h is $`0.24`$ M (Warner 1995a) and as V630 Sgr was a very fast nova we expect $`M(1)1.2`$ M, giving a mass ratio $`q0.20`$. These figures are not quite compatible with the theoretical $`\mathrm{\Delta }\psi ,q,i`$ relationship (Horne 1985), but they serve to show that $`i85^{}`$. The presence of superhumps shows V630 Sgr to be a high $`\dot{M}`$ system which will have a disc radius at the 3:1 resonant radius, which is $`r_d=0.45a`$. The Roche lobe radius for the secondary is $`R(2)=0.25a`$ for $`q=0.2`$. Therefore, even at large inclination, the eclipses are partial, as is seen in the eclipse profiles in Fig. 5.

The observed eruption range of 17.5 mag is 3.5 mag greater than expected for a $`t_2=4`$ d nova (see Fig. 5.4 of Warner 1995a). Using the correction formula (Eqn. 2.63 of Warner 1995a) for disc inclinations, removes the 3.5 mag excess if $`i88^{}`$. We conclude that V630 Sgr is a very high inclination system.

### 2.6 YY Telescopium

YY Tel is a dwarf nova with a range V $`14.419.3`$ (Downes et al. 1997) suspected to be an SU UMa star of long outburst interval (O’Donoghue et al. 1991). The last reported outburst was in October 1998, and no high speed photometry has been reported, either in outburst or at quiescence. The candidate identified at quiescence (Downes et al. 1997) was observed spectroscopically by Zwitter & Munari (1996) and showed only a continuum with no emission lines characteristic of a CV. Their photometry gave V = 19.26, (B-V) = 0.47, (U-B) = -0.28, placing the candidate near the black body line at $`T6700`$ K in the two colour diagram. This is also a region frequented by late F type subdwarfs.

Our observation on 5 June 2000 showed the YY Tel candidate to be constant in brightness (see Fig. 8). Re-examination nearly one year later again showed absence of any rapid variation but the star changed in brightness by 0.15 mag over three to five days. A nearby, fainter star, marked in Fig. 9 (YY Tel is the brighter of the two stars that are marked), shows marginal evidence for variation and should be examined spectroscopically.

### 2.7 CQ Velorum

CQ Vel was Nova Velorum 1940, reaching V $`9.0`$. It has been identified with a 21st magnitude star (Duerbeck 1987) but no confirmatory spectra or photometry have been obtained. In our observations on 10 March 2000 we found the candidate nova remnant to be of constant brightness as were all the stars in the CCD frame, except for a star slightly fainter than the candidate, 9′′ to the northwest (Fig. 10). This star has the flickering characteristic of a CV, see Fig. 11. These observations provide a correct identification for CQ Vel.

### 2.8 CE-315

CE-315 is a recently discovered member of the AM CVn class of CVs: systems transferring helium rather than hydrogen (Ruiz et al. 2001). It is only the second of this class to be found in a long-term state of low $`\dot{M}`$, the other being GP Com (alias G61-29). Ruiz et al. found a period of 65.1 $`\pm `$ 0.7 min from an S-wave in the spectra, which by analogy with GP Com is thought to be the orbital period.

We have obtained two light curves for CE-315, which are displayed in Fig. 12. There is no sign in the light curve, or its Fourier Transform, of the orbital periodicity. CE-315 resembles GP Com in the general appearance of its light curve (see Warner 1972; Harrop-Allin 1996): much short time scale flaring or flickering set on a slowly varying background. In GP Com there are occasions when a few cycles at the 46.5 min orbital period appear in the light curve (Harrop-Allin 1996); if the same occurs in CE-315 it would need extensive photometry to detect.

The fact that the two longest period AM CVn stars, GP Com and CE-315, are low $`\dot{M}`$ systems (equivalent to dwarf novae in hydrogen transferring CVs) is in accord with the evolutionary model where orbital angular momentum is the result of radiation of gravitational waves (Warner 1995a,b). According to Eqn. 9.61 of Warner (1995a) the mass transfer rates should be $`\dot{M}1.0\times 10^{11}`$ M y-1 for GP Com and $`\dot{M}3\times 10^{12}`$ M y-1 for CE-315. GP Com has X-rays modulated at the orbital period (Beuermann & Thomas 1993). At B $`17.2`$, CE-315 is $`1.4`$ mag fainter than GP Com, but this is accounted for by the lower $`\dot{M}`$, so they are at similar distances and it would be worth seeking an X-ray detection.

## 3 Conclusions

The results presented here cover a miscellany of CV types and demonstrate the rewards waiting among the fainter members of this class. For those objects where new phenomena have been found (e.g., the possible quiescent superhumps in V359 Cen), or where rarity is an incentive (e.g., the eclipsing superhumping nova remnant V630 Sgr), further observations with larger telescopes are required. V351 Pup may show circular polarization or other magnetic signatures when it has descended the nova light curve further.

## Acknowledgments

We thank Dr. D. O’Donoghue for the use of his EAGLE program for Fourier analysis of the light curves. PAW is funded partly through strategic funds made available to BW by the University of Cape Town. BW’s research is funded entirely by that university.


text:
# Elastic critical behaviour in a 3𝑑 model for polymer gels

## I Introduction

The elastic response in polymer gels is due to the presence of the polymeric network formed in the gelation process: this random structure makes the system able to respond to external stresses. Due to the large possibility of conformational changes which characterizes this macromolecular phase, there is an entropic term contributing to the elastic properties of the system. This picture is typical of gelling systems and may fit a wide range of very different materials. Then because of the different role played by the entropic term and the particular features of the polymeric networks, a rich phenomenology is produced.
The viscoelastic behaviour of gelling materials and the mechanisms at its origin certainly represent a central issue in soft matter physics and have a high relevance to a wide range of applications from food processing to materials science. On the other hand these materials show many typical features of complex systems and their study is connected to some fundamental problems in soft matter physics, from the entropic elastic behaviour, which we are interested in here, to glassy dynamics. As a consequence these systems are intensively investigated both experimentally and by means of statistical mechanics models, bringing up a lively debate.
In particular here we study the critical behaviour of the elastic response of the system as the gelation process takes place producing the polymeric network: in a gelling system the elastic modulus starts growing at the gelation threshold with a power law behaviour, usually expressed in terms of the polymerization degree. In the experiments performed on different gelling materials the value of the critical exponent $`f`$ describing the critical behaviour of the elastic modulus appears to be close to either $`f2`$, or $`f3`$, or else $`f4`$. The value $`f2`$ is actually observed in experiments on agarose gel , gelatin gels , some silica gels (TEOS) . Within statistical mechanics models it corresponds to the prediction based on the de Gennes’ analogy between the elasticity of a percolating random network of Hookean springs and the conductivity in a random percolating network of resistors . The critical exponent $`f3.`$ has been observed in diisocianate/triol gel , in epoxy resins , TMOS silica gels , polyesters gels . It is very close to the value $`f=d\nu `$, where $`d`$ is the dimensionality of the system and $`\nu `$ is the critical exponent of the connectedness length $`\xi `$, predicted in ref. which in $`3d`$ gives $`f2.64`$ (in mean field both $`d\nu `$ and the random resistor network exponent are equal to $`3`$). Following the argument in ref. this exponent would then correspond to the case of a dominating entropic term, evaluated by means of scaling arguments for the percolating network: the elastic energy of the system is due to the entropic elasticity of the macrolink whose length is of the order of the connectedness correlation length. Finally the value $`f4`$ can be linked to the prediction of the bond-bending model $`f=d\nu +1`$  or alternatively $`f=t+2\nu `$ , which in $`3d`$ give $`f3.7`$. This would imply that within the elastic energy describing the system there is a bending term playing a relevant role in the elastic response of the network. This critical behaviour is typical of some colloidal gels . The clustering of the experimental values for the elasticity critical exponent in gels around discrete values suggests the possibility of correspondently individuating different universality classes, which should be characterized by some intrinsic features of the networks formed in the materials.
Recent numerical studies via molecular dynamics  and Monte Carlo simulations  of percolating networks of tethered particles with no hard core interactions have shown that the shear modulus critical exponent is $`1.3`$ in $`d=2`$ and $`2.`$ in $`d=3`$. Monte Carlo simulations of two and three-dimensional percolating networks of tethered particles with hard core repulsion  find consistent results for the shear modulus critical exponent. These results agree with the de Gennes prediction and with some experimental results.
Within the numerical studies we have approached the study of this problem introducing a percolation dynamic model, and directly investigating the dynamic viscoelastic properties as the percolation transition takes place. Our model introduces in the percolation model the bond-fluctuation dynamics, which takes into account the conformational changes of the polymer molecules and the excluded volume interactions. This model has been translated in a lattice algorithm and studied via numerical simulations on hypercubic lattices. Actually it presents many fundamental features of the gelation phenomenology and has already allowed to study the critical behaviour of the viscoelastic properties and the relaxation process in gelling systems in the sol phase . Numerical simulations of the model in $`d=2`$  have shown that the elasticity critical exponent is $`2.7`$, a value which agrees with the prediction $`f=d\nu `$ in . The study of the model in $`d=3`$ is of great interest to check the elasticity critical behaviour and to eventually compare the findings with the experimental results. We here present the results of large scale $`3d`$ simulations: the data show that $`f2.5\pm 0.1`$, again in agreement with the prediction $`f=d\nu `$, coherently with the findings in $`d=2`$.
In the next section we present the model and the details of the numerical simulations, in section III the elastic response of the percolating cluster is discussed and a scaling behaviour is obtained; in section IV the results of the numerical simulations are presented and discussed; section V contains concluding remarks.

## II Model and numerical simulations

We consider a solution of tetrafunctional monomers of concentration $`p`$. The monomers interact via excluded volume interactions, i.e. a monomer occupies a lattice elementary cell and two occupied cells cannot have common sites. Nearest neighbours or next-nearest neighbours are instantaneously linked by a permanent bond with probability $`p_b`$. In terms of these two parameters the percolation line can be determined via the critical behaviour of the percolation properties, individuating the sol and the gel phase in the system. Actually in the simulation we fix $`p_b=1`$ and study the system varying $`p`$ . The percolation quantities critical exponents agree with the random percolation predictions  (e.g. $`\gamma 1.8\pm 0.05`$ and $`\nu 0.89\pm 0.01`$ in $`3d`$ ). The monomers free or linked in clusters diffuse via random and local movements on the lattice according to the bond-fluctuation dynamics , which is ruled by the possibility of varying the bond lengths within a set of values determined by the excluded volume interactions and the SAW condition. This produces a high number of different bond vectors and we consider the case of permanent bonds, which corresponds to the strong gelation process. In $`3d`$ the allowed bond lengths on the cubic lattice are $`l=2,\sqrt{5},\sqrt{6},3,\sqrt{10}`$ and in Fig.1 an example of different allowed configurations for a polymer molecule is shown.
We here present the results of extended numerical simulations of the model in the gel phase to study the elastic response in the system. It is worth noticing that there is no elastic potential energy for the bond vectors and then the elastic behaviour is purely entropic implying that our study is performed at finite temperature.
The data presented here refer to lattice sizes $`L`$ ranging from $`12`$ to $`32`$, and have been averaged over $`30`$ different realizations. The simulations have been performed on the $`CRAYT3E`$ system at CINECA taking more than $`30000`$ hours/node of CPU time.

## III Elastic response in the gel phase

We study the elastic response in the gel phase in terms of the macroscopic elastic constant of the system $`K`$, which is experimentally defined as the ratio between an applied external force and the deformation. In a simple elongation experiment if $`l_0`$ is the undeformed length and $`\delta =(ll_0)`$ is the deformation in the system, within the linear response approximation the elastic free energy $`FK\delta ^2`$. In terms of the Young elastic modulus $`E`$ the free energy per unit volume is $`F/VE\delta ^2/ł_0^2`$. Then $`KEV/l_0^2`$ and for a cube of size $`L`$, $`KEL^{d2}`$, expressing the fact that the elastic modulus has the dimensions of an energy per unit volume and is an intensive quantity whereas $`K`$ depends on the system size $`L`$.
In the gel, since $`E`$ vanishes at $`p_c`$ as $`\xi ^{\stackrel{~}{f}}`$ (where $`\stackrel{~}{f}=f/\nu `$) one has

$$KL^{d2}\xi ^{\stackrel{~}{f}}$$
(1)
Following eq.(1) the macroscopic elastic constant presents the corresponding scaling behaviours as function of the system size $`L`$ and of the distance from the percolation threshold $`(pp_c)`$:

$$\begin{array}{ccccc}p>p_c\hfill & & & & KL^{d2}\hfill \\ p=p_c\hfill & & & & KL^{\stackrel{~}{z}}\hfill \\ \text{fixed}L\hfill & & & & K(pp_c)^f\hfill \end{array}$$
(2)
where $`\stackrel{~}{z}=\stackrel{~}{f}(d2)`$.
An alternative way to obtain these scaling relations is to consider the percolating cluster as a network of nodes connected by macrolinks of linear size $`\xi `$ . Each macrolink can be considered as a spring with an effective elastic constant. At the percolation threshold there is only one macrolink spanning the system and its effective elastic constant coincides with the system macroscopic elastic constant $`K`$.
In order to evaluate $`K`$ we notice that for a spring of elastic constant $`K`$ in a thermal bath at temperature $`T`$, the mean fluctuation in the energy $`U`$ is $`\mathrm{\Delta }U=\frac{1}{2}Kx^2`$, where $`x`$ is the spring elongation. From the energy equipartition $`Kx^2=k_BT`$, so that at the equilibrium the elastic constant $`K`$ is related to the the fluctuations in the spring length . This result can be more generally obtained by means of the Fokker-Plank equation for the probability distribution of the spring elongation .
Therefore the macroscopic elastic constant of the gel phase is related to the fluctuations of the linear size of the infinite cluster, i.e. the squared fluctuations in the gyration radius of the percolating cluster $`\mathrm{\Delta }R_g^2=(R_gR_g)^2`$, $`K1/\mathrm{\Delta }R_g^2`$.
An alternative way to calculate the macroscopic elastic constant in the gel phase by means of the fluctuations in the unperturbed system is to consider the center of mass of the percolating cluster as a brownian particle, subject to a restoring force responsible for the elastic behaviour of the system. The restoring force introduces a limitation on the diffusive process of a brownian particle. Using the same argument based on the energy equipartition, the asymptotic equilibrium value $`\mathrm{\Delta }`$ of its displacement fluctuations $`\mathrm{\Delta }R^2(t)`$ is inversely proportional to the elastic constant, i.e. $`\mathrm{\Delta }1/K`$.

## IV Results

In order to numerically study the elastic response we have used both the approaches mentioned before namely we have calculated the average fluctuation of the gyration radius of the percolating cluster $`\mathrm{\Delta }R_g^2`$ and the asymptotic value $`\mathrm{\Delta }`$ of the mean square displacement of its center of mass $`\mathrm{\Delta }R^2(t)`$ .
In the first approach the average fluctuation in the gyration radius $`\mathrm{\Delta }R_g^2`$ of the percolating cluster has been computed at the percolation threshold in system of different size $`L`$ with periodic boundary conditions. In Fig.2 $`\mathrm{\Delta }R_g^2`$ at $`p_c`$ is shown in a log-log plot as function of the lattice size $`L`$. In the considered range the data are well fitted by a power law behaviour according to eq.(2), giving a critical exponent $`\stackrel{~}{z}1.9\pm 0.1`$, i.e. $`z1.7\pm 0.1`$. As $`\stackrel{~}{z}=\stackrel{~}{f}(d2)`$, $`f=z+(d2)\nu `$ and this result gives $`f2.6\pm 0.1`$.
In the second approach we have calculated $`\mathrm{\Delta }R^2(t)`$ at different steps of the gelation process, i.e. as $`p`$ grows above the percolation threshold $`p_c0.718\pm 0.005`$ , and at $`p_c`$ for different lattice sizes. In these simulations hard-wall boundary conditions have been used . In the gel at the gelation transition the percolating cluster is a quite loose network, the center of mass is rather free and the elastic response is weak. As the gelation process goes on, the network tightens becoming more rigid, the elastic constant of the system increases and the center of mass motion is progressively constrained. This would be then the physical mechanism producing the critical behaviour of the elastic response for a critical gel. In agreement with our picture $`\mathrm{\Delta }R^2(t)`$ grows with time up to a limiting plateau value $`\mathrm{\Delta }`$ as it is shown in Fig.3. This quantity is inversely proportional to the elastic constant of the system $`\mathrm{\Delta }1/K`$ and increases as the percolation threshold is approached from above.
In Figures 4,5 and 6 the scaling behaviour obtained for $`\mathrm{\Delta }`$ is presented.
In Fig.4 $`\mathrm{\Delta }(L,p)`$ for $`p>>p_c`$ ($`p=0.85`$) is shown in a double logarithimic plot: being far from the percolation threshold the system is reasonably homogeneous and in fact the scaling behaviour as $`\mathrm{\Delta }L^{0.99\pm 0.1}`$ is observed, in agreement with the behaviour $`KL^{d2}`$.
In Fig.5 $`\mathrm{\Delta }(L,p=p_c)`$ is shown in a double logaritmic plot as a function of the lattice size $`L`$: the data exhibit a behaviour $`\mathrm{\Delta }L^{\stackrel{~}{z}}`$ with $`\stackrel{~}{z}2.0\pm 0.1`$. This result again gives $`f2.6\pm 0.1`$.
Finally in Fig.6 $`\mathrm{\Delta }(pp_c)`$ is plotted for the lattice size $`L=32`$: we fit the data with a power law behaviour (eq.(2)) and obtain a critical exponent $`f2.5\pm 0.1`$.
It is straightforward to notice that all these numerical results can be coherently interpreted in terms of the scaling relations obtained for $`K`$.
The value of the critical exponent $`f2.6`$ is in good agreement with the prediction $`f=d\nu `$ of ref., therefore supporting the picture there proposed, and consistent with the value obtained in the $`2d`$ study of the model .
Due to the limited extension of the critical parameter $`(pp_c)`$ and of $`L`$ here investigated the eventual occurrence of a crossover to a different exponent cannot be excluded.

## V Conclusions

The numerical results of Figures 2,5 and 6 show that $`\stackrel{~}{z}=2.`$, coherently agreeing with the prediction $`f=d\nu `$. On the whole they support the scaling picture we propose and the argument of ref.. They also agree with some experimental results -. This result has been obtained via two independent calculations giving consistent numerical values and is also consistent with the value previously obtained in the $`2d`$ study.
On the other hand the recent numerical works on entropic elastic models of refs. find a good agreement with the de Gennes prediction. These results, together with the experimental data, seem to suggest the possibility that there are two distinct universality classes characterized one by an exponent $`f=d\nu `$ and another by the electrical analogy exponent $`f=t`$, which in $`3d`$ are respectively $`2.64`$ and $`2.`$. However since the models in the different numerical studies are rather similar the possibility of a crossover between different dynamic regimes, as it is observed in some experiments , cannot be completely excluded. In both cases these results give a hint for the interpretation of the experimental data and indicate the aspects which should be further investigated.
This work has been partially supported by the European TMR Network-Fractals under contract No.FMRXCT980183, by INFM PRA-HOP 99, by MURST-PRIN-2000, by the INFM Parallel Computing Initiative and by the European Social Fund.


text:
# References

Recently, there has been active interest in the possibility that our four-dimensional universe may be a three-brane embedded in higher-dimensional spacetime , as such brane world scenarios might provide possible solutions to the hierarchy and the cosmological constant problems. In the scenario proposed by Randall and Sundrum (RS) , the extra spatial dimension needs not be compact for reproducing four-dimensional gravity on the brane due to the strongly warped bulk spacetime and the fine-tuned brane tension is required for solving the hierarchy and the cosmological constant problems.

In particular, the brane universes in the bulk of the five-dimensional AdS black holes  and higher-dimensional branes  in string theories have been studied. In such models, the motion of the probe brane in the bulk of the source brane is responsible for the expansion of the brane universe on the probe brane. When the equation of state for the matter field on the brane takes a special form, the static brane configuration, in which the probe brane remains static in the bulk, is possible. Such static brane configurations may be regarded as generalizations of the static RS model  to different bulk spacetimes. In Refs. , it was found out that such static brane configuration is possible even without matter fields on the probe brane for some bulk spacetimes. Such result was generalized to the case of the topological AdS black holes in Ref. .

In this paper, we study the static brane configurations in the bulk of the topological black holes in asymptotically flat spacetime. Contrary to the previous belief that black holes in asymptotically flat spacetime should have spherical horizon , new class of dilatonic $`p`$-brane solutions in asymptotically flat spacetime where horizon is the maximally symmetric space with any values of the curvature parameter was constructed . In particular, uncharged 0-brane solution with flat horizon can be transformed to the self-tuning flat domain wall solution  when the dilaton coupling parameter takes a special value. We find that the static brane configuration is possible for any horizon geometries, unlike the case of the topological AdS black holes studied in Ref. . We obtain the fine-tuned value of the brane tension and the stabilized location of the probe brane in terms of the parameters of the bulk black hole solution. We also study massless scalar field in such bulk background.

We begin by discussing the new class of dilatonic $`p`$-brane solutions constructed in Ref. . For the purpose of putting the solutions into more convenient forms, we choose to parameterize the action for the $`p`$-brane solutions in the following form:

$$S_p=\frac{1}{2\kappa _D^2}d^Dx\sqrt{G}\left[\frac{4}{D2}(\varphi )^2\frac{1}{2(p+2)!}e^{2a_p\varphi }F_{p+2}^2\right],$$
(1)
where $`\kappa _D^2`$ is the $`D`$-dimensional gravitational constant, $`\varphi `$ is the dilaton field, $`F_{p+2}`$ is the field strength of a $`(p+1)`$-form potential $`A_{p+1}`$ and $`a_p`$ is the dilaton coupling parameter. The most general form of the action with arbitrary coefficients in front of the kinetic terms, considered in Ref. , is related to this action just through rescaling of $`\varphi `$ and $`A_{p+1}`$ with constant factors. In terms of such new parametrization, the solution for the topological non-extreme dilatonic $`p`$-brane with the longitudinal coordinates $`𝐱=(x^1,\mathrm{},x^p)`$, constructed in Ref. , takes the following form:

$`G_{MN}dx^Mdx^N`$ $`=`$ $`h_{}^{\frac{4(n1)}{(p+n)\mathrm{\Delta }_p}1}h_+dt^2+h_{}^{\frac{4(n1)}{(p+n)\mathrm{\Delta }_p}}d𝐱d𝐱+h_{}^{\frac{4(p+1)}{(p+n)\mathrm{\Delta }_p}\frac{n3}{n1}}h_+^1dr^2`$ (3)
$`+h_{}^{\frac{4(p+1)}{(p+n)\mathrm{\Delta }_p}+\frac{2}{n1}}r^2ds_{n,k}^2,`$
$`e^{2\varphi }`$ $`=`$ $`h_{}^{\frac{(p+n)a_p}{\mathrm{\Delta }_p}},A_{tx^1\mathrm{}x^p}={\displaystyle \frac{2}{\sqrt{\mathrm{\Delta }_p}}}{\displaystyle \frac{m\mathrm{sinh}\alpha \mathrm{cosh}\alpha }{r^{n1}}},`$ (4)
$`h_+`$ $`=`$ $`s(r)\left(1{\displaystyle \frac{m\mathrm{cosh}^2\alpha }{r^{n1}}}\right),h_{}=\left|k{\displaystyle \frac{m\mathrm{sinh}^2\alpha }{r^{n1}}}\right|,`$ (5)
$`\mathrm{\Delta }_p`$ $`=`$ $`{\displaystyle \frac{(p+n)a_p^2}{2}}+{\displaystyle \frac{2(p+1)(n1)}{p+n}},`$ (6)
where $`nDp2`$, $`s(r)\mathrm{sgn}(km\mathrm{sinh}^2\alpha /r^{n1})`$ and $`ds_{n,k}^2`$ is the metric for the maximally symmetric $`n`$-dimensional space with constant curvature, which we parametrize as

$$ds_{n,k}^2=\left(1+\frac{k}{4}\delta _{ij}y^iy^j\right)^2[(dy^1)^2+\mathrm{}+(dy^n)^2],$$
(7)
and $`k=0,\pm 1`$ is the curvature parameter. With our parametrization it is manifest that the curvature singularity at $`r=r_{}=(m\mathrm{sinh}^2\alpha )^{1/(n1)}`$ for the $`k=1`$ case is always covered by a horizon at $`r=r_+=(m\mathrm{cosh}^2\alpha )^{1/(n1)}`$, provided $`m>0`$. In this paper, we consider the $`p=0`$ case, for which the solution takes the form:

$`G_{MN}dx^Mdx^N`$ $`=`$ $`h_{}^{\frac{4(n1)}{n\mathrm{\Delta }}1}h_+dt^2+h_{}^{\frac{4}{n\mathrm{\Delta }}\frac{n3}{n1}}h_+^1dr^2+h_{}^{\frac{4}{n\mathrm{\Delta }}+\frac{2}{n1}}r^2ds_{n,k}^2,`$ (8)
$`e^{2\varphi }`$ $`=`$ $`h_{}^{\frac{na}{\mathrm{\Delta }}},A_t={\displaystyle \frac{2}{\sqrt{\mathrm{\Delta }}}}{\displaystyle \frac{m\mathrm{sinh}\alpha \mathrm{cosh}\alpha }{r^{n1}}},`$ (9)
$`h_+`$ $`=`$ $`s(r)\left(1{\displaystyle \frac{m\mathrm{cosh}^2\alpha }{r^{n1}}}\right),h_{}=\left|k{\displaystyle \frac{m\mathrm{sinh}^2\alpha }{r^{n1}}}\right|,`$ (10)
$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{na^2}{2}}+{\displaystyle \frac{2(n1)}{n}},`$ (11)
where now $`D=n+2`$ and $`a:=a_0`$ is the dilaton coupling parameter for the $`U(1)`$ gauge field $`A_M`$. We may regard this topological 0-brane as describing bulk spacetime of a brane world with the worldvolume coordinates $`(t,y^1,\mathrm{},y^n)`$ and the transverse coordinate $`r`$. However, unlike the RS model , such brane world model generally has asymmetrically warped spacetime, meaning that the space and the time coordinates have different warp factors, and as a result the Lorentz invariance is violated . However, for the uncharged 0-brane solution with $`k=0`$, the bulk metric can become $`(D1)`$-dimensional Poincaré invariant. The uncharged solution in Ref.  takes the following form in our parametrization:

$`G_{MN}dx^Mdx^N`$ $`=`$ $`h_{}^{\frac{4(n1)}{n\mathrm{\Delta }}1}h_+dt^2+h_{}^{\frac{4}{n\mathrm{\Delta }}\frac{n3}{n1}}h_+^1dr^2+h_{}^{\frac{4}{n\mathrm{\Delta }}+\frac{2}{n1}}r^2ds_{n,k}^2,`$ (12)
$`e^{2\varphi }`$ $`=`$ $`h_{}^{\frac{na}{\mathrm{\Delta }}},h_+=\mathrm{sgn}h,h_{}=|h|,h=k{\displaystyle \frac{l}{r^{n1}}},`$ (13)
where $`l`$ is an integration constant. Note, this uncharged $`p`$-brane solution cannot be obtained by taking the zero charge limit of the charged $`p`$-brane solution (6) 222This fact does not imply that our convenient parametrization (6) for the charged topological $`p`$-brane solution is wrong. This can be easily seen by the fact that the usual uncharged black $`p`$-brane solution (with $`k=1`$) cannot be coordinate-transformed to take the form (13) with $`k=1`$.. With a choice of $`l=1`$ and $`k=0`$, the metric (13) has the $`(D1)`$-dimensional Poincaré invariance, when the temporal and the spatial warp factors have the same $`r`$-dependence:

$$r^{(n1)\left[\frac{4(n1)}{n\mathrm{\Delta }}1\right]}=r^2r^{(n1)\left[\frac{4}{n\mathrm{\Delta }}+\frac{2}{n1}\right]}\mathrm{\Delta }=4.$$
(14)
For such case, the uncharged solution (13) becomes the self-tuning domain wall solution constructed in Refs. , after the redefinition of the transverse coordinate $`r`$ . So, we see that the self-tuning domain wall solution can be embedded in string theories as an uncharged 0-brane solution with flat ($`k=0`$) horizon and $`\mathrm{\Delta }=4`$. Such uncharged 0-brane solution can be obtained by compactifying, for example, uncharged D6-brane, for which $`\mathrm{\Delta }_6=4`$, on a compact Ricci flat manifold along its longitudinal directions.

In this paper, we study the brane world on a probe $`n`$-brane in the bulk background of the topological dilatonic 0-brane (11). For such purpose, it is convenient to redefine the transverse coordinate in the following way :

$$dz=h_{}^{\frac{2}{n\mathrm{\Delta }}\frac{n3}{n1}}h_+^{\frac{1}{2}}dr,$$
(15)
so that the bulk metric (11) takes the form that resembles the RS domain wall metric except that the warp factor is asymmetric:

$$G_{MN}dx^Mdx^N=A(z)^2dt^2+B(z)^2\sigma _{ij}(x)dx^idx^j+dz^2,$$
(16)
where $`A(z)`$ and $`B(z)`$ are implicitly given by

$$A(z)=h_{}^{\frac{2(n1)}{n\mathrm{\Delta }}1}h_+^{\frac{1}{2}},B(z)=h_{}^{\frac{2}{n\mathrm{\Delta }}+\frac{1}{n1}}r,$$
(17)
with $`r`$ related to $`z`$ in the way defined in Eq. (15). We consider the following form of the $`n`$-brane action:

$$S_\sigma =d^{n+1}x\sqrt{g}f(\varphi ),$$
(18)
where $`g`$ is the determinant of the induced metric $`g_{\mu \nu }=\delta _\mu ^M\delta _\nu ^NG_{MN}`$ on the $`n`$-brane and $`f(\varphi )`$ is an arbitrary function of $`\varphi `$.

We consider the static brane configuration without matter fields on the $`n`$-brane. From the total action $`S_{p=0}+S_\sigma `$, we obtain the following Einstein’s equations and equation of motion for the dilaton field:

$`𝒢_{MN}`$ $`=`$ $`\frac{2}{D2}G_{MN}(\varphi )^2+\frac{4}{D2}_M\varphi _N\varphi +\frac{1}{2}e^{2a\varphi }(F_{MP}F_N^P\frac{1}{4}G_{MN}F_2^2)`$ (20)
$`\kappa _D^2\sqrt{{\displaystyle \frac{g}{G}}}\delta _M^\mu \delta _N^\nu g_{\mu \nu }f(\varphi )\delta (zz_b),`$
$$\frac{4}{D2}_M\left[\sqrt{G}G^{MN}_N\varphi \right]=\frac{a}{4}e^{2a\varphi }F_2^2+\kappa _D^2\sqrt{g}f^{}(\varphi )\delta (zz_b),$$
(21)
where $`𝒢_{MN}=_{MN}\frac{1}{2}G_{MN}`$ is the Einstein tensor for the bulk metric $`G_{MN}`$, the constant $`z_b`$ is the location of the $`n`$-brane and $`D=n+2`$. After the explicit expression for the bulk metric (16) is substituted into these equations of motion, the $`(z,z)`$-, $`(t,t)`$\- and $`(i,j)`$-components of the Einstein’s equations and the dilaton equation of motion respectively take the following forms:

$`{\displaystyle \frac{(D2)(D3)}{2}}{\displaystyle \frac{k}{B^2}}+(D2){\displaystyle \frac{A^{}B^{}}{AB}}+{\displaystyle \frac{(D2)(D3)}{2}}{\displaystyle \frac{B^{\mathrm{\hspace{0.17em}2}}}{B^2}}`$ (22)
$`={\displaystyle \frac{2}{D2}}(_z\varphi )^2{\displaystyle \frac{3}{8}}{\displaystyle \frac{1}{A^2}}e^{2a\varphi }F_{tz}^2,`$ (23)
$`{\displaystyle \frac{(D2)(D3)}{2}}{\displaystyle \frac{k}{B^2}}(D2){\displaystyle \frac{B^{\prime \prime }}{B}}{\displaystyle \frac{(D2)(D3)}{2}}{\displaystyle \frac{B^{\mathrm{\hspace{0.17em}2}}}{B^2}}`$ (24)
$`={\displaystyle \frac{2}{D2}}(_z\varphi )^2+{\displaystyle \frac{3}{8}}{\displaystyle \frac{1}{A^2}}e^{2a\varphi }F_{tz}^2+\kappa _D^2f(\varphi )\delta (zz_b),`$ (25)
$`{\displaystyle \frac{D3}{2}}\left[(D4){\displaystyle \frac{k}{B^2}}+(D4){\displaystyle \frac{B^{\mathrm{\hspace{0.17em}2}}}{B^2}}+2{\displaystyle \frac{B^{\prime \prime }}{B}}+2{\displaystyle \frac{A^{}B^{}}{AB}}\right]+{\displaystyle \frac{A^{\prime \prime }}{A}}`$ (26)
$`={\displaystyle \frac{2}{D2}}(_z\varphi )^2\kappa _D^2f(\varphi )\delta (zz_b),`$ (27)
$$\frac{4}{D2}_z(AB^n_z\varphi )=\frac{a}{4}\frac{1}{A^2}e^{2a\varphi }F_{tz}^2+\kappa _D^2AB^nf^{}(\varphi )\delta (zz_b),$$
(28)
where the prime denotes the derivative w.r.t. the transverse coordinate $`z`$.

We now construct the brane world model with an orbifold structure $`S^1/𝐙_2`$ just like the RS1 model . Following the procedure proposed in Ref. , we cut off the transverse space at the black hole horizon $`z=z_H`$ (corresponding to $`r=r_+=(m\mathrm{cosh}^2\alpha )^{1/(n1)}`$) and at the location $`z=z_b`$ of the $`n`$-brane. The spacetime between these boundaries is referred to as the bulk. We introduce a new transverse coordinate $`\overline{z}z_bz`$ and identify the corresponding two points under the $`𝐙_2`$-symmetry $`\overline{z}\overline{z}`$. So, in this new coordinate the $`n`$-brane is located at $`\overline{z}=0`$ and the black hole horizon is located at orbifold points $`\overline{z}=\pm |z_bz_H|`$, which are identified under the $`𝐙_2`$-symmetry $`\overline{z}\overline{z}`$. We rewrite the bulk metric (16) in terms of the new coordinates as

$$G_{MN}dx^Mdx^N=A(|\overline{z}|)^2dt^2+B(|\overline{z}|)^2\sigma _{ij}(x)dx^idx^j+d\overline{z}^2,$$
(29)
where $`A`$ and $`B`$ are now functions of $`|\overline{z}|`$ due to the imposed $`𝐙_2`$-symmetry and $`|z_bz_H|\overline{z}|z_bz_H|`$. Although the metric components are continuous at $`\overline{z}=0`$, their first and second derivatives w.r.t. $`\overline{z}`$ are discontinuous and have $`\delta `$-function singularity at $`\overline{z}=0`$, respectively.

By applying the boundary conditions on the first derivatives at $`z=z_b`$ resulting from the $`\delta `$-function terms in Eqs. (25,27,28), we can fix the parameters in the $`n`$-brane action (18) and the (stabilized) location $`r_b`$ of the $`n`$-brane. For the phenomenological relevance, from now on we consider the $`k=0`$ case, only 333The nontrivial static brane configuration is possible also when $`k=\pm 1`$. The below equations (30) and (31) continue to hold even for the $`k=\pm 1`$ cases. However, the stabilized brane location $`r_b`$ take more complicated forms $`r_b^{n1}=m(2n2+n\mathrm{\Delta }2\mathrm{\Delta })\mathrm{cosh}^2\alpha \mathrm{sinh}^2\alpha /[(2n2+n\mathrm{\Delta }2\mathrm{\Delta })\mathrm{cosh}^2\alpha 2n+2+\mathrm{\Delta }]`$ and $`r_b^{n1}=m(2n2+n\mathrm{\Delta }2\mathrm{\Delta })\mathrm{cosh}^2\alpha \mathrm{sinh}^2\alpha /[(2n2n\mathrm{\Delta })\mathrm{cosh}^2\alpha 2n+2+\mathrm{\Delta }]`$, respectively for $`k=1`$ and $`k=1`$.. By applying the boundary conditions resulting from Eqs. (25,28) we obtain the following constraints on the parameter(s) of $`f(\varphi )`$:

$$\kappa _D^2f(\varphi _b)=\frac{4(n1)}{\mathrm{\Delta }}r_b^1h_{}^{\frac{2}{n\mathrm{\Delta }}+\frac{n3}{n1}}(r_b)h_+^{\frac{1}{2}}(r_b),$$
(30)
$$\kappa _D^2f^{}(\varphi _b)=\frac{4(n1)a}{\mathrm{\Delta }}r_b^1h_{}^{\frac{2}{n\mathrm{\Delta }}+\frac{n3}{n1}}(r_b)h_+^{\frac{1}{2}}(r_b),$$
(31)
respectively, where $`\varphi _b\varphi (r_b)`$. We note that the following choice of $`f(\varphi )`$ trivially satisfies these constraints:

$$f(\varphi )=\sigma e^{a\varphi },$$
(32)
where $`\sigma `$ is the $`n`$-brane tension. The consistency of the boundary conditions resulting from Eqs. (25,27) fixes the stabilized location of the $`n`$-brane to be

$$r_b^{n1}=\frac{4\mathrm{\Delta }}{42\mathrm{\Delta }}m\mathrm{cosh}^2\alpha .$$
(33)
From this, we see that a real-valued brane position $`0r_b<r_+`$ is possible when $`0<\mathrm{\Delta }<2`$ or $`\mathrm{\Delta }4`$ with $`m>0`$. This result is in contrast to the case of the uncharged topological AdS black hole, in which the static brane configuration with $`k=0`$ is possible only when the black hole mass is zero . We briefly discuss the limiting cases. When $`\mathrm{\Delta }=4`$, the stable position of the $`n`$-brane is at $`r_b=0`$, namely at the singularity. When $`\mathrm{\Delta }=0`$, the stable position of the $`n`$-brane is at the horizon $`r_b^{n1}=m\mathrm{cosh}^2\alpha `$, meaning that the range of the transverse coordinate $`\overline{z}`$ is zero, i.e., this case is degenerate. When $`\mathrm{\Delta }=2`$ the static configuration is not possible, since $`r_b=\mathrm{}`$ for such case. For the bulk background of the uncharged 0-brane (13) with $`k=0`$, the consistency of the boundary conditions resulting from Eqs. (25,27) requires that $`\mathrm{\Delta }=2(n1)`$, meaning that as long as $`\mathrm{\Delta }=2(n1)`$ the static brane configuration for any values of the $`n`$-brane location $`r_b`$ is possible. So, the static brane configuration in the background of the self-tuning solution, satisfying (14), is possible only when $`n=3`$. The following fine-tuned value of the $`n`$-brane tension is obtained by substituting the expression for $`r_b`$ in Eq. (33) into Eq. (30) along with Eq. (32):

$$\sigma =\frac{1}{\kappa _D^2}\frac{4(n1)}{\mathrm{\Delta }}r^1h_{}^{\frac{2(\mathrm{\Delta }n+1)}{(n1)\mathrm{\Delta }}}h_+^{\frac{1}{2}}|_{r=r_b},$$
(34)
where $`h_{}(r_b)=\frac{42\mathrm{\Delta }}{4\mathrm{\Delta }}\mathrm{tanh}^2\alpha `$ and $`h_+(r_b)=\frac{\mathrm{\Delta }}{\mathrm{\Delta }4}`$. From this explicit expression for the brane tension, we see that real-valued nonzero brane tension is possible when $`\mathrm{\Delta }>4`$ or $`\mathrm{\Delta }<0`$. So, in order to have the real-valued $`n`$-brane location and tension, the dilaton coupling parameter $`a`$ should be such that $`\mathrm{\Delta }>4`$. When $`\mathrm{\Delta }=4`$, for which the $`n`$-brane is located at the singularity $`r=0`$, the $`n`$-brane tension becomes infinite. When $`\mathrm{\Delta }=0`$, for which the $`n`$-brane is located at the event horizon $`r=r_+`$, the $`n`$-brane tension is also infinite.

Finally, we study a massless bulk scalar field $`\mathrm{\Phi }`$ satisfying the following equation of motion:

$$\frac{1}{\sqrt{G}}_M(\sqrt{G}G^{MN}_N\mathrm{\Phi })=0.$$
(35)
We take the following ansatz for the solution:

$$\mathrm{\Phi }(x^\mu ,z)=u_m(x^\mu )f_m(z),A^2_t^2u_m+B^2\delta ^{ij}_i_ju_m=m^2u_m.$$
(36)
Then, Eq. (35) reduces to the Sturm-Liouville equation:

$$_z[AB^n_z]f_m=m^2AB^nf_m.$$
(37)
In terms of the new $`z`$-dependent function $`\stackrel{~}{f}_m=A^{\frac{1}{2}}B^{\frac{n}{2}}f_m`$, this equation takes the following zero energy eigenvalue Schrödinger equation:

$$\frac{d^2\stackrel{~}{f}_m}{dz^2}+V(z)\stackrel{~}{f}_m=0,$$
(38)
with the potential

$$V(z)=\frac{1}{2}\frac{A^{\prime \prime }}{A}+\frac{n}{2}\frac{B^{\prime \prime }}{B}+\frac{n}{2}\frac{A^{}}{A}\frac{B^{}}{B}\frac{1}{4}\left(\frac{A^{}}{A}\right)^2+\frac{n(n2)}{4}\left(\frac{B^{}}{B}\right)^2+m^2.$$
(39)
This equation can be rewritten in the following convenient form:

$$\left[\overline{Q}Q+m^2\right]\stackrel{~}{f}_m=0,$$
(40)
where

$$Q=\frac{d}{dz}+\frac{d\mathrm{ln}A^{\frac{1}{2}}B^{\frac{n}{2}}}{dz},\overline{Q}=\frac{d}{dz}+\frac{d\mathrm{ln}A^{\frac{1}{2}}B^{\frac{n}{2}}}{dz}.$$
(41)
So, the zero mode ($`m=0`$) annihilated by $`Q`$ is

$$\stackrel{~}{f}_0=A^{\frac{1}{2}}B^{\frac{n}{2}},$$
(42)
from which we see that the Kaluza-Klein zero mode $`\mathrm{\Phi }(x^\mu ,z)=u_0(x^\mu )f_0(z)=u_0(x^\mu )A^{\frac{1}{2}}B^{\frac{n}{2}}\stackrel{~}{f}_0`$ of the bulk scalar field is independent of $`z`$, just as in the symmetrically warped brane world case. It is straightforward to show that such result continues to hold for other bosonic fields. In order for the eigenvalue problem for the Sturm-Liouville equation to be well-defined, appropriate boundary conditions have to be satisfied. For the Sturm-Liouville equation of the form (37), the eigenvalue $`m^2`$ is real and the eigenfunctions $`f_m`$ with different eigenvalues are orthogonal w.r.t. the weighting function $`w(z)=AB^n`$, provided the operator $`=_z[AB^n_z]`$ is self-adjoint. The operator $``$ is self-adjoint, in general when the boundary condition $`AB^n(f_m^{}f_m^{}f_m^{}^{}f_m)|_{z_b}^{z_H}=0`$ is satisfied. There are many ways in which this boundary condition can be satisfied. One of possible ways would be to first write the eigenfunction as $`f_n(z)=F_m(m^2z)`$ and require that $`F_m(m^2z)`$ vanishes at the boundaries. This requirement determines the eigenvalues $`m^2`$ satisfying the required boundary condition.