text:
# Zipf’s Law in Importance of Genes for Cancer Classification Using Microarray Data

## 1 Introduction

Not all of the more than 30,000 human genes and perhaps a multiple of that number of proteins \[Venter, et al., 2001; IHGSC, 2001\] are expected to be useful for all situations. For a specialized phenotype, a particular genetic disease, or a biological process of interest, only a portion of the genes are involved. Even for the involved genes, their contribution to that phenotype (or disease, or process) varies: some genes can be more important than others. The number of genes involved in a phenotype/disease/process also differs from one case to another: some phenotypes require a large number of genes to contribute (polygenic), others may only need a few (oligogenic). This theme has been discussed in the study of human genetic diseases: those caused by mutation in one gene are called Mendalian or simple diseases, whereas those caused by multiple mutations in many genes, or due to an interplay between genetic susceptibility and environmental factors, are called complex genetic diseases (see, e.g., \[Haines, et al., 1998\]).

Once a phenotype, a disease, or a biological process has been chosen, we may examine each gene’s contribution to the phenotype/disease/process under study, and its “importance” can be ranked. This ranking is a very common practice in quantitative analyses of natural and social phenomena. For example, cities can be ranked by their population, companies can be ranked by revenue or profit, ecological disasters can be ranked by the population of species being destroyed, etc. Plotting the measurement by which the ranking is determined, versus the ranking number (1 for the best or the largest, etc., 2 for the second-best, $`\mathrm{}`$) is called a Zipf’s plot, originated from the work of George Zipf \[Zipf, 1935, 1949\].

Among the large number of Zipf’s plots drawn by Zipf, especially that of the frequency of word occurrence in human language, it was observed that they tend to follow a power-law (algebraic) function: $`y1/r^\alpha `$, where $`r`$ is the rank, $`y`$ is the measurement that determines the rank, and $`\alpha 1`$ the scaling exponent. This power-law behavior is called Zipf’s law. Power-law Zipf’s plot where the value of $`\alpha `$ is not close to 1 may be called a generalized form of Zipf’s law. Examples of Zipf’s law, besides those mentioned above, include magnitude of earthquakes \[Sornett, et al., 1996\], popularity of webpage visits \[Crovella & Bestavros, 1997\], usage of library books, frequency of key word search \[Chen & Leimkuhler, 1986\], etc. More references on Zipf’s law can be found at an online resource at http://linkage.rockefeller.edu/wli/zipf.

This paper addresses the question of whether human genes can be ranked by their “importance” in their contribution to a particular phenotype such as cancer. Once a measurement of the importance of a gene is defined, we ask whether this measure of importance as a function of the rank follows an inverse power-law function, i.e. Zipf’s law. There are two aspects of this study. Qualitatively, Zipf’s plot will show whether or not the importance of genes, at least for a particular phenotype classification task, follows a continuous spectrum. A discontinuous gap would separate “important” genes from “unimportant” ones. Quantitatively, how fast the importance of genes falls off as a function of the rank provides information on how many genes should be kept for this phenotypic classification.

The invention of the microarray (DNA chip) was a breakthrough in biotechnology (see, e.g., \[Southern, 1996; Ekins & Chu, 1999\]), which makes it possible to monitor the expression level of thousands of gene products such as mRNA simultaneously. A comparison of these gene expression levels in cancerous as well as normal tissues will point to genes relevant to oncogenesis. When each gene is examined individually, the more its expression is different in cancerous tissues from that in normal tissues, the more relevant this gene is to cancer development, and the more important it is for cancer classification. In this paper, we use microarray data to rank the importance of genes, and Zipf’s plot of this measure of importance is studied. Section 2 introduces the mathematical formula and quantities used in this paper; section 3 illustrates the likelihood-rank plots (Zipf’s plot) for the simplest cancer classifications (two-type); section 4 shows the Zipf’s plot for multiple-type classifications; section 5 shows the Zipf’s plot for paired samples classification; section 6 discusses the issue concerning the exponent in the power-law plot; and section 7 contains other discussions.

## 2 Measure of importance of genes in cancer classification: likelihood

The measurement we use for the importance of genes in cancer classification is the maximized likelihood, which is proportional to the probability of observing the data when a model is given, and when the parameters in the model are adjusted to give the maximum value. Mathematically, likelihood $`L`$ is \[Edwards, 1972\]:

$$L=cP(D|M,\theta ),$$
(1)
where $`D`$ is all data points in a data set, $`M`$ is a model, $`\theta `$ represents all parameters in the model, and $`c`$ is a proportional constant (often set to be 1). The maximized likelihood is

$$\widehat{L}=\underset{\theta }{\mathrm{max}}L=cP(D|M,\widehat{\theta }),$$
(2)
where the $`\widehat{}`$ represents a maximization/estimation procedure. The microarray data sets have the form of $`\{x_{1i},x_{2i},\mathrm{},x_{pi},y_i\}`$ ($`i=1,2,\mathrm{}N`$). Each sample point $`i`$ (one microarray experiment, one tissue sample) contains measurements of (logarithm of ) mRNA expression level of $`p`$ genes ($`p`$ is typical of the order of thousands), and $`y`$ is a categorical label indicating a condition (e.g. cancerous or normal). The raw data in a microarray experiment could be more complicated: one has to consider background noise, normalization, and controls. These considerations depend on the type of chips: Stanford/cDNA array with two fluorescence dyes \[Schena, et al., 1995; Shalon et al., 1996\] or Affymetrix/oligonucleotide arrays with only one image intensity to measure but multiple oligonucleotide probes \[Fodor, et al., 1991\]. In this paper, only the filtered/processed data are used and the subtle issue of scaling/normalization of the raw data is not discussed.

The model $`M`$ is a classification model (classifier, predictor, discriminator, supervised learning machine), that classifies the label $`y`$ by the gene expression level $`\{x_j\}`$ (logarithm of an image intensity from the DNA chip). We use a particularly simple classifier, the single-gene logistic regression:

$$P(y_i=1|x_{ji})=\frac{1}{1+e^{a_jb_jx_{ji}}}.j=1,2,\mathrm{}p\text{, and }i=1,2,\mathrm{}N.\text{ }$$
(3)
In other words, the probability of a sample being in one class is a “sigmoid” or “logistic” function of the (log) expression level. If the coefficient $`b`$ is positive, larger expression levels lead to higher probabilities of being the $`y=1`$ label; if $`b<0`$, larger expression level corresponds to the $`y=0`$ label.

The likelihood of the whole data set is the product of these model-based probabilities (for a given gene):

$$L_j=\underset{i=1}{\overset{N}{}}\left[P(y=1|x_{ji})𝐈_{y_i=1}+(1P(y=1|x_{ji}))𝐈_{y_i=0}\right],j=1,2,\mathrm{}p\text{ }$$
(4)
where $`𝐈`$ is the indicator function (1 or 0 depending on whether the condition is true or not). Eq.(4) is maximized by adjusting the parameters in the model. The maximized likelihood $`\widehat{L}`$ can then be used to rank all genes: the larger the $`\widehat{L}`$, the higher the ranking.

## 3 Zipf’s plot of classification likelihood

The binomial (two-class) logistic regression Eq.(3) has been applied to two microarray data sets. The first is colon cancer data from Princeton University \[Alon, et al., 1999\]. The expression levels of 2000 genes were available for 62 tissue samples: 40 cancerous and 22 normal. The second data set is leukemia data from Whitehead Institute/MIT \[Golub, et al., 1999\]. Expression levels of 7129 genes were measured on 72 leukemia samples: 47 obtained from tissues of one subtype of leukemia, acute lymphoblastic leukemia (ALL), and 25 from tissues of another subtype of leukemia, acute myeloid leukemia (AML). Although the number of labels is 2 for both datasets, we distinguish cancerous from normal tissues in the first data set, whereas distinguish one cancer subtype from another in the second data set.

For the colon cancer data, the single-gene maximized likelihood Eq.(4) for each gene was calculated and ranked. The likelihood-rank plot (Zipf’s plot) is shown in Fig.1 in log-log scale. The power-law behavior of the curve is clearly visible. Fitting the first 600 genes using a generalized form of Zipf’s law:

$$\widehat{L}_r1/r^\alpha ,$$
(5)
leads to exponent $`\alpha 2.18`$, whereas the fitting of the top 1000 genes leads to $`\alpha 2.10`$. The genes from roughly rank 1000 to 2000 do not seem to follow the same power-law decay of the likelihood. As will be discussed more later in this paper, the exponent $`\alpha `$ is not an intrinsic quantity for our likelihood-rank plots. The reason is that the likelihood is a product of probabilities of $`N`$ sample points $`L_r(p_1p_2\mathrm{}p_N)_r\overline{p}_r^N`$; if the per-sample averaged likelihood $`\overline{p}_r`$ of rank-r gene does not change with the sample size $`N`$, the exponent $`\alpha `$ is then a function of $`N`$: $`\alpha N\mathrm{log}(\overline{p}_r)/log(r)`$.

Fig.2 shows the Zipf’s plot of the second data set. Among the 72 samples, 38 are from bone marrow tissues, and were separated as a training set in \[Golub, et al., 1999\]. These 38 samples were considered to be more homogeneous, while the rest of the samples were from various sources or other tissue types such as peripheral blood \[Golub, et al., 1999\], and may not be homogeneous. From Fig.2, the Zipf’s plot obtained from the training set seems to follow a generalized form of Zipf’s law with a fitting exponent of $`\alpha 2.56`$ from the top 900 genes. The Zipf’s plot of all sample points seems to deviate from a power-law trend around rank 100-200. As mentioned earlier, the exponent $`\alpha `$ from a bigger data set is indeed larger than the exponent from a small data set.

When the top genes are examined, it was found that the top-ranking genes for the training set \[Li & Yang, 2001\] and those for the overall set \[Li, et al., 2001\] may not be identical. For example, the top performing gene for the 38-sample training set is No.4847, zyxin, a gene encoding the LIM domain protein used in cell adhesion in fibroblasts \[Golub, et al., 1999\]. On the other hand, the best performing gene for the 72-sample combined set is No. 1834, CD33 antigene which encodes cell surface proteins commonly found in AML leukemia cells (see, e.g., \[Lauria, et al., 1994\]). The zyxin becomes the rank-4 gene for the combined dataset. This difference reflects a certain degree of inhomogeneity between samples in the training set and those in the remaining set (testing or validation set). Note from Fig.1 that for the training set, genes from rank 3 to 9 exhibits similar likelihood, and form a flat step on the Zipf’s plot. Such steps are a dominant feature in the Zipf’s plot of the frequency of word occurrence in randomly generated texts \[Li, 1992\].

## 4 Classification of multiple cancer classes

The logistic regression Eq.(3) for a two-class data set can be generalized to multiple classes: multinomial logistic regression \[Agresti, 1996\]:

$$P(y_i=I|x_{ji})=\frac{e^{a_Ib_Ix_{ji}}}{_{K=1}^Ce^{a_Kb_Kx_{ji}}}j=1,2,\mathrm{}p\text{}I=1,2,\mathrm{}C\text{}i=1,2,\mathrm{}N\text{,}$$
(6)
where the label $`y`$ can be in one of the $`C`$ classes, and there are two parameters for each class (though only $`C1`$ class probabilities are independent). A gene is important (higher maximized likelihood) if it is more able to distinguish all $`C`$ classes simultaneously.

For this multiple-class analysis, we use the microarray data for lymphoma (Stanford University and National Cancer Institute \[Alizadeh, et al., 2000\]). There are a total 96 tissue samples, with 66 cancerous and 30 normal. Within the cancer samples, there are 46 diffuse large B-cell lymphoma (DLBCL), 9 follicular lymphoma (FL), 11 chronic lymphocyte leukemia (CLL). These 3 cancer subtypes plus the normal are the 4 classes to be distinguished. In \[Alizadeh, et al., 2000\], it is also recommended that two new subclasses of DLBCL can be defined based on the microarray data. These are the germinal centre B-like DLBCL (GC-DLBCL) and the activated B-like DLBCL (A-DLBCL). The GC-DLBCL, A-DLBCL subclasses plus the normal can be the 3 classes to be distinguished.

Fig.3 shows the Zipf’s plot of multinomial logistic regression likelihood for both the 4-class and the 3-class classification. The Zipf’s plot for the 3-class multinomial logistic regression follows a perfect power-law function for all genes. No deviation from power-law behavior was observed even for the low-ranking genes. The Zipf’s plot for the 4-class multinomial logistic regression, on the other hand, does not seem to follow a power-law function. The fall off near the rank-200 gene is perhaps due to the computer roundoff error since the value of likelihood at rank-200 is already as low as $`10^{39}`$ (the likelihood is multipled by $`10^{10}`$ in Fig.3). The fitting of the top 200 genes by a power-law function is, however, very good.

## 5 Cancer treatment effect

Another interesting situation is provided by the data set from Stanford University and the Norwegian Radium Hospital \[Perou, et al., 2000\]. Part of the data is the expression level of 8102 genes measured before and after a 16-week course of doxorubicin chemotherapy \[Perou, et al., 2000\] on 20 patients. Naively, one may consider these 40 microarray experiments as an example of binary logistic regression. However, logistic regression in Eq.(3) does not apply because it requires samples to be independent of each other, whereas microarray experiments done on the same patient are clearly not independent. This situation can be handled by the logistic regression of paired case-controls (case is a sample with label 1, and control is a sample with label 0) \[Breslow & Day, 1980\]:

$$P(x(case)_{ji}x(control)_{ji})=\frac{1}{1+e^{b_j(x(case)_{ji}x(control)_{ji})}},j=1,2,\mathrm{}p\text{}i=1,2,\mathrm{}N\text{}$$
(7)
Note that the probability is no longer that for observing $`y`$ – when a case sample and a control sample are paired, their $`y`$ value is fixed at 1 and 0 – but that of observing the difference of $`x`$’s. Also note that the first parameter $`a`$ is now zero.

Fig.4 shows the Zipf’s plot for genes in the breast cancer microarray data. The top three genes were able to perfectly identify the chemotherapy effect – the expression level is always higher (or lower) before the treatment than after in all 20 samples. The likelihood of a perfect fit is equal to 1. Unlike the previous three plots, the likelihood of genes in Fig.4 does not follow a power-law function, or even a smooth function. There seems to be a gap from the ranking of 19 to the ranking of 20-25. Fitting the two segments (from 4 to 19, and from 25 to 500) by power-law functions leads to different exponents (2.51 versus 1.77).

Plots like Fig.4 can be good news for a microarray data analysis, because there seems to be a separation between “relevant” and “irrelevant” genes. Irrelevant genes are not important for distinguishing samples before and after chemotherapy, and may be discarded for further analysis. Of course, this is only a rough description of the gene set. A more systematic approach can be based on the framework of model selection (see, e.g., \[Burnham & Anderson, 1998\]) or model averaging (see, e.g., \[Geisser, 1993\]). With model selection, the expression level of different genes can be combined, and adding one gene increases the number of parameters in the model by 1. The appropriate number of genes to be included in a classification is the model with the best “model selection criterion” such as Akaike information criterion \[Burnham & Anderson, 1998\] or Bayesian information criterion \[Raftery, 1995\], with a best balance between a larger likelihood value and fewer numbers of parameters \[Li & Yang, 2001, Li, et al., 2001\]. With model averaging, there is in principle no limitation on the number of genes to be included, but irrelevant genes have smaller weights in an averaged classifier \[Golub, et al., 1999, Li & Yang, 2001\]. This makes the effective number of genes used much smaller than the apparent number.

## 6 Scaling exponent

As mentioned earlier, the exponent of the inverse power-law fitting function for Zipf’s plot (Figs. 1-4) depends on the sample size. The reason for this is that the measurement for the importance of a gene, the sample classification likelihood, is a function of the sample size. Here we ask the question whether one can define a normalized exponent. Since $`\widehat{L}_r1/r^\alpha `$, we can obviously draw Zipf’s plot for the per-sample likelihood: $`\overline{p}_r=\widehat{L}_r^{1/N}1/r^{(\alpha /N)}`$. There is another consideration for classification of more than two labels: it is unfair to compare per-sample classification likelihoods when the number of classes differs. Just by random guess, the probability of classifying the label correctly for a binary-class case is 0.5, whereas that for classifying $`C`$ labels is $`1/C`$. We can normalize the per-sample likelihood $`\overline{p}_r`$ by the random-guess probability: $`\widehat{l}_r=\overline{p}_r/(1/C)C/r^{(\alpha /N)}`$ with the same scaling exponent.

Fig.5 shows the Zipf’s plot of normalized per-sample likelihood $`\widehat{l}_r`$ for all data sets analyzed so far. These curves can be compared in two ways. First, the performance of the top-ranking genes can be compared by their $`\widehat{l}_r`$’s. It is clear from Fig.5 that for the leukemia data, the performance obtained from the training set is better than that obtained from the whole set, indicating a possible heterogeneity in the data set. The performance of the top genes for the breast cancer data set is better (in terms of classification) than that for the leukemia data set, which in turn is better than that for the colon cancer data set.

It is usually difficult to compare performance between a two-label classification and a multiple-label classification, because it depends on the base-line expectation. Two base-line expectations (called null models) were defined in \[Li & Yang, 2001, Li, et al., 2001\]: one is to randomly guess all classes with equal probability, and another is to guess the class by the proportion of samples with this class in the data set. In Fig.5, the first base-line expectation is used as the normalization factor. The performance of the 4-label classification relative to this base-line expectation for the lymphoma data is considered to be better than that of the 3-label classification, partly due to its low expectation.

Different data sets in Fig.5 can also be compared by the rate of falling of likelihood. The fall-off rate, as measured by the scaling exponent $`\alpha `$, ranges from 3 to 9 per 100 samples (i.e., $`\alpha /N0.030.09`$). It is interesting that scaling exponents obtained from the leukemia data (both the whole data set and the training set) and from the lymphoma data (two DLBCL subtypes plus normal) are almost identical (around 65 - 67 per 100 samples). The breast cancer data is different from other data sets for falling faster in the likelihood-rank plot.

## 7 Discussion

Results from Figs. 1-4 established that when microarray data is used to classify cancer tissues, the classification likelihood of individual genes follows a generalized form of Zipf’s law. The reason that these power-law functions in the Zipf’s plot do not have a scaling exponent equal to 1 is partly because exponent $`\alpha `$ in Eq.(5) depends on the sample size $`N`$. Besides the issue of the exponent value, the overall power-law trend is excellent: a perfect power-law function for all points in Fig.3 (3-label classification, in 3.5 decades), a partial power-law fitting in a range of 2.5 decades (Fig.1), 3 decades (Fig.2), and 2 decades (Fig.3, 4-label classification), respectively. The only poor fitting by a power-law function is Fig.4. The Zipf’s plot is flattern out for low-ranking genes in Figs. 1 and 2, while it drops off in Fig.3 (4-label classification). Which functional form is more generic for low-ranking genes is not clear, perhaps our results are sample-size sensitive.

Our Zipf’s plot can be compared to those obtained from biomolecular sequences. In \[Gamow & Yc̆as, 1955\], a Zipf’s plot for 20 amino acids usage in protein sequences was presented (in linear-linear scale). Redrawing their data in log-log scale does not show any power-law behavior (result not shown). Recently, it was claimed that oligonucleotide frequencies in DNA sequences follow Zipf’s law \[Mantegna, et al., 1994\]. This paper drew criticism on its strong claim concerning the connection between Zipf’s law and human language \[Martindale & Konopka, 1995; Israeloff, et al., 1996, Bonhoeffer, et al., 1996a,1996b, Voss, 1996\]: one of the best counter-examples is Zipf’s law in money-typing texts \[Li, 1992\]. Also, the paper did not show convincingly that the Zipf’s plot for oligonucleotide usage was better fitted by a power-law function \[Martindale & Konopka, 1996\]: the deviation from the power-law fitting function can be gradual and systematic, an indication that the power-law function is not the best choice of fitting function. Finally, the scaling exponent in the power-law fitting function in \[Mantegna, et al., 1994\] is much smaller than 1 \[Li, 1996\]. In comparison, our Zipf’s plots in Figs 1-5 are a much better example of a generalized form of Zipf’s law than those in \[Gamow & Yc̆as, 1955\] and \[Mantegna, et al., 1994\].

One may wonder whether the power-law behavior in Figs 1-5 can be derived by a simple random model. In \[Gamow & Yc̆as, 1955\], the Zipf’s plot of the frequency of amino acids usage was compared to a “random partition of a unit length” model. Suppose a unit interval is randomly partitioned into $`p`$ segments (e.g. $`p=20`$ for 20 amino acids). These segments are ranked by their size $`L_{(1)}=max_i(L_i),L_{(2)}=max_i(L_iLL_{(1)})L_{(p)}=min_i(L_i)`$, etc. If this random partition is repeated, the mean value of the ranked size can be shown to be \[Gamow & Yc̆as, 1955\]:

$$<L_{(r)}>=\frac{1}{p}\underset{i=1}{\overset{p+1r}{}}\frac{1}{p+1i}\text{ }r=1,2,\mathrm{}p$$
(8)
Drawing $`<L_{(r)}>`$ from Eq.(8) versus rank $`r`$ shows that it is a straight line in linear(y)-log(x) plot, and not a straight line in either log-log or log(y)-linear(x) plots (results not shown). Actually this analytic result may not be applicable to our case, because we keep track of the performance of a given gene on all samples, whereas in Eq.(8) the longest interval (or any given ranking interval) is averaged over all random simulation. In any case, the power-law behavior in Fig.1-5 does not seem to be explainable by a simple random model.

We may ask whether Zipf’s plot has any practical application for microarray data analysis. In information retrieval \[van Rijsbergen, 1975, Salton, 1988\] and library/documentation science \[Egghe & Rousseau, 1990\], Zipf’s law is an important foundation that many applications are based upon. It is one of the “bibliometric laws” \[White & McCain, 1989\] concerning regularities in bibliographies, lists of authors, citation lists, etc. For the purpose of finding relevant, content-bearing words (“keywords”), common (highest-ranking) and rare (lowest-ranking) words should be avoided \[Luhn, 1957,1958\].

Do we have a similar situation where the highest-ranking genes may not be interesting for cancer classification? (Lowest-ranking genes are obviously not interesting.) Our ranking system is not really the same as for word usage, since a discrimination or classification ability has already been included in our definition, whereas it is not included in the word usage example. In this sense, our top ranking genes are in fact most efficient for the purpose of cancer classification. On the other hand, if one is more interested in subtle gene effects, not in the known main/dominant effect, it is perhaps useful to remove the well-known genes from the list and examine other genes in a future study. This idea was not tried in our previous analysis \[Li & Yang, 2001, Li, et al., 2001\].

In conclusion, a generalized form of Zipf’s law was observed in microarray data for the likelihood on cancer classification using single-gene logistic regression. We suspect that this power-law behavior is generic rather than an exception. A rank-likelihood plot (Zipf’s plot) can be a useful quantitative tool for discriminant microarray data analysis.

Acknowledgements: The author would like to thank Yaning Yang for discussion on the random partition of unit interval, Ronald Rousseau and Heting Chu for suggesting references on library sciences, Victoria Haghighi and Joanne Edington for help with the microarray data. The work was supported by NIH grants K01HG00024 and HG00008.

## References

Agresti A (1996), An Introduction to Categorical Data Analysis (Wiley & Sons).

Alizadeh AA, Eisen MB, et al. (2000), “Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling”, Nature, 403:503-511.

Alon U, Barkai N, Notterman DA, Gish K, Ybarra S, Mack D, Levine AJ (1999), “Broad patterns of gene expression revealed by clustering analysis of tumor and normal colon tissues probed by oligonucleotide arrays”, Proceedings of National Academy of Sciences, 96(12):6745-6750.

Bonhoeffer S, Herz AVM, Boerlijst MC, Nee S, Nowak MA, May RM (1996a) “Explaining ’linguistic features’ of noncoding DNA”, Science, 271(5245):14-15.

Bonhoeffer S, Herz AVM, Boerlijst MC, Nee S, Nowak MA, May RM (1996b), “No signs of hidden language in noncoding DNA” (letters), Physical Review Letters, 76(11):1977.

Breslow NE, Day NE (1980), Statistical Methods in Cancer Research. I - The Analysis of Case-Control Studies (International Agency for Research on Cancer, Lyon).

Burnham KP, Anderson DR (1998), Model Selection and Inference (Springer).

Chatzidimitriou-Dreismann CA, Streffer RMF, Larhammar D (1996), “Lack of biological significance in the ’linguistic features’ of noncoding DNA-a quantitative analysis”, Nucleic Acids Research, 24(9):1676-1681.

Chen YS, Leimkuhler FF (1987), “Analysis of Zipf’s law: an index approach”, Information Processing and Management, 23:71-182.

Crovella ME, Bestavros A (1997), “Self-similarity in world wide web traffic: evidence and possible causes”, IEEE/ACM Transactions on Networking, 5(6):835-846.

Edwards AWF (1972), Likelihood (Cambridge Univ Press).

Egghe L, Rousseau R (1990), Introduction to Informetrics: Quantitative Methods in Library, Documentation and Information Science (Elsevier).

Ekins R, Chu FW (1999), “Microarrays: their origin and applications”, Trends in Biotechmology, 17:217-218.

Fodor SP, Read JL, Pirrung MC, Stryer L, Lu AT, Solas D (1991), “Light-directed, spatially addressable parallel chemical synthesisLight-directed, spatially addressable parallel chemical synthesis”, Science, 251:767-773.

Gamow G, Yc̆as M (1955), “Statistical correlation of protein and ribonucleic acid composition”, Proceedings of National Academy of Sciences, 41(12):1011-1019.

Geisser S (1993), Predictive Inference: An Introduction (Chapman & Hall).

Golub TR, Sonim DK, Tamayo P, Huard C, Gassenbeek M, Mesirov JP, Coller H, Loh ML, Downing JR, Caligiuri MA, Bloomfield CD, Lander ES (1999), “Molecular classification of cancer: class discovery and class prediction by gene expression monitoring”, Science, 286:531-536.

IHGSC (International Human Genome Sequencing Consortium) (2001), “Initial sequencing and analysis of the human genome”, Nature, 409:860-921.

Israeloff NE, Kagalenko M, Chan K (1996), “Can Zipf distinguish language from noise in noncoding DNA?” (letters) Physical Review Letters, 76(11):1976.

Haines JL, Pericak-Vance MA (1998), eds. Approaches to Gene Mapping in Complex Human Diseases (Wiley-Liss).

Lauria F, Raspadori D, et al. (1994), “Increased expression of myeloid antigen markers in adult acute lymphoblastic leukemia patients: diagnostic and prognostic implications”, British Journal of Haematology, 87:286-292.

Li W (1992), “Random texts exhibit Zipf’s-law-like word frequency distribution”, IEEE Transactions on Information Theory, 38:1842-1845.

Li W (1996), “Comments on ’Bell curves and monkey languages’ ” (letters to the editor), Complexity, 1(6):6.

Li W, Yang Y (2001), “How many genes are needed for a discriminant microarray data analysis?”, Proceedings of the Critical Assessment of Microarray Data Analysis Workshop (CAMDA 2000).

Li W, Yang Y, Edington J, Haghighi F (2001), “Determining the number of genes needed for cancer classification using microarray data”, submitted.

Luhn HP (1957), “A statistical approach to mechanized encoding and search of literature information”, IBM Journal of Research and Development, 1:309-317.

Luhn HP (1958), “The automatic creation of literature abstract”, IBM Journal of Research and Development, 2:159-165.

Mantegna RN, Buldyrev SV, Goldberger AL, Havlin S, Peng CK, Simon M, Stanley HE (1994), “Linguistic features of noncoding DNA sequences”, Physical Review Letters, 73:3169-3172.

Martindale C, Konopka AK (1995), “Noncoding DNA, Zipf’s law, and language” (letters), Science, 268(5212):789.

Martindale C, Konopka AK (1996), “Oligonucleotide frequencies in DNA follow a Yule distribution”, Computers and Chemistry, 20:35-38.

Perou CM, Sorlie T, et al. (2000), “Molecular portraits of human breast tumors”, Nature, 406:747-752.

Raftery AE (1995), “Bayesian model selection in social research”, in Sociological Methodology, ed. PV Marsden (Blackwells), pages 185-195.

Salton G (1988), Automatic Text Processing: The Transformation, Analysis, and Retrieval of Information (Addison-Wesley).

Schena M, Shalon D, Davis RW, Brown PO (1995), “Quantitative monitoring of gene expression patterns with a complementary DNA microarray”, Science, 270:467-470.

Shalon D, Smith SJ, Brown PO (1996), “A DNA microarray system for analyzing complex DNA samples using two-color fluorescent probe hybridization”, Genome Research, 6:639-645.

Sornette D, Knopoff L, Kagan YY, Vanneste C (1996), “Rank-ordering statistics of extreme events: application to the distribution of large earthquakes”, Journal of Geophysical Research, 101(B6):13883-13894.

Southern EM (1996), “DNA chips: analysing sequence by hybridization to oligonucleotides on a large scale”, Trends in Genetics, 12:110-115.

van Rijsbergen CJ (1975), Information Retrieval (Butterworths).

Venter JC (2001), et al. “The sequence of the human genome”, Science, 291:1304-1351.

Voss RF (1996), “Linguistic features of noncoding DNA sequences - Comment” (letters), Physical Review Letters, 76(11):1978.

White HD, McCain KW (1989), “Bibliometrics”, Annual Review of Information Science and Technology, 24:119-186.

Zipf GF (1935), Psycho-Biology of Languages (Houghton-Mifflin)

Zipf GF (1949), Human Behavior and the Principle of Least Effort (Addison-Wesley).


format: text

text:
# I Introduction

## I Introduction

Understanding why the vacuum energy density is essentially zero is a fundamental challenge of contemporary physics. The vacuum energy density is a free parameter of nature unprotected from large quantum corrections and, therefore, an exquisitely small value for this parameter is an unexplained fine-tuning.

One approach to dealing with this problem  relies on braneworld theories with infinite volume extra dimensions. The first attempt of such models appeared in ; however, we consider theories based on a different approach . These theories exhibit the following unique properties. First, gravity becomes higher-dimensional at large distances; and second, the infinite volume extra dimensions allow for exact bulk supersymmetry (compatible with SUSY broken on the brane) which can control the value of the bulk cosmological term. As a result, it is natural to expect that the large scale cosmology is that of the higher-dimensional theory with a cosmological term that can be naturally zero due to unbroken bulk SUSY. Even so, there are questions that must be addressed before such a solution can be regarded as viable:

* What is the mechanism that makes gravity appear 4-dimensional on the brane?
* Can the Universe evolve to the desired (present) state without sacrificing successful features of standard Friedmann–Lemaître–Robertson–Walker (FLRW) cosmology during the earlier stages of the evolution?

To provide an answer to the first question we work in the framework of models with an induced intrinsic curvature term on the brane . In this framework gravity is 4-dimensional up to astronomical distances due to the large graviton kinetic term on the brane. Explicit cosmological solutions confirming large distance crossover behavior within this class of theories were already found in.

In the present paper we discuss a particular class of those cosmological scenarios where the cosmological constant on the brane is negative. Irrespective of the magnitude of this constant, the 4-dimensional metric on the brane automatically evolves to a static Minkowski spacetime for generic initial conditions with positive energy density. This attractive scenario does not requires fine-tuning; however, it cannot be regarded as a solution of the cosmological constant problem since phenomenological considerations require a small brane cosmological constant. Nevertheless, the example is rather remarkable since the meaning of cosmological constant problem is changed. Unlike the usual 4-dimensional scenario where the smallness of the vacuum energy is required for the observed flatness of the Universe, in the present framework the Universe is automatically flat for arbitrarily large negative vacuum energy. The constraint on the brane cosmological constant comes from completely different considerations, such as ultralarge distance gravity measurements and 4-dimensional FLRW–cosmological history. In other words, in our scenario we ask not: “Why is the Universe flat?” but rather: “Why does the Universe have FLRW–type history?”

Finally let us note that a byproduct of our scenario is small Lorentz-violating effects,*** The metric on the brane is Minkowskian, even though the energy momentum tensor is nonzero and does not respect Lorentz symmetry. The metric of the 5-dimensional space, which is supported solely by the energy-momentum on the braneworld, also does not respect 4-dimensional Lorentz symmetry. Any slice of the 5-dimensional spacetime parallel to the braneworld is indeed Minkowski, however, the speed of light varies from slice to slice. Gravity can probe between slices and is capable of explicitly reflecting this 4-dimensional Lorentz symmetry violation. for which there is a growing interest in both particle phenomenology contexts (see for example ), as well as in braneworld scenarios .

We first discuss the ingredients necessary to construct our model, as well as describe its cosmological evolution and the global 5-dimensional structure of spacetime. We then go on to evaluate the propagation of perturbations on the braneworld itself and phenomenological constraints. We discuss these limitations in detail in the concluding remarks.

## II The Solution, Dynamics and Cosmology

Consider a three-brane embedded in a 5-dimensional spacetime. The bulk is empty; all energy-momentum is isolated on the brane. The action is

$$S_{(5)}=\frac{1}{2}M_{(5)}^3d^5x\sqrt{\stackrel{~}{g}}\stackrel{~}{R}+d^4x\sqrt{g}_m.$$
(1)
The first term in Eq. (1) corresponds to the Einstein-Hilbert action in five dimensions for a 5-dimensional metric $`\stackrel{~}{g}_{AB}`$ (bulk metric) with Ricci scalar $`\stackrel{~}{R}`$. This action is sufficient to generate our flattening cosmological solution. However, in order to obtain 4-dimensional gravitation on the braneworld as well as 4-dimensional cosmology, we also consider an intrinsic curvature term which is generally induced by radiative corrections by the matter density on the brane :

$$\frac{1}{2}M_{(4)}^2d^4x\sqrt{g}R.$$
(2)
Similarly, Eq. (2) is the Einstein-Hilbert action for the induced metric $`g_{cd}`$ on the brane, $`R`$ being its scalar curvature. The induced metricThroughout this article, we adopt the following convention for indices: upper case Latin letters $`A,B,\mathrm{}`$ denote 5D indices: $`0,1,2,3,5`$; lower case Latin letters from the beginning of the alphabet: $`c,d,\mathrm{}`$ denote 4-dimensional indices parallel to the brane, lower case Latin letters from the middle of the alphabet: $`i,j,\mathrm{}`$ denote space-like 3D indices parallel to the brane. $`g_{cd}`$ is defined as usual from the bulk metric $`\stackrel{~}{g}_{AB}`$ by

$$g_{cd}=_cX^A_dX^B\stackrel{~}{g}_{AB},$$
(3)
where $`X^A(x^c)`$ represents the coordinates of an event on the brane labeled by $`x^c`$. We neglect higher-derivative terms in the bulk and worldvolume actions as they lead to the modification of gravity at undetectably short distances.

For a brane embedded in a Minkowski spacetime with such an action, it is shown in Ref.  that the usual 4-dimensional Newton’s law for static point like sources on the brane is recovered at observable distances. At ultralarge cosmological distances the gravitational force is given by the 5-dimensional $`1/r^3`$ force law. The crossover length scale between the two different regimes is given by

$$r_0=\frac{M_{(4)}^2}{2M_{(5)}^3}.$$
(4)
We consider a system where there is a negative cosmological constant as well as some arbitrary matter localized on the brane. Another mass scale that arises in our discussion is

$$\mu ^2=\frac{\lambda }{M_{(4)}^2},$$
(5)
where $`\lambda `$ is the (negative) cosmological constant on the brane.

There exists a static solution for this system with a metric determined by the line element

$$ds^2=\left(1+c|y|\right)^2d\tau ^2+dx_idx^i+dy^2,$$
(6)
where the total energy momentum on the brane itself is given by $`\mathrm{diag}(0,p_b,p_b,p_b)`$ and $`c=p_b/2M_{(5)}^3`$. The value of the pressure as a function of the cosmological constant, $`\lambda `$, is determined by the equation of state of the matter component on the brane. The energy density of that matter component is exactly balanced by the negative energy density of the cosmological constant. Note that the metric of the full spacetime explicitly breaks 4-dimensional Lorentz invariance, though if one is confined to the braneworld, the spacetime appears Minkowskian.

Although this static solution may appear to be fine-tuned, we will show that for generic initial states with net positive energy density, the system dynamically asymptotes to this static solution. The general time-dependent line element under consideration is of the form

$$ds^2=N^2(\tau ,y)d\tau ^2+A^2(\tau ,y)\gamma _{ij}dx^idx^j+B^2(\tau ,y)dy^2,$$
(7)
and, the metric components are given by

$`N`$ $`=`$ $`1|y|\ddot{a}\left(\dot{a}^2+k\right)^{1/2}`$ (8)
$`A`$ $`=`$ $`a|y|\left(\dot{a}^2+k\right)^{1/2}`$ (9)
$`B`$ $`=`$ $`1,`$ (10)
where $`k=1,0`$ or $`1`$ is the intrinsic spatial curvature parameter. We take the total energy-momentum tensor which includes matter and the cosmological constant on the brane to be

$$T_B^A|_{_{\mathrm{brane}}}=\delta (y)\mathrm{diag}(\rho _b,p_b,p_b,p_b,0).$$
(11)
When the matter content on the brane is specified, the induced scale factor $`aA(\tau ,y=0)`$ is determined by the Friedmann equations.

### A Friedmann equations without intrinsic curvature

For simplicity, let us ignore for the moment the intrinsic curvature term Eq. (2) in the action. The Friedmann equations derived in  read

$`H^2+{\displaystyle \frac{k}{a^2}}`$ $`=`$ $`{\displaystyle \frac{1}{36M_{(5)}^6}}\rho _b^2`$ (12)
$`{\displaystyle \frac{\ddot{a}}{a}}`$ $`=`$ $`{\displaystyle \frac{1}{36M_{(5)}^6}}\rho _b\left(2\rho _b+3p_b\right),`$ (13)
where $`H`$, the Hubble parameter on the brane, is defined by the usual expression $`H=\frac{\dot{a}}{a}`$. Energy-momentum conservation leads to

$$\dot{\rho }_b+3\frac{\dot{a}}{a}(\rho _b+p_b)=0.$$
(14)
Using Eqs. (1213), the metric components in Eqs. (9) can be written as

$`N`$ $`=`$ $`1+{\displaystyle \frac{1}{6M_{(5)}^3}}|y|\left(2\rho _b+3p_b\right)`$ (15)
$`A`$ $`=`$ $`a\left(1{\displaystyle \frac{1}{6M_{(5)}^3}}|y|\rho _b\right)`$ (16)
$`B`$ $`=`$ $`1.`$ (17)
We may now specify the matter content of our model. Consider a spatially flat brane ($`k=0`$), and a brane energy-momentum tensor given by the sum of that of a cosmological constant $`\lambda `$ and a matter energy momentum tensor, such that

$$\rho _b=\lambda +\rho ,p_b=\lambda +p,$$
(18)
where $`\rho `$ is the energy density of matter. We assume that matter obeys the usual equation of state of the form $`p=w\rho `$. In this case the brane energy density is given by

$$\rho _b=\frac{\rho _0}{a^q}+\lambda ,$$
(19)
where $`q=3(1+w)`$ and $`\rho _0`$ is the initial state matter density. Implicitly, we have chosen $`a(t=t_0)=1`$. The Friedmann equations Eqs. (1213) demand that the scale factor (for $`\lambda 0`$) takes the form

$$\rho _0+\lambda a^q=\left(\rho _0+\lambda \right)\mathrm{exp}\left[\lambda q(tt_0)/(6M_{(5)}^3)\right].$$
(20)
Consider the case when the cosmological constant is negative. Initially, if the energy density is much larger than the magnitude of the cosmological constant, the evolution proceeds in the 5-dimensional analog of the big bang scenario: the scale factor increases as a power of time, the energy density decreases as an inverse-power with respect to time. However, as the energy density decreases, it crosses the threshold of order $`|\lambda |`$. The time-dependence of the scale factor changes. When the net energy density is much smaller than the scale of the cosmological constant, the energy density asymptotes exponentially to zero. Similarly, the scale factor asymptotes exponentially to a constant value. The brane pressure at late time is simply given by $`|\lambda |(1+w)`$ so that one sees from Eq. (16) that the metric asymptotes to the solution given in Eq. (6). The behavior described is generic. It is independent of the initial state, so long as the net energy density of that state is positive.

### B Friedmann equations including intrinsic curvature

The dynamical behavior just described continues to holds true even in the presence of an intrinsic curvature term on the brane. This term is needed to recover 4-dimensional gravity and cosmology on the brane. The equations governing the dynamics are modified simply by replacing the brane energy density $`\rho _b`$ (respectively pressure $`p_b`$) by the sum of the brane energy density $`\rho _b`$ (pressure $`p_b`$) and the effective energy density $`\rho _{curv}`$ (effective pressure $`p_{curv}`$) due to the presence of the brane intrinsic curvature term. The quantities $`\rho _{curv}`$ and $`p_{curv}`$ are given by

$`\rho _{curv}`$ $`=`$ $`3M_{(4)}^2\left\{{\displaystyle \frac{\dot{a}^2}{a^2}}+{\displaystyle \frac{k}{a^2}}\right\}`$ (21)
$`p_{curv}`$ $`=`$ $`M_{(4)}^2\left\{{\displaystyle \frac{\dot{a}^2}{a^2}}+2{\displaystyle \frac{\ddot{a}}{a}}+{\displaystyle \frac{k}{a^2}}\right\}.`$ (22)
For example Eq. (12) now reads

$$\sqrt{H^2+\frac{k}{a^2}}=\frac{1}{6M_{(5)}^3}\rho _br_0\left(H^2+\frac{k}{a^2}\right),$$
(23)
which can, in turn, be rewritten

$$\sqrt{H^2+\frac{k}{a^2}}=\frac{1}{2r_0}\left[1+\sqrt{1+\frac{2r_0}{3M_{(5)}^3}\rho _b}\right].$$
(24)
The effect of including the intrinsic curvature term is only felt when the Hubble parameter is larger than some critical threshold determined by the parameter $`r_0^1`$. Then the cosmological behavior is 4-dimensional. Eventually, the Hubble parameter decreases below this threshold, and the system behaves as though the intrinsic curvature term were not present. The system inevitably asymptotes to the static solutions under consideration. This can be seen solving explicitly the brane Friedmann equation, Eq. (24).

Assuming again that the brane content is given by Eqs. (18) and taking $`k=0`$, one can solve explicitly for the scale factor. Let us define some useful parameters:

$`B`$ $`=`$ $`1{\displaystyle \frac{4}{3}}r_0^2\mu ^2`$ (25)
$`x`$ $`=`$ $`1+{\displaystyle \frac{4}{3}}r_0^2\mu ^2\left({\displaystyle \frac{\rho _b}{|\lambda |}}\right),`$ (26)
where we have used the mass scale $`\mu `$ defined in the last section. Three distinct cases emerge. If $`B<0`$,

$$\frac{q}{2r_0}(1B)(tt_0)=\mathrm{ln}\left|\frac{xB}{x1}\frac{1+\sqrt{x}}{1\sqrt{x}}\right|2\sqrt{B}\mathrm{tan}^1\left(\frac{\sqrt{x}}{\sqrt{B}}\right).$$
(27)
If $`B>0`$ and $`B1`$,

$$\frac{q}{2r_0}(1B)(tt_0)=\mathrm{ln}\left|\frac{xB}{x1}\frac{1+\sqrt{x}}{1\sqrt{x}}\left(\frac{\sqrt{B}\sqrt{x}}{\sqrt{B}+\sqrt{x}}\right)^\sqrt{B}\right|.$$
(28)
Finally, if $`B=1`$,

$$\frac{q}{2r_0}(tt_0)=\frac{\sqrt{x}}{1x}+\frac{1}{2}\mathrm{ln}\left|\frac{1+\sqrt{x}}{1\sqrt{x}}\right|+\frac{1}{x1}.$$
(29)
When the brane cosmological constant is negative, one can check that one recovers the asymptotic behavior mentioned previously. In this particular case one has $`B<1`$, and $`B0`$, the late time behavior is then given by considering the limit $`x1`$ in the above equations. On the other hand, the early time limiting behavior is given by $`x\mathrm{}`$ and can easily be seen to match with standard cosmology.

The three cases correspond to when the energy scale $`\mu `$ is larger than the crossover energy scale ($`r_0^1`$), when energy scale of the cosmological constant is smaller than the crossover scale, and when the cosmological constant is zero, respectively. If the crossover energy scale $`r_0^1`$ is much smaller than $`\mu `$, we see that a system with very large initial positive energy density evolves with a conventional power-law behavior. As the energy density decreases past the scale of the cosmological constant, the system evolves as if to bounce and recollapse in finite time (the acceleration approaches a constant). This behavior would occur if the Hubble parameter were dominated by the second term under the square root in Eq. (24). However, before that bounce can happen, the energy density reaches the scale $`r_0^1M_{(5)}^3`$, the Hubble parameter depends linearly on the energy density, and the system asymptotically approaches the static solution, and the bounce is averted.

### C Scalar field evolution

$$V(\varphi )=\frac{1}{2}m^2\left(\varphi ^2\varphi _0^2\right).$$
(30)
This asymptotic dynamics is generic for systems with a negative cosmological constant and matter on a braneworld. As another example of a system that has these elements, consider a scalar field, $`\varphi `$, living exclusively on the brane with the potential One can envision that this scalar field may be an inflaton, quintessence or other such cosmological field. The constant offset acts as a negative cosmological constant, with $`\lambda =m^2\varphi _0^2/2`$. Fluctuations in the field around the true vacuum state act as the matter. A more general potential may always be considered, but this example captures the qualitative features of interest. Consider the brane homogeneously filled with a scalar field possessing such a potential. The equations describing the scalar field evolutions is

$$\ddot{\varphi }+3H\dot{\varphi }+m^2\varphi =0,$$
(31)
and the Hubble parameter is determined by Eq. (24) when

$$\rho _b=\frac{1}{2}\dot{\varphi }^2+V(\varphi ).$$
(32)
Once again, the precise evolution of the system is dependent on two parameters: the crossover scale, $`r_0`$, and the quantity $`r_0^2\mu ^2`$, and the qualitative features are the same. Figures 1 and 2 depict the evolution of the scalar field, the scale factor, and the Hubble parameter of a typical case. With a positive energy density, the field oscillates with a damping proportional to the Hubble parameter. Eventually, the energy density is drawn toward zero, and the scalar field exponentially asymptotes to a dissipationless oscillation.

One expects that the oscillations in the scalar field amplify spatially inhomogeneous perturbations in that field. However, if the self-coupling of that field is small, such an amplification does not grow without bound. Also, unlike the matter considered earlier, the asymptotic metric of the spacetime is not static. Although the energy density vanishes, the pressure is a rapidly oscillating function.

## III Global structure of spacetime

As shown explicitly in  (see also) the bulk spacetime in Eq. (9) is two identical pieces of 5-dimensional Minkowski spacetime glued across the braneworld worldvolume. The same holds true for the asymptotic solution, Eq. (6), where the bulk is now a piece of Rindler spacetime. In this latter case the $`y=\mathrm{constant}`$ surfaces in the $`(\tau ,y)`$ plane are hyperbolas $`_y`$ in the two dimensional Minkowski space time, so that the brane space-time (6) is the hyperbolic cylinder $`𝒞`$ defined by $`_0\times E_3`$, where $`E_3`$ is the 3 dimensional Euclidian plane. The whole spacetime reflected in Eq. (6) is then simply given by gluing two copies of one side of $`𝒞`$ along $`𝒞`$.

This picture of the global structure of spacetime may be verified explicitly using a new set of coordinates $`Y^A`$ ($`A=0,1,2,3,5`$) of 5-dimensional Minkowski spacetime. The line element may be written

$$ds^2=(dY^0)^2+(dY^1)^2+(dY^2)^2+(dY^3)^2+(dY^5)^2,$$
(33)
where the coordinate change need to arrive at this line element is determined by equations given in .

We are mainly interested in the coordinate change for late times when the braneworld system approaches a static state. In this case, whether one considers an intrinsic curvature term Eq. (2) or not, we may use the solution in Eq. (20) for the scale factor. Equation (20) may be rewritten as

$$a=a_{\mathrm{}}(1ϵ)^{1/q},$$
(34)
with $`a_{\mathrm{}}`$ a constant given by

$$a_{\mathrm{}}=\left(\frac{\rho _0}{\lambda }\right)^{1/q}$$
(35)
is the $`t\mathrm{}`$ limit of the scale factor if $`a(t=t_0)=1`$. The quantity $`ϵ`$ defined by

$$ϵ\left(1+\frac{\lambda }{\rho _0}\right)e^{(tt_0)/T},$$
(36)
which is small for asymptotically late times, and where the time scale $`T`$ is given by

$$T6M_{(5)}^3/\lambda q.$$
(37)
Using the formulas given in , one finds at leading order in $`ϵ`$

$`Y^0+Y^5{\displaystyle \frac{a_{\mathrm{}}r^2}{2}}`$ $`=`$ $`ϵ^1{\displaystyle \frac{qT^2}{2a_{\mathrm{}}}}\left(1+{\displaystyle \frac{|y|}{T}}\right),`$ (38)
$`Y^0Y^52a_{\mathrm{}}`$ $`=`$ $`ϵ{\displaystyle \frac{2a_{\mathrm{}}}{q}}\left(1+{\displaystyle \frac{|y|}{T}}\right),`$ (39)
$`Y^i`$ $`=`$ $`a_{\mathrm{}}x^i,`$ (40)
where $`r=x^ix_i`$. One then easily sees that up to terms of order $`ϵ`$ (and a change in the origin of the Y-coordinate system), one has for a given $`r=\mathrm{constant}`$ slice

$$(Y^0)^2(Y^5)^2=T^2\left(1+\frac{|y|}{T}\right)^2.$$
(41)
For asymptotically late times, the world line of the braneworld follows this hyperbola. At specific time slices, the braneworld coordinates, $`Y^0`$ and $`Y^5`$ evolve toward infinity asymptoting the light cone, and the values of those coordinates grow exponentially with cosmological (braneworld) time $`t`$. Figure 3 depicts the global 5-dimensional Minkowski space for a given spatial slice.

## IV Four-dimensional gravity

### A Propagator on the brane

We can now study the propagation of gravitons in this asymptotic static space. The line element reads as follows:

$$ds^2=\left(1+r_0\mu ^2|y|\right)^2d\tau ^2+dx_idx^i+dy^2.$$
(42)
For simplicity we treat the graviton as a scalar particle. As was emphasized in , the scalar part of the graviton propagator reduces to that of a minimally coupled spin-0 particle in various backgrounds which we study below.

Let us start with the action:

$$d^4x𝑑y\left(\frac{M_{(5)}^3}{2}\sqrt{\stackrel{~}{g}}\stackrel{~}{g}^{AB}_A\varphi _B\varphi +\frac{M_{(4)}^2}{2}\delta (y)\sqrt{g}g^{ab}_a\varphi _b\varphi \right).$$
(43)
The first term here is a counterpart of the bulk Einstein action and the second one accounts for the induced curvature term. The equation for the Green function takes the form:

$$\left(M_{(5)}^3_A\sqrt{\stackrel{~}{g}}^A+\delta (y)M_{(4)}^2_b\sqrt{g}^b\right)G(x_a,y)=M_{(5)}^3\sqrt{g}\delta ^{(4)}(x)\delta (y).$$
(44)
Using (42) we find:

$`\left\{{\displaystyle \frac{1}{(1+r_0\mu ^2|y|)}}_\tau ^2+(1+r_0\mu ^2|y|)_i^2+_y(1+r_0\mu ^2|y|)_y\right\}G(x_a,y)`$ $`+`$ (45)
$`2r_0\delta (y)\left\{{\displaystyle \frac{1}{(1+r_0\mu ^2|y|)}}_\tau ^2+(1+r_0\mu ^2|y|)_i^2\right\}G(x_a,y)`$ $`=`$ $`\delta ^{(4)}(x)\delta (y).`$ (46)
To find the Green’s function it is convenient to perform Fourier transform to momentum space with respect to the four worldvolume coordinates $`x_a`$,

$$G(\tau ,x_i,y)\frac{d\omega d𝐪_i}{(2\pi )^4}e^{i\omega \tau }e^{i𝐪𝐱}\stackrel{~}{G}(\omega ,𝐪_i,y).$$
(47)
As a result, the equation takes the form:

$`\left\{{\displaystyle \frac{\omega ^2}{(1+r_0\mu ^2|y|)}}(1+r_0\mu ^2|y|)𝐪_i^2+_y(1+r_0\mu ^2|y|)_y\right\}\stackrel{~}{G}(\omega ,q,y)`$ $`+`$ (48)
$`2r_0\delta (y)\left\{\omega ^2𝐪_i^2\right\}\stackrel{~}{G}(\omega ,q,y)`$ $`=`$ $`\delta (y).`$ (49)
Let us introduce the following notation:

$$\stackrel{~}{G}\frac{1}{\sqrt{1+r_0\mu ^2|y|}}\chi .$$
(50)
Making use of this relation, the resulting equation reads as follows:

$$\chi ^{\prime \prime }+\left[\frac{\omega ^2+(r_0\mu ^2)^2/4}{(1+r_0\mu ^2|y|)^2}𝐪_i^2\right]\chi +2r_0\delta (y)\left\{\omega ^2𝐪_i^2\right\}\chi =\delta (y)\left[1+r_0\mu ^2\chi (0)\right].$$
(51)
To find the solution of this equation with appropriate boundary conditions let us introduce the following substitution:

$$\chi (\omega ,𝐪_i,y)=B(\omega ,𝐪_i)\mathrm{\Sigma }(y),$$
(52)
where $`\mathrm{\Sigma }`$ satisfies the following equation:

$$\mathrm{\Sigma }^{\prime \prime }+\left[\frac{\omega ^2+(r_0\mu ^2)^2/4}{(1+r_0\mu ^2|y|)^2}𝐪_i^2\right]\mathrm{\Sigma }=\delta (y).$$
(53)
From this definition, as well as from Eq. (52), one sees that the expression for the Green’s function on the brane ($`y=0`$) takes the form:

$$\stackrel{~}{G}|_{y=0}=\left[2r_0q^2r_0\mu ^2+\mathrm{\Sigma }(0)^1\right]^1,$$
(54)
where we have defined $`q^2=\omega ^2+𝐪_i^2`$, the invariant 4-momentum squared.

The solution to Eq. (53) with asymptotically vanishing behavior as $`|y|\mathrm{}`$ is

$$\mathrm{\Sigma }(y)\sqrt{u}K_\nu (u/2),$$
(55)
with

$$u=\frac{1+r_0\mu ^2|y|}{r_0\mu ^2}2|𝐪_i|,$$
(56)
and $`K_\nu `$ is the MacDonald function of order $`\nu =\pm i\omega /r_0\mu ^2`$. The normalization of $`\mathrm{\Sigma }(y)`$ may be established by the delta function in the right-hand side of Eq. (53), indicating that

$$\mathrm{\Sigma }^{}(y=0^+)=\frac{1}{2}.$$
(57)
One obtains

$$\mathrm{\Sigma }(0)^1=2|𝐪_i|\left[\frac{12\nu }{u_0}\frac{K_{\nu 1}(u_0/2)}{K_\nu (u_0/2)}\right],$$
(58)
where we define

$$u_0u(y=0)=\frac{2|𝐪_i|}{r_0\mu ^2}.$$
(59)
Using Eq. (54), one finds the exact expression of the propagator on the brane:

$$\stackrel{~}{G}|_{y=0}=\frac{1}{2r_0}\left[q^2\frac{i\omega }{r_0}\frac{|𝐪_i|}{r_0}\frac{K_{1\pm i\omega /r_0\mu ^2}\left(|𝐪_i|/r_0\mu ^2\right)}{K_{\pm i\omega /r_0\mu ^2}\left(|𝐪_i|/r_0\mu ^2\right)}\right]^1.$$
(60)
In the above expression the first term on the right hand side represents the usual 4-dimensional Lorentz invariant propagator, whereas the two other terms are responsible for deviation from 4-dimensional behavior as well as for Lorentz violating effects.Because the MacDonald function is even with respect to its order, one can indeed verify that one arrives at the same result regardless of the sign taken in front of $`\omega `$ in Eq. (60).

### B Asymptotic developments

We wish to elaborate on the propagator Eq. (60) in some special limits. We see that the propagator is 4-dimensional Lorentz symmetry violating in general. Let us first examine the limit in which the spatial momentum $`𝐪_i0`$, with $`\omega 0`$. Taking the propagator with outgoing radiative boundary conditions, at leading order

$$\stackrel{~}{G}|_{y=0}^1=2r_0\omega ^22i\omega .$$
(61)
We see then a restoration of the propagator found in  in the same limit, indicating there is no massive graviton state stationary with respect to the cosmological frame.

Next we wish to take the limit where the propagator yields the static potential, $`\omega 0`$.

$$\stackrel{~}{G}|_{y=0}^1=2r_0𝐪_i^22|𝐪_i\frac{K_1(|𝐪_i|/r_0\mu ^2)}{K_0(|𝐪_i|/r_0\mu ^2)}.$$
(62)
Note that this form is not covariant with respect to the form found in Eq. (61). For large distances ($`|𝐪_i|0`$), the propagator reduces to

$$\stackrel{~}{G}|_{y=0}=\frac{1}{2r_0\mu ^2}\mathrm{ln}\frac{|𝐪_i|}{2r_0\mu ^2}.$$
(63)
This form corresponds to a $`r^3`$ static potential. For short distances, ($`|𝐪_i|\mathrm{}`$), the propagator reduces to

$$\stackrel{~}{G}|_{y=0}=\frac{1}{2r_0𝐪_i^2},$$
(64)
corresponding to a $`r^1`$ static potential. Crossover behavior varies depending on the value of the dimensionless parameter, $`r_0\mu `$.

In order to see the restoration of 4-dimensional Lorentz invariance, we must go to large values of $`\omega `$ and $`𝐪_i`$. Taking $`r_0\mu ^20`$ while holding $`\omega `$ and $`𝐪_i`$ fixed, the propagator becomes

$$\stackrel{~}{G}|_{y=0}^1=2r_0q^22qr_0\mu ^2(1+\omega ^2/q^2).$$
(65)
Note that when $`r_0\mu ^20`$, explicit 4-dimensional Lorentz invariance is recovered and the propagator reduces to the form found in . Equation (65) expresses the dominant off-shell behavior for the scalar graviton propagator.

Alternatively, when $`q^20`$, but in a manner such that $`\omega ^2r_0^2\mu ^4`$ with $`q/\omega <(r_0\mu ^2/\omega )^{1/3}`$, then

$$\stackrel{~}{G}|_{y=0}^1=2r_0q^22\omega \frac{\mathrm{\Gamma }(2/3)}{\mathrm{\Gamma }(1/3)}\left(\frac{6r_0\mu ^2}{\omega }\right)^{1/3}.$$
(66)
The pole in the propagator in this regime indicates the following dispersion relationship for an on-shell particle:

$$\omega ^2=𝐪_i^2+\frac{\mathrm{\Gamma }(2/3)}{\mathrm{\Gamma }(1/3)}\left(6r_0\mu ^2\right)^{1/3}(𝐪_i^2)^{1/3}.$$
(67)
As $`|𝐪_i|\mathrm{}`$, the first term in the expression is dominant and the expected relativistic dispersion relationship is recovered.

## V Concluding Remarks

Our system has the property of having energy density evolve to zero from a generic set of initial conditions with positive energy density. The pressure of this asymptotic state is non-zero, so the full 5-dimensional spacetime explicitly breaks 4-dimensional Lorentz invariance, even while the induced metric on the braneworld is Minkowskian. The graviton propagator on such a braneworld also reflects a lack of 4-dimensional Lorentz symmetry, since that graviton can probe the full spacetime.

A static Minkowski braneworld coupled to a braneworld energy-momentum which explicitly violates 4-dimensional Lorentz symmetry is possible because the larger metric of the 5-dimensional space does not respect 4-dimensional Lorentz symmetry. Any slice of the 5-dimensional spacetime parallel to the braneworld is indeed Minkowski. However, the speed of light varies from slice to slice. If standard model particles trapped on the braneworld do not couple to the dark matter generating the 5-dimensional spacetime, then there should be no observable Lorentz symmetry violation in the standard model, except through graviton exchange.

One can ask how well this model complies with known phenomenology. We have two parameters that may be adjusted: the crossover scale, $`r_0=M_{(4)}^2/2M_{(5)}^3`$, which reflects the relative sizes of the 4-dimensional Planck mass and the 5-dimensional Planck mass, as well as $`\mu ^2=\lambda /M_{(4)}^2`$, which reflects the relative size of the brane cosmological constant to the 4-dimensional Planck scale.

We may place constraints on the system parameters at various stages. Let us first address solely the constraints from gravitational force law observations. The exchange of the tower of Kaluza-Klein (KK) modes may be understood from the 4-dimensional point of view as an exchange of a metastable graviton . The decay width is controlled by the larger of the two energy scales $`r_0^1`$ or $`\mu `$. Thus, both energy scales must be tuned to some sufficiently small value to be consistent with a stable, massless graviton. We know that the graviton is stable and massless up until $`100`$ parsecs (for a discussion, see ). We may choose the galactic size as well, but for here we settle for a strict confirmation of the $`1/r^2`$ force law. Then both $`r_0`$, and $`\mu ^1`$ must be as large as $`100`$ parsecs. This places constraints on $`M_{(5)}<100\mathrm{GeV}`$ and $`|\lambda |<(10\mathrm{eV})^4`$.

We may also constrain our system parameters based on the standard cosmological model. We found that early cosmology only follows the conventional hot big bang scenario when the energy density on the brane is much larger than $`|\lambda |`$ as well as a similar constraint between the Hubble parameter and the crossover scale $`r_0^1`$. Cosmology constraints imply $`M_{(5)}<100\mathrm{MeV}`$ and $`|\lambda |<(10^3\mathrm{eV})^4`$.

Having such a small fundamental Planck scale, $`M_{(5)}`$, does not contradict any astrophysical or particle physics bounds as long as there is a 4-dimensional induced Planck scale on the brane (see ), assuming quantum gravity becomes soft above this energy scale. However, in the present context this issue deserves clarification, due to the presence of a Lorentz symmetry violating metric which is not considered in . Since the Lorentz violation is due to a tiny cosmological constant, the prior analysis should be valid, and any new effect should be a small perturbation. On the other hand, the warp factor (which parametrizes Lorentz symmetry violation) diverges as $`M_{(5)}0`$. While this statement is true, Lorentz symmetry violation in the propagator on the brane is controlled by the second and thirds terms on the right-hand side of Eq. (60). For small $`M_{(5)}`$ (when $`M_{(4)}`$ and $`\lambda `$ are fixed) those terms are suppressed relative to the leading 4-dimensional contribution to the propagator. In other words, although the Lorentz violating warp factor becomes steep when $`M_{(5)}`$ is small, it cannot strongly affect brane observers since probing the bulk becomes more difficult. This can be also understood in the KK-mode expansion. In this language, the intrinsic curvature term repels KK-modes heavier than $`1/r_0`$, suppressing their wavefunction on the brane and shielding the brane from bulk gravity effects, in particular Lorentz symmetry violation.

All these observations imply that one cannot start with a large cosmological constant on the brane, and therefore this model cannot yield a resolution of the cosmological constant problem. The mechanism that so generously compels the total energy density to zero is equally unforgiving. Any combination of matter whose total energy density is on the order of magnitude of the cosmological constant evolves exponentially fast to a flat and perfectly static universe, in a time comparable to $`M_{(5)}^3/|\lambda |`$.

We have proposed a model which has some intriguing properties. It dynamically evolves to a zero energy-density state and provides a novel mechanism for producing a static, Lorentz symmetry violating universe. Perhaps these mechanisms may be applied to more phenomenologically promising models in the future.

###### Acknowledgements.

We would like to thank C. Grojean and D. Hogg for discussions. C. D. and G. G. thank the LBL Theory Group for its hospitality where part of this work was done. This work is sponsored in part by NSF Award PHY-9803174 and PHY-0070787, David and Lucille Packard Foundation Fellowship 99-1462, the Alfred P. Sloan Foundation Fellowship, and DOE Grant DE-FG02-94ER408.


text:
# Comment on “Long-range electrostatic interactions between like-charged colloids: Steric and confinement effects”

## Abstract

In a recent study \[Phys. Rev. E 60, 6530 (1999)\], Trizac and Raimbault showed that the effective pair interaction between like charged colloids immersed in a cylindrically confined electrolyte remains repulsive even when the size of the micro-ions or the finite longitudinal extension of the confining cylinder are taken into account. Contrary to their claim, we argue that the case of finite longitudinal confinement doesn’t always generate repulsive interactions and to illustrate this point we also provide a simple example.

There have been a number of recent theoretical efforts towards understanding the mechanism behind the long range attraction that is experimentally observed  when like charged colloids are immersed in a confined electrolyte. The fact that this interaction returns to a purely repulsive one as soon as the confining surface is removed , suggested that the attractive part could be obtained within the framework of the Poisson-Boltzmann theory by properly taking into account the influence of the confining surface. However, Neu  and Sader and Chan , using a general model of cylindrical confinement, have rigorously proven that this is not the case.

In their paper, Trizac and Raimbault  extended the aforementioned studies  to include the influence brought upon the sign of the interaction by the size of the micro-ions and the finite longitudinal extension of the confining cylinder (case hereafter referred to as complete confinement). Although we agree with their conclusion that the interaction remains repulsive when the specific size of the micro-ions is taken into account through a modified Poisson-Boltzmann equation, we don’t agree with a similar conclusion for the case of complete confinement. In this comment, we show that this controversy is brought about by an inaccuracy in their calculation and that in general, colloids immersed in a completely confined electrolyte don’t always repel. We also provide a simple example to illustrate this point.

To be more specific, the issue has to do with the use of Green’s theorem in two dimensions ($`2D`$). If $`W`$ is a closed domain in $`𝐑^\mathrm{𝟐}`$ (for example in the $`xy`$ plane) and $`\varphi (x,y)`$ and $`\psi (x,y)`$ are scalar fields defined on $`𝐑^\mathrm{𝟐}`$, then Green’s theorem in $`2D`$ is:

$$_W𝑑x𝑑y\left[\left(_{xy}\varphi \right)\left(_{xy}\psi \right)+\varphi _{xy}^2\psi \right]=_W𝑑\mathrm{}\varphi \left(\widehat{𝐧}_{xy}\psi \right),$$
(1)
where $`\widehat{𝐧}`$ (which is contained in the $`xy`$ plane) represents the outer unit normal to the boundary $`W`$ of $`W`$ and $`_{xy}`$ is the $`2D`$ gradient operator.

With this in mind and making use of the same notation and assumptions as in , we obtain a different expression for Eqs. (14) and (15) in that paper, where Green’s theorem was inaccurately applied. For clarity, let us first separate the $`z`$ component of the gradient in the equation that generates the results $`(14)`$ and $`(15)`$ in :

$`{\displaystyle _{Oxy}}𝑑x𝑑y𝐄_{z=L}\left(𝐃_{z=0}𝐃_{z=L}\right)`$ (2)
$`={\displaystyle _{Oxy}}𝑑x𝑑y\left(\psi \right)_{z=L}\left[ϵ\left(\psi \right)_{z=0}ϵ\left(\psi \right)_{z=L}\right]`$ (3)
$`={\displaystyle _{Oxy}}𝑑x𝑑y\left\{ϵ\left(_{xy}\psi \right)_{z=L}\left[\left(_{xy}\psi \right)_{z=0}+\left(_{xy}\psi \right)_{z=L}\right]ϵ\left({\displaystyle \frac{\psi }{z}}\right)_{z=L}^2\right\},`$ (4)
where $``$ is the $`3D`$ gradient operator and $`\psi `$ possesses mirror symmetry with respect to the plane $`z=0`$, as in . The reader should remember that the notation $`_{Oxy}`$ actually stands for the two integrals over the domains of the $`xy`$ plane that are interior and exterior, respectively, to the confining cylinder. The dimensionality of these integrals requires the use of the $`2D`$ version of Green’s theorem, unlike it was inaccurately done in , where the $`3D`$ version was employed instead. This issue is clarified if we consider the more explicit representation shown in Eq. (4), where we can correctly apply the $`2D`$ version of Green’s theorem \[as stated in Eq. (1)\] to its first term to get:

$`{\displaystyle _\mathrm{\Sigma }}𝑑\mathrm{}\left(D_nD_n^{}\right)_{z=L}\left(\psi _{z=o}+\psi _{z=L}\right){\displaystyle _{Oxy}}𝑑x𝑑y\left[ϵ\left(_{xy}^2\psi \right)_{z=L}\left(\psi _{z=o}\psi _{z=L}\right)+ϵ\left({\displaystyle \frac{\psi }{z}}\right)_{z=L}^2\right]`$ (5)
$`=\sigma {\displaystyle _\mathrm{\Sigma }}𝑑SE_z{\displaystyle _{Oxy}}𝑑x𝑑y\left[\left(\psi _{z=o}\psi _{z=L}\right){\displaystyle \frac{P}{\psi }}\left(\psi _{z=L}\right)\left(\psi _{z=o}\psi _{z=L}\right)ϵ\left({\displaystyle \frac{^2\psi }{z^2}}\right)_{z=L}+ϵ\left(E_z\right)_{z=L}^2\right],`$ (6)
where in going from Eq. (5) to Eq. (6) we used the dielectric boundary condition $`D_n^{}D_n=\sigma `$ (and the fact that $`_\mathrm{\Sigma }𝑑\mathrm{}\left(\psi _{z=0}\psi _{z=L}\right)=_\mathrm{\Sigma }𝑑SE_z`$) and also formed the $`3D`$ Laplacian and employed Eq. (10) in .

In fact, Eqs. (5) and (6) are the correct forms for Eqs. (14) and (15), respectively, in . After adding all the terms, the correct expression for the $`z`$ component of the force (Eq. (16) in ) should read:

$`F_z=`$ $`{\displaystyle _{Oxy}}𝑑x𝑑y\left[P\left(\psi _{z=0}\right)P\left(\psi _{z=L}\right)\left(\psi _{z=0}\psi _{z=L}\right){\displaystyle \frac{P}{\psi }}\left(\psi _{z=L}\right)\right]`$ (8)
$`+{\displaystyle _{Oxy}}𝑑x𝑑y{\displaystyle \frac{ϵ}{2}}\left(𝐄_{z=0}𝐄_{z=L}\right)^2+{\displaystyle _{Oxy}}𝑑x𝑑y\left(\psi _{z=0}\psi _{z=L}\right)ϵ\left({\displaystyle \frac{^2\psi }{z^2}}\right)_{z=L}.`$
Although one can use the same argument as in  to prove that the first two integrals in Eq. (8) are always positive, the third one (which is different from that obtained by Trizac and Raimbault) can in fact be negative and sometimes overcome the first two. Before we show how this can happen, let us mention another equivalent expression for $`F_z`$, which can be obtained by recasting in a different manner the second integrand of Eq. (13) in , namely $`\left(𝐃𝐄\right)_{z=0}\left(𝐃𝐄\right)_{z=L}=ϵ\left(𝐄_{z=0}𝐄_{z=L}\right)^22𝐄_{z=0}\left(𝐃_{z=L}𝐃_{z=0}\right)`$ and then following the same steps used in deriving Eq. (8):

$`F_z=`$ $`{\displaystyle _{Oxy}}𝑑x𝑑y\left[P\left(\psi _{z=L}\right)P\left(\psi _{z=0}\right)\left(\psi _{z=L}\psi _{z=0}\right){\displaystyle \frac{P}{\psi }}\left(\psi _{z=0}\right)\right]{\displaystyle _{Oxy}}𝑑x𝑑y{\displaystyle \frac{ϵ}{2}}\left(𝐄_{z=0}𝐄_{z=L}\right)^2`$ (10)
$`+{\displaystyle _{Oxy}}𝑑x𝑑yϵ\left(E_z\right)_{z=L}^2{\displaystyle _{Oxy}}𝑑x𝑑y\left(\psi _{z=L}\psi _{z=0}\right)ϵ\left({\displaystyle \frac{^2\psi }{z^2}}\right)_{z=0}.`$
The reader can easily check that Eqs. (8) and (10) are equivalent, by subtracting them from one another and then using Eq. (10) in , Green’s theorem in $`2D`$ and the dielectric boundary condition to prove that the result is zero.

We can now propose a simple example in order to prove that the complete confinement doesn’t always generate repulsive interactions. As before, we assume that the particles are immersed in an electrolyte of permittivity $`ϵ`$, confined by an infinite and uniformly charged cylindrical surface (of arbitrary cross section), while the medium outside the cylinder has a permittivity $`ϵ^{}`$ and also contains an electrolyte (most general case). The only requirement is that $`\psi `$ possesses mirror symmetry with respect to the plane $`z=0`$. In our example the particles $`𝐒_1`$ and $`𝐒_2`$ are taken to be flat surfaces with the same surface charge density $`\sigma `$ as the confining cylinder while the micro-ions are able to equilibrate at all times (by being in contact with a reservoir and possibly permeating through the colloidal particles). The walls A and B that realize the complete confinement are surrounded on both sides by electrolyte and also have the same uniform surface charge density $`\sigma `$ as the confining cylinder (this way the complete confinement is realized by a surface of constant $`\sigma `$). To compute the z component of the force on the colloid $`𝐒_2`$, we follow the same procedure outlined in  only that instead of choosing the surface $`\mathrm{\Sigma }_L`$ to integrate over in Eq. (11) that defines $`F_z`$ in , we choose the surface $`\mathrm{\Sigma }_d`$, midway between the colloid $`𝐒_2`$ and the wall B (see Fig. 1). Keeping in mind the correction made in this comment, the result ends up to be identical to Eq. (10), with the only exception that all the indices $`z=L`$ are replaced by $`z=d`$:

$`F_z=`$ $`{\displaystyle _{Oxy}}𝑑x𝑑y\left[P\left(\psi _{z=d}\right)P\left(\psi _{z=0}\right)\left(\psi _{z=d}\psi _{z=0}\right){\displaystyle \frac{P}{\psi }}\left(\psi _{z=0}\right)\right]{\displaystyle _{Oxy}}𝑑x𝑑y{\displaystyle \frac{ϵ}{2}}\left(𝐄_{z=0}𝐄_{z=d}\right)^2`$ (12)
$`+{\displaystyle _{Oxy}}𝑑x𝑑yϵ\left(E_z\right)_{z=d}^2{\displaystyle _{Oxy}}𝑑x𝑑y\left(\psi _{z=d}\psi _{z=0}\right)ϵ\left({\displaystyle \frac{^2\psi }{z^2}}\right)_{z=0}.`$
It is easy to see that the last two integrals in Eq. (12) approach zero as the separation $`2L`$ between the walls is increased to infinity and the distance $`2(Ld)`$ between each particle and the adjacent wall is kept constant. Indeed, in this situation, $`\psi `$ becomes symmetric with respect to the plane $`z=d`$ (due to our choice for the shape of the colloidal particle) and therefore $`\left(E_z\right)_{z=d}0`$, while at $`z=0`$ (infinitely far away from the particles and the walls), the problem of solving for $`\psi `$ becomes essentially a $`2D`$ one (independent of $`z`$) and therefore $`^n\psi /z^n0`$ for any $`n1`$. At the same time, as argued in , the first two integrals in Eq. (12) are always positive (the first integrand is positive due to the global convexity of $`P(\psi )`$) and moreover, reach some finite asymptotic values in this limit case. Since both integrals contribute with a minus sign in Eq (12), it becomes clear that $`F_z<0`$ when $`L\mathrm{}`$ and $`Ld`$ is constant.

Therefore, within the approximations of our example, we conclude that for any fixed distance $`2(Ld)`$ between the colloidal particle and the adjacent wall, there is a finite, maximum separation between the walls $`2L_0`$, at which the particles are in equilibrium ($`F_z=0`$) and beyond which $`F_z`$ becomes negative and remains so as $`L`$ increases. Although this is in contradiction with the result obtained earlier by Trizac and Raimbault  that $`F_z`$ should always be positive under complete confinement, we have shown that the controversy stems from an inaccuracy in their calculation. It is therefore possible that completely confined colloidal particles don’t always repel, mainly due to the stronger interaction with the walls (or conceivably other particles ) that are responsible for the finite longitudinal confinement and break the translational symmetry of the confined space. Although this type of calculation is not suitable to yield an estimate on $`L_0`$, it would nevertheless be of high interest if this could be done through computer simulations.

###### Acknowledgements.

It is a pleasure to acknowledge stimulating discussions with M. N. Tamashiro, C. Jeppesen and P. Pincus. This work was supported by the DMR Program of the National Science Foundation under Awards DMR99-72246 and DMR96-32716.