text:
# Self-similarity limits of genomic signatures

## I Introduction

Along with an increasing amount of DNA sequences extracted from experiments, it is important to develop methods for extracting meaningful information from the one-dimensional symbolic sequences composed of the four letters ‘A’, ‘C’, ‘G’ and ‘T’ (or ‘U’). To detect similarity in DNA sequences, scatter plots are introduced to make classification of cytochromes and illustrate a dendrogram. From a comparison of a pair of duplicated genes by a distance matrix, evolutionary relationship of three primary kingdoms of life is inferred. Due to investigating relative abundances of short oligonucleotides in subsequences, genomic signature phenomenon and derivation of partial-ordering relationships among bacterial genomes are proposed. The genomic signature describes that the difference of dinucleotide relative abundance values within a single genome is larger than that between distinct genomes. Chaos game representation (CGR), which generates a two-dimensional square from a one-dimensional sequence, provides a technique to visualize the composition of DNA sequences. By composing the CGR and short-sequence representation methods, the evolution of species-type specificity in mitochondral genomes is analyzed. In terms of the CGR method, it is shown that the main characteristics of the whole genome can be exhibited by its subsequences. The genomic signature is extended to describe characteristics of CGR images. By making a Euclidean metric between two CGR images, classification of species in three primary kingdoms is discussed.

Recently, metric representation (MR), which is borrowed from the symbolic dynamics, makes an ordering of subsequences in a plane. The MR method is an extension of CGR. Suppression of certain nucleotide strings in the DNA sequences leads to a self-similarity of pattern seen in the MR of DNA sequences. In this paper, first, we show that the MR is one-to-one. Due to the MR method, we determine suppression of nucleotide strings in DNA sequences. Then, eliminating effects of finite sequences on suppression of nucleotide strings, we give an optimal string length to display genomic signature. Moreover, we plot information function versus string lengths to determine self-similarity limits in MR images. Using the method, we present self-similarity limits of bacteria complete genomic signatures.

## II Suppression of nucleotide strings

For a given DNA sequence, we have a one-dimensional symbolic sequence $`s_1s_2\mathrm{}s_i\mathrm{}s_N`$ ($`s_i\{A,C,G,T\}`$). In a two-dimensional MR, we take the correspondence of symbol $`s_i`$ to number $`\mu _i`$ or $`\nu _i\{0,1\}`$ and calculate the values ($`\alpha `$, $`\beta `$) of all subsequences $`\mathrm{\Sigma }_m=s_1s_2\mathrm{}s_m`$ ($`1mN`$). The number $`\alpha `$ represented in base 3, between 0 and 1, is defined as

$$\alpha =2\underset{j=1}{\overset{m}{}}\mu _{mj+1}3^j+3^m=2\underset{i=1}{\overset{m}{}}\mu _i3^{(mi+1)}+3^m,$$
(1)
where $`\mu _i`$ is 0 if $`s_i\{A,C\}`$ or 1 if $`s_i\{G,T\}`$. Similarly, the number $`\beta `$ is defined as

$$\beta =2\underset{j=1}{\overset{m}{}}\nu _{mj+1}3^j+3^m=2\underset{i=1}{\overset{m}{}}\nu _i3^{(mi+1)}+3^m,$$
(2)
where $`\nu _i`$ is 0 if $`s_i\{A,T\}`$ or 1 if $`s_i\{C,G\}`$. According to (1) and (2), the one-dimensional symbolic sequence $`s_1s_2\mathrm{}s_N`$ is partitioned into 4 kinds of subsequences, which correspond to points in 4 fundamental zones $`A,C,G`$ and $`T`$ of Fig. 1. Under left or right shift operators, each zone can be further shrunk to less zones with a factor of $`1/3^2`$. For an infinite sequence, this procedure can be defined as a fractal, which has a self-similarity. The subsequences with the same ending $`k`$-nucleotide string are labeled by $`\mathrm{\Sigma }_{(k)}`$. All subsequences $`\mathrm{\Sigma }_{(k)}`$ correspond to points in the zone encoded by the $`k`$-nucleotide string.

Lemma 1 $`(\alpha ,\beta )\{S(\mathrm{\Sigma }_m)\}=2(\mu _{m+1},\nu _{m+1})/3+(\alpha ,\beta )\{\mathrm{\Sigma }_m\}/3`$. $`S`$ is a left shift operator.

PROOF: Note that for the left shift operator, $`S(\mathrm{\Sigma }_m)=\mathrm{\Sigma }_ms_{m+1}`$. From the definition (1) and (2), we can immediately obtain the result.

Lemma 2 $`(\alpha ,\beta )\{\mathrm{\Sigma }_m\}=(\alpha ,\beta )\{G^{\mathrm{}}\mathrm{\Sigma }_m\}`$.

PROOF: When $`m=1`$, $`\mathrm{\Sigma }_1=s_1`$ and $`G^{\mathrm{}}\mathrm{\Sigma }_1=S(G^{\mathrm{}})`$. By the Lemma 1, we can obtain $`(\alpha ,\beta )\{G^{\mathrm{}}\mathrm{\Sigma }_1\}=2(\mu _1,\nu _1)/3+(\alpha ,\beta )\{G^{\mathrm{}}\}/3=2(\mu _1,\nu _1)/3+(1,1)/3=(\alpha ,\beta )\{\mathrm{\Sigma }_1\}`$. Suppose when $`m=i`$, we have $`(\alpha ,\beta )\{\mathrm{\Sigma }_i\}=(\alpha ,\beta )\{G^{\mathrm{}}\mathrm{\Sigma }_i\}`$. For $`m=i+1`$, we have $`\mathrm{\Sigma }_{i+1}=\mathrm{\Sigma }_is_{i+1}=S(\mathrm{\Sigma }_i)`$ and $`G^{\mathrm{}}\mathrm{\Sigma }_{i+1}=S(G^{\mathrm{}}\mathrm{\Sigma }_i)`$. By the Lemma 1, we obtain $`(\alpha ,\beta )\{\mathrm{\Sigma }_{i+1}\}=2(\mu _{m+1},\nu _{m+1})/3+(\alpha ,\beta )\{\mathrm{\Sigma }_i\}/3`$ and $`(\alpha ,\beta )\{G^{\mathrm{}}\mathrm{\Sigma }_{i+1}\}=2(\mu _{m+1},\nu _{m+1})/3+(\alpha ,\beta )\{G^{\mathrm{}}\mathrm{\Sigma }_i\}/3`$. So, using the supposition $`(\alpha ,\beta )\{\mathrm{\Sigma }_i\}=(\alpha ,\beta )\{G^{\mathrm{}}\mathrm{\Sigma }_i\}`$, we can lead to $`(\alpha ,\beta )\{\mathrm{\Sigma }_{i+1}\}=(\alpha ,\beta )\{G^{\mathrm{}}\mathrm{\Sigma }_{i+1}\}`$.

By the Lemma 2, each finite subsequence $`\mathrm{\Sigma }_m`$ has a correspondent infinite sequence $`G^{\mathrm{}}\mathrm{\Sigma }_m`$. Here, we define a set of the infinite sequences as $`\mathrm{\Sigma }`$.

Theorem 1 $`(\alpha ,\beta ):\mathrm{\Sigma }\mathrm{\Lambda }`$ is one-to-one. $`\mathrm{\Lambda }`$ is a set of points in the ($`\alpha ,\beta `$) plane.

This means that given $`\mathrm{\Sigma }^1`$, $`\mathrm{\Sigma }^2\mathrm{\Sigma }`$, if $`\mathrm{\Sigma }^1\mathrm{\Sigma }^2`$, then $`(\alpha ,\beta )\{\mathrm{\Sigma }^1\}(\alpha ,\beta )\{\mathrm{\Sigma }^2\}`$. We give a proof by contradiction. Suppose $`(\alpha ,\beta )\{\mathrm{\Sigma }^1\}=(\alpha ,\beta )\{\mathrm{\Sigma }^2\}`$ and is marked as $`P`$ in the the ($`\alpha ,\beta `$) plan. For the zone including the point $`P`$, we encode it as two subsequences $`\mathrm{\Sigma }_1^1`$ and $`\mathrm{\Sigma }_1^2`$ with the same mononucleotide. Then, enlarge the zone by a area factor of $`3^2`$, we can obtain two encoding subsequences $`\mathrm{\Sigma }_2^1`$ and $`\mathrm{\Sigma }_2^2`$ with the same dinucleotide. Each enlarging process provides a right shift to two subsequences. At the same time, the point $`P`$ is only included in one of four enlarged zones. So, two shifting subsequences are the same. Following the enlarging process in an infinite step, we can obtain $`\mathrm{\Sigma }^1=\mathrm{\Sigma }^2`$, contradicting our original assumption. This contraction is due to the fact that we have assumed $`(\alpha ,\beta )\{\mathrm{\Sigma }^1\}=(\alpha ,\beta )\{\mathrm{\Sigma }^2\}`$; thus, $`\mathrm{\Sigma }^1\mathrm{\Sigma }^2`$, then $`(\alpha ,\beta )\{\mathrm{\Sigma }^1\}(\alpha ,\beta )\{\mathrm{\Sigma }^2\}`$.

For the DNA sequence, some zones in CGR are replenished by points, so that a pattern appears. In CGR, there exists an correspondence of more subsequences with different ending $`k`$-nucleotide strings to the same points in bounds of zones. For examples, subsequences $`G^{\mathrm{}}A`$ in the zone $`A`$, $`T^{\mathrm{}}C`$ in the zone $`C`$, $`A^{\mathrm{}}G`$ in the zone $`G`$ and $`C^{\mathrm{}}T`$ in the zone $`T`$ have the same points in CGR (1/2,1/2). Under left shift operators, the corresponding relation between points and subsequences can preserve in zones with small enough lengths. For example, subsequences $`G^{\mathrm{}}AC`$ in the zone $`AC`$, $`T^{\mathrm{}}C^2`$ in the zone $`C^2`$, $`A^{\mathrm{}}GC`$ in the zone $`GC`$ and $`C^{\mathrm{}}TC`$ in the zone $`TC`$ have the same points in CGR (1/4,3/4). In MR of DNA sequences, each zone in CGR is shrunk and clearly divided by four bands. There exists a one-to-one correspondence between zones and ending $`k`$nucleotide strings of subsequences. Frequency of points in the zone can be determined by using MR method as follows. In order to compute frequencies in zones encoded by $`k`$-nucleotide strings, we need to determine partition lines of MR in Fig. 1. For mononucleotides, there exist $`2\times 2`$ zones in the MR. We have $`n_1(=3)`$ partition lines $`b_0^1=0`$, $`b_1^1=1/2`$ and $`b_2^1=1`$ along the $`\alpha `$ axis. For denucleotides, there exist $`4\times 4`$ zones in the MR. We have $`n_2(=5)`$ partition lines $`b_0^2=b_0^1=0`$, $`b_1^2=b_1^1/3=1/6`$, $`b_2^2=b_1^1=1/2`$, $`b_3^2=1b_1^2=5/6`$ and $`b_4^2=1b_0^2=1`$ along the $`\alpha `$ axis. In general, for $`k1`$-nucleotide strings, if knowing $`n_{k1}(=2^{k1}+1)`$ partition lines $`b_i^{k1}(i=0,1,\mathrm{},n_{k1}1)`$ along the $`\alpha `$ axis, we can obtain $`n_k(=2^k+1=2n_{k1}1)`$ partition lines $`b_i^k(i=0,1,\mathrm{},n_k1)`$ for $`k`$-nucleotide strings as follows. For the $`k`$-nucleotide strings, there exist $`2^k\times 2^k`$ zones in the MR. The left half $`(0in_{k1}1)`$ of partition lines along the $`\alpha `$ axis are described as follows

$$\begin{array}{cc}b_i^k=b_{i/2}^{k1},\hfill & \mathrm{if}i\%2=0;\hfill \\ b_i^k=b_i^{k1}/3,\hfill & \mathrm{if}i\%2=1.\hfill \end{array}$$
(3)
From (3), the right half $`(n_{k1}in_k1)`$ of partition lines along the $`\alpha `$ axis can be determined immediately

$$b_i^k=1b_{n_k1i}^k.$$
(4)
For example, for trinucleotides, 9 partition lines along the $`\alpha `$ axis are $`0,\frac{1}{18},\frac{1}{6},\frac{5}{18},\frac{1}{2},\frac{13}{18},\frac{5}{6},\frac{17}{18}`$ and 1. We can obtain 17 partition lines $`0,\frac{1}{54},\frac{1}{18}`$, $`\frac{5}{54},\frac{1}{6},\frac{13}{54},\frac{5}{18},\frac{17}{54}`$, $`\frac{1}{2},\frac{37}{54},\frac{13}{18},\frac{41}{54},\frac{5}{6}`$, $`\frac{49}{54},\frac{17}{18},\frac{53}{54}`$ and 1 along the $`\alpha `$ axis for tetranucleotides. Partition lines along the $`\beta `$ axis are the same to those along the $`\alpha `$ axis. Each zone in the MR can thus be surrounded by the combined partition lines along the $`\alpha `$ and $`\beta `$ axes.

Using the MR method, we determine suppression of $`k`$-nucleotide strings in HUMHBB (human $`\beta `$-region, chromosome 11) with 73308 bases and YEAST1 (yeast chromosome 1) with 230209 bases in Table I, respectively. In order to check efficiency of the method, we also determine the number of disappearing strings in all strings for a giving string length in HUMHBB and YEAST1, respectively. The results are identical with those in Table I, respectively. So, the MR method is effective to determine suppression of nucleotide strings in DNA sequences.

In CGR of DNA sequences, self-similarity patterns change more obscurely as lengths of sequences increase. A grey plot describes frequency values in small zones, which sizes ($`2^k\times 2^k`$) can be given by lengths of strings encoding the zones ($`k`$). Along with increase of string lengths, the self-similarity patterns in CGR are more clear. A high and low frequent zones are redivided to smaller and described by a grey scale. Some empty zones may appear in the patterns of CGR, i.e., some nucleotide strings are suppressed in the sequences. In the procedure of decreasing zone sizes, more and more empty zones emerge in the patterns of CGR. For example, evolution of a self-similarity pattern in CGR of the archaebacteria Archeoglobus fulqidus is shown in Fig. 1 of Ref. . If DNA sequences are infinite, the compositional structure can be displayed in small enough zones. Empty zones are a part of the global feature in CGR. But the DNA sequences are finite. A finite sequence, even a random sequence, may also lead to suppression of strings. Along with increase of string length, more and more strings are suppressed in the finite sequences.

In Table I, we compare the suppression of nucleotide strings between DNA and random sequences with the same length. Suppression of nucleotide strings for HUMHBB starts at $`k`$=5. For a random sequence with the same length, which is given by using a random number generator, suppression of nucleotide strings is delayed to start at $`k`$=7. The number of suppressed nucleotide strings for the random number is 5.67% of that for HUMHBB. The finite sequence of HUMHBB effects on the suppression of 7-nucleotide strings. Along with increase of $`k`$, numbers of suppressed nucleotide strings for the random number more increase and approach those for HUMHBB. At $`k=`$10, the number of suppressed nucleotide strings for the random number is 99.3% of that for HUMHBB. In this case, suppression of nucleotide strings in HUMHBB is mainly caused by the finite length of sequence. Moreover, suppression of nucleotide strings for YEAST1 starts at $`k`$=7. For a random sequence with the same length, which is given by using a random number generator, suppression of nucleotide strings is delayed to start at $`k`$=8. The number of suppressed nucleotide strings for the random number is 22.7% of that for YEAST1. The finite sequence of YEAST1 effects on the suppression of 8-nucleotide strings. At $`k`$=10, the number of suppressed nucleotide strings for the random number is 97.5% of that for YEAST1. Due to the comparison of suppression of nucleotide strings, we can thus obtain that HUMHBB and YEAST1 have shorter suppressed nucleotide strings than random sequences with the same lengths, respectively. Along with increase of string lengths, the finite sequences take stronger effects on suppression of nucleotide strings.

In order to display genomic signature, we must eliminate effects of finite sequences on suppression of nucleotide strings. For a DNA sequence, we take the longest string length before suppression of nucleotide strings in a random sequence with the same lengths as an optimal option of string lengths. According to the definition, string lengths 6 and 7 can be chosen as optimal options for genomic signatures of HUMHBB and YEAST1, respectively.

## III Limits of self-similarity scales

Suppression of certain nucleotide strings in the DNA sequences leads to a fractal pattern seen in the MR of DNA sequences. To quantify the fractal feature in the MR of DNA sequences, we introduce information dimension. For a given length $`k`$ of nucleotide strings, we have $`M(=Nk+1)`$ subsequences $`\mathrm{\Sigma }_i(i=k,k+1,\mathrm{},N)`$, which end with $`M`$ $`k`$-nucleotide strings. The subsequences are corresponding to $`M`$ points in a MR. In the MR, the length of a zone and the total number of zones are $`ϵ=3^k`$ and $`Z=4^k`$, respectively. The numbers of points falling in the $`i`$-th zone and of non-empty zones are labeled by $`m_i(ϵ)`$ and $`Z(ϵ)`$, respectively. Dividing the number $`m_i(ϵ)`$ by the total point number $`M`$ yields a probability $`p_i(ϵ)`$ for the $`i`$-th zone. Information function and dimension for the points in MR are respectively defined as

$$I(ϵ)=\underset{i=1}{\overset{Z(ϵ)}{}}p_i\mathrm{log}p_i,$$
(5)
and

$$D_1=\underset{ϵ0}{lim}\frac{I(ϵ)}{\mathrm{log}(1/ϵ)}.$$
(6)
The information function $`I(ϵ)`$ during a range of $`\mathrm{log}(1/ϵ)`$ has a scaling region. The scaling region reflects the self-similarity of pattern in the MR. The information dimension $`D_1`$ can be found from the slope in scaling region $`I(ϵ)`$ versus $`\mathrm{log}(1/ϵ)`$. When the length $`ϵ`$ of a zone in MR increases from $`3^k`$ to $`2^k`$, MR of DNA sequences changes to CGR. Information dimension in CGR can thus be determined as $`(\mathrm{log}_23)D_1`$. We compute information function $`I(ϵ)`$ with different sizes $`ϵ`$ for HUMHBB and draw in Fig. 2. A linear part of the curve $`I(ϵ)`$ versus $`\mathrm{log}(1/ϵ)`$ emerges between $`\mathrm{log}(1/ϵ)=\mathrm{log3}=1.10`$ and $`\mathrm{log}(1/ϵ)=6\mathrm{l}\mathrm{o}\mathrm{g}3`$=6.59. A fitting line is also draw in Fig. 2. The point for $`\mathrm{log}(1/ϵ)=7\mathrm{l}\mathrm{o}\mathrm{g}3`$=7.69 is started leaving from the line. Along with the decrease of $`\mathrm{log}(1/ϵ)`$, farther and farther the points leave from the line. Since points in the zones correspond to $`k`$-nucleotide strings, we can obtain that the self-similarity of pattern in the MR preserves approximately from mononucleotides to $`6`$-nucleotide strings, as well as the suppression of many nucleotide strings emerges at $`7`$-nucleotide strings. Using the least-squares fit method for the liner part, we determine its slope, i.e., information dimension $`D_1`$, to 1.20. It is less than the information dimension 1.26 for random sequence with the same length. Moreover, in Fig. 3, we draw information function $`I(ϵ)`$ versus $`\mathrm{log}(1/ϵ)`$ for YEAST1. A linear part of the curve $`I(ϵ)`$ versus $`\mathrm{log}(1/ϵ)`$ exists between $`\mathrm{log}(1/ϵ)=\mathrm{log3}=1.10`$ and $`\mathrm{log}(1/ϵ)=7\mathrm{l}\mathrm{o}\mathrm{g}3`$=7.69. We can obtain that the suppression of many nucleotide strings in YEAST1 emerges from $`8`$-nucleotide strings. Using the least-squares fit method for the liner part, we also plot a fitting line in Fig. 3 and determine its slope, i.e., information dimension $`D_1`$, to 1.22. It is less than the information dimension 1.26 for random sequence with the same length. The limits of self-similarity in MR of HUMHBB and YEAST1 are equivalent to the optimal string lengths for genomic signatures, respectively. Thus, for presenting genomic signature, a self-similarity limit as an optimal string length can be determined in computing information dimension.

Using the MR method, we determine suppression of $`k`$-nucleotide strings of bacteria complete genomes in Table II, where we put suppression of $`k`$-nucleotide strings in the order of decrease. For each of the bacteria complete genomes, a linear part exists in the plot of information function $`I(ϵ)`$ versus $`\mathrm{log}(1/ϵ)`$. From the linear parts, we determine self-similarity limits of genomic signatures in Table II. Keeping in the order, we find the suppression of bacteria complete genomes does not necessarily depend on the lengths of sequences. The common optimal string length for the bacteria complete genomic signatures can be chosen as 7.

## IV Conclusion

In summary, we have shown MR of DNA sequences is one-to-one. Due to the MR method, suppression of nucleotide strings in the DNA sequences is determined. For a DNA sequence, an optimal string length to display genomic signature is obtained by eliminating effects of the finite sequence. The optimal string length is further shown as a self-similarity limit in computing information dimension. By using the method, self-similarity limits of bacteria complete genomic signatures are further determined.

###### Acknowledgements.

This work was supported in part by the National Key Program for Developing Basic Science G1999032801-11.
FIGURES
Fig. 1 Metric representation of HUMHBB. Its boundary and partition lines are labeled by solid lines and dash lines, respectively.
Fig. 2 A plot of information function $`I(ϵ)`$ versus $`\mathrm{log}(1/ϵ)`$ labeled by dots and its fitting line for HUMHBB.
Fig. 3 A plot of information function $`I(ϵ)`$ versus $`\mathrm{log}(1/ϵ)`$ labeled by dots and its fitting line for YEAST1.


format: text

text:
# Probabilistic prediction in scale-free networks: Diameter changes

\[Phys. Rev. Lett. 91, 058701 (2003)\]

## Abstract

In complex systems, responses to small perturbations are too diverse to predict how much they would be definitely, and then such diverse responses can be predicted in a probabilistic way. Here we study such a problem in scale-free networks, for example, the diameter changes by the deletion of a single vertex for various in silico and real world scale-free networks. We find that the diameter changes are indeed diverse and their distribution exhibits an algebraic decay with an exponent $`\zeta `$ asymptotically. Interestingly, the exponent $`\zeta `$ is robust as $`\zeta 2.2(1)`$ for most scale-free networks, insensitive to the degree exponents $`\gamma `$ as long as $`2<\gamma 3`$. However, there is another type with $`\zeta 1.7(1)`$ and its examples include the Internet and its related in silico model.

A complex system consists of many constituents, generating emerging behaviors through diverse interactions nature ; science . One of the powerful ways of examining the intrinsic nature of a complex system is to observe how such emerging patterns change by small perturbation applied to the system. In complex systems, such a change or response is so sensitive to the details of the perturbation that it is extremely diverse. In such a case, it is not adequate to predict how much the change would be definitely. Recently, Parisi has argued parisi  that the prediction for the responses to small perturbations in complex systems can be made in a probabilistic way. He showed examples of protein structures in biological systems and spin glasses in physical systems. In case of proteins, subject to small external perturbations such as the change in pH or the substitution of a single amino acid, they would fold to a completely different 3D structure but with practically the same free energy. In case of the disordered magnetic systems, each spin responds to a slowly varying external field by changing its orientation, forming a series of bursts, known as Barkhausen noise crackle . The number of spins burst depends on the disorder strength of the system, following a power-law distribution at a critical strength of disorder. The prediction of the number of spins burst in this case can only be probabilistic. The stock market is another example of complex systems. Stock prices are determined as a result of the complicated interplay between numerous investors, and the price changes were also found to exhibit a power-law distribution price . All these examples aptly illustrate how the concept of probabilistic prediction may apply as a new paradigm in modern science. Other examples can also be found in as diverse fields as the meteorology and the geology buchanan .

Recently, there are many works to describe complex systems in terms of graphs rev ; porto , where vertices represent constituents and edges interactions between constituents. An interesting feature emerging in such complex networks is the emergence of a power-law behavior in the degree distribution, $`P_\mathrm{d}(k)k^\gamma `$ ba , where the degree $`k`$ is the number of edges incident upon a given vertex. Such complex networks are called scale-free (SF) networks.

In this Letter, we study how SF networks respond to small perturbations and check if the concept of probabilistic prediction can be applied. For this purpose, we investigate a simple problem of diameter change when a single vertex is removed from the system. Diameter, defined as the average distance between every pair of vertices in a network, is a simple yet fundamental quantity of SF networks to characterize the small-world nature, and can be thought of as a measure reflecting the efficiency of a network. Our main interest is how much the efficiency of a network would be affected by the removal of a single vertex. When a vertex is removed, each pair of remaining vertices whose shortest pathway had passed through the removed vertex should find detours, resulting in the rearrangement of shortest pathways over the network. Thus the diameter change occurs in a collective manner. From the extensive numerical calculations for a number of SF network models and real-world examples, we find that the diameter changes indeed are very diverse and crucially depend on the degree of the removed vertex. When a vertex with a few number of connections is removed, the diameter changes little. However, when a vertex with a large number of connections is removed, the diameter change is drastic, exhibiting a power-law distribution with an exponent $`\zeta `$,

$$P_\mathrm{c}(\mathrm{\Delta })\mathrm{\Delta }^\zeta $$
(1)
for large $`\mathrm{\Delta }`$, where $`\mathrm{\Delta }`$ is the dimensionless relative diameter change defined as the diameter change caused by the removal of a certain vertex divided by the original diameter before the removal, and $`P_\mathrm{c}(\mathrm{\Delta })`$ is its distribution. Moreover the exponent $`\zeta `$ turns out to be robust for various SF networks, insensitive to the degree exponent $`\gamma `$ for $`2<\gamma 3`$.

To be specific, we consider a undirected SF network with finite number of vertices $`N`$ and measure the diameter of the network. Note that we limit our interest to undirected networks only in this work. Next we remove a certain vertex $`i`$ and measure the diameter $`d_i`$ of the rest of the network. Measuring a dimensionless quantity, $`\mathrm{\Delta }_i=(d_id_0)/d_0`$ for all $`i`$, where $`d_0`$ is the diameter of the original unperturbed network, we obtain the distribution of $`\mathrm{\Delta }`$ for the network. Note that our case is different from the previous study of the robustness of SF networks archilles ; lai  where vertices are removed successively. In our case, on the other hand, only a single vertex is removed each time. When a certain vertex is removed, the network may disintegrate into more than one cluster. In such cases, $`d_i`$ is calculated only within the largest cluster. The diameter can be measured via a simple breadth-first search algorithm. To obtain the distribution of the diameter changes, we need the computation time of order $`𝒪(N^3)`$.

To look into more details, we consider the static model load . It is constructed by connecting $`mN`$ pairs of vertices $`i`$ and $`j`$ with probability proportional to $`(ij)^\alpha `$, where $`N`$ is the vertex number and $`\alpha `$ is a parameter. We use $`m=2`$. Its degree distribution follows a power law, $`P_d(k)k^\gamma `$ with $`\gamma =1+1/\alpha `$. Thus, tuning the parameter $`\alpha `$ in $`[0,1)`$, we obtain a continuous spectrum of the exponent $`\gamma `$ in the range $`2<\gamma <\mathrm{}`$.

The diameter change $`\mathrm{\Delta }_i`$ by the removal of a certain vertex $`i`$ in such SF networks can be positive or negative, and the histogram of the diameter changes is highly centralized around $`\mathrm{\Delta }=0`$ (Fig. 1). However, it exhibits a fat tail for $`\mathrm{\Delta }>0`$ (the inset of Fig. 1). For the static model with $`\gamma =3`$, for example, the case of small diameter changes in the range $`|\mathrm{\Delta }|<2\times 10^4`$ occurs with frequency as high as 96%. Thus the effect of a vertex removal usually is negligible as a whole, which is manifested by the exponentially bounded fluctuations of the diameter around its original value. We estimate the $`N`$-dependence of such small diameter changes in a mean-field-type approach. It is known that the diameter $`d_0`$ depends on the number of vertices as $`d_0\mathrm{ln}N`$ for random graph and $`d_0\mathrm{ln}N/\mathrm{ln}\mathrm{ln}N`$ ultra  for the Barabási-Albert model ba  with $`\gamma =3`$. When a vertex is removed, the diameter may be reduced as $`d\mathrm{ln}(N1)`$ or $`d\mathrm{ln}(N1)/\mathrm{ln}\mathrm{ln}(N1)`$, both leading to $`\mathrm{\Delta }1/N\mathrm{ln}N`$ for large $`N`$. Thus when $`N=10^4`$, $`\mathrm{\Delta }𝒪(10^5)`$, which is comparable to numerical values of the central part in Fig. 1. On the other hand, substantial (about 4%) vertices have a serious impact on the system’s efficiency and they indeed contribute to the positive tail of the histogram, showing the power-law behavior, Eq. (1). We find that such large diameter changes are mainly due to the removal of a vertex with large degree. This feature is reminiscent of the percolation problems on the SF networks percol ; havlin .

Let us investigate the power-law behavior for large $`\mathrm{\Delta }`$ in details. The exponent $`\zeta `$ seems to be robust as $`\zeta 2.2(1)`$ as long as $`2<\gamma 3`$ for the static model as shown in Fig. 2. Similar behaviors are found in other model networks (ii)–(vii) listed in Table 1. These include the SF networks showing nontrivial degree-degree correlations assort . For $`\gamma >3`$, on the other hand, as $`\gamma `$ increases, the power-law behavior sets in only for larger values of $`\mathrm{\Delta }`$ and the exponent $`\zeta `$ increases with $`\gamma `$. Eventually the diameter change distribution for the Erdös-Rényi random networks decays exponentially as shown in Fig. 2.

To see such universal behavior of $`\zeta `$ in real world, we consider a couple of real-world networks, the protein interaction networks (PIN) and the Internet. For the PIN of the yeast Saccharomyces cerevisiae goh3 , we also find a power law in the diameter change distribution with an exponent $`\zeta 2.3(1)`$ (Fig. 3), consistent with the one obtained for various model networks, including the one proposed as its own in silico model (vii) sole\_pin .

For the Internet at the autonomous systems level nlanr , the diameter change distribution again follows a power law, however, with a different exponent $`\zeta 1.7(1)`$ (Fig. 4). The smaller exponent $`\zeta `$ indicates that the effect of the removal of vertices contributing to the tail of the distribution is much more severe than the previous cases with $`\zeta 2.2`$ \[(i)–(viii) in Table 1\]. To confirm the novel value of $`\zeta `$ for the Internet, we perform the same calculations for its in silico model, called the adaptation model goh\_internet , and indeed obtain $`\zeta 1.7`$ for it, too. The two different behaviors of the diameter change distribution are rooted from distinct topological features of shortest pathways of each case, which will be discussed later.

Recently, it was proposed that the SF networks with $`2<\gamma 3`$ can be classified into two classes load ; class , following the power-law behavior of the betweenness centrality (BC) distribution freeman ; newman1 . The BC $`g_k`$ of a vertex $`k`$ is the accumulated sum of the fraction of shortest pathways passing through $`k`$ and its distribution follows a power law, $`P_g(g)g^\eta `$ for SF networks. The BC exponent $`\eta `$ turns out to be robust as either $`\eta 2.2(1)`$ (class I) or $`\eta 2.0(1)`$ (class II) as long as $`2<\gamma 3`$ load ; class . Interestingly, the networks (i)–(viii) in Table 1 having the diameter change exponent $`\zeta 2.2`$ belong to the class I, and the values of $`\zeta `$ and $`\eta `$ coincide with each other within our numerical resolutions, while they are different for the class II. Empirically, the rank of a vertex in $`g`$ and that in $`\mathrm{\Delta }`$ are likely to be the same for vertices with large degrees. If then, the relation $`P_g(g)dgP_c(\mathrm{\Delta })d\mathrm{\Delta }`$ would hold asymptotically, leading to

$$\mathrm{\Delta }(g)g^{(\eta 1)/(\zeta 1)},$$
(2)
for large $`g`$. This type of relation also holds between degree and BC load . Indeed, the slopes in the double logarithmic scale in the upper insets of Figs. 3 and 4 are $`1.1(1)`$ for the PIN and $`1.4(1)`$ for the Internet, respectively, consistent with the predictions from the formula, Eq. (2). Thus the two classes, the classes I and II, are also categorized by the diameter change distribution and the distinction between them can be observed more clearly through it.

Our finding that the diameter change distribution is also classified into the classes I and II following those for BC distribution may be rooted from the fact that both quantities, diameter and BC, depend on universal features of the shortest pathways topology between a vertex pair in networks. When the sum rule sumrule , i.e., $`_kg_kd`$ is applied, one can see immediately that the diameter change distribution is the same as the total BC change distribution. On the other hand, the networks belonging to the class II are more sparse and ramified than those in the class I, so that the Internet is more fragile by the removal of a single vertex than the PIN. We compare the distribution of the size change $`\delta S`$ of the largest cluster for the PIN and the Internet by a single vertex removal. As shown in the lower insets of Figs. 3 and 4, the giant cluster in the Internet becomes much smaller than in the PIN. Thus the number of vertex pairs connected after the removal becomes much smaller in the Internet than in the PIN. Consequently, the difference of the exponent $`\zeta `$ between the two classes appears much larger than that of the exponent $`\eta `$ in the class II. However, it is not clear how the power-law behavior in $`P_c(\mathrm{\Delta })`$ arises and what determines its exponent.

It would be interesting to generalize our study to the case of having more than one vertex removed. For simplicity, we consider the case of two-vertex removal, in particular, one is the hub and the other is an arbitrary vertex. Interestingly we find that the diameter change distribution also exhibits a fat-tail behavior with the same exponent $`\zeta `$. We cannot check if the fat-tail behavior still holds for more general cases due to huge amount of computing time. Meanwhile, it has been studied that SF network is robust against random failures. To show this, the diameter change due to removal of a finite fraction of vertices was measured by taking average over a few samples, not over the whole ensemble. If the diameter change distribution after those removals still possesses a power law distribution with $`\zeta 3`$, as is the case here, then we could say that the average diameter change cannot reflect its intrinsic nature because its variance diverges. Thus the diameter changes after random failures should also be described in a probabilistic way.

In summary, we have studied the diverse behavior in response to a small perturbation, a deletion of a single vertex in SF networks. The diameter change $`\mathrm{\Delta }`$ by a removal of a vertex is very diverse, exhibiting a power-law distribution with an exponent $`\zeta `$ for large $`\mathrm{\Delta }`$. Moreover, the diameter change exponent $`\zeta `$ is robust as $`\zeta 2.2`$ for most SF networks with $`2<\gamma 3`$, or $`\zeta 1.7`$ for the Internet as an exception.

###### Acknowledgements.

This work is supported by the KOSEF Grant No. R14-2002-059-01000-0 in the ABRL program.


text:
# Cosmological parameter estimation and Bayesian model comparison using VSA data

## 1 Introduction

In the past two years, a number of experiments have produced accurate measurements of the power spectrum of anisotropies in the cosmic microwave background (CMB) radiation on a range of angular scales (Hanany et al. 2002; Netterfield et al. 2002; Halverson et al. 2002; Sievers et al. 2002; Benoit et al. 2002). These data, together with other cosmological observations, have been used to place increasingly tight constraints on the values of cosmological parameters in current models of the formation and evolution of structure in the Universe.

In this letter, we repeat this process with the inclusion of the latest observations from the Very Small Array (VSA) in its extended configuration. Results from the VSA in its compact configuration have already been presented in Watson et al. (2002), Taylor et al. (2002), Scott et al. (2002) and Rubiño-Martin et al. (2002) (hereafter Papers I - IV). In Grainge et al. (2002, hereafter Paper V) these data are combined with the new extended configuration observations to produce a combined power spectrum with 16 spectral bins spanning the range $`\mathrm{}=1601400`$. This joint set of observed band-powers provides powerful new constraints on cosmological parameters. In this letter we extend the traditional likelihood approach of previous analyses to a fully Bayesian treatment, including careful consideration of our knowledge of cosmological parameters prior to the inclusion of any data, and the calculation of Bayesian evidences to perform model comparisons.

## 2 Models, methods and priors

We restrict our attention to cosmological models in which the initial density fluctuations are adiabatic with a simple power-law spectrum; such perturbations are naturally produced in the standard single-field inflationary model. We assume the contents of the Universe to consist of three components: baryonic matter, dark matter and vacuum energy, with (present day) densities denoted by $`\mathrm{\Omega }_\mathrm{b}`$, $`\mathrm{\Omega }_{\mathrm{dm}}`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ respectively, measured as a fraction of the critical density required to make the Universe spatially-flat (with $`\mathrm{\Omega }_\mathrm{b}+\mathrm{\Omega }_{\mathrm{dm}}+\mathrm{\Omega }_\mathrm{\Lambda }=\mathrm{\Omega }_{\mathrm{tot}}`$).

### 2.1 Model parameterisation and priors

The parameterisation of the cosmological model can be performed in numerous ways, although it is generally preferable to use physically-motivated parameters along principal degeneracy directions. To this end, in the most general case, we describe the cosmological model using the following 11 parameters: the Hubble parameter $`h`$ (defined as $`H_0=h\times 100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$); the physical baryon density $`\omega _\mathrm{b}\mathrm{\Omega }_\mathrm{b}h^2`$; the physical dark matter density $`\omega _{\mathrm{dm}}\mathrm{\Omega }_{\mathrm{dm}}h^2`$; the curvature density $`\mathrm{\Omega }_\mathrm{k}=1\mathrm{\Omega }_{\mathrm{tot}}`$; the fraction $`f_\nu `$ of dark matter in the form of massive neutrinos; the parameter $`w`$ describing the equation of state of the dark energy ($`p=w\rho `$), the redshift of (instantaneous) reionisation $`z_{\mathrm{re}}`$; the amplitude of scalar modes $`A_s`$; the spectral index of scalar modes $`n_s`$; the amplitude ratio $`R`$ of tensor to scalar modes and the spectral index of the tensor modes $`n_t`$.

We consider four separate models of increasing complexity, in which a successively larger number of the above parameters are are allowed to vary. Model A assumes spatial flatness, no massive neutrinos and no tensor modes; Model B includes the possibility of non-zero curvature; Model C additionally allows for the presence of tensor modes; and finally Model D also allows for a contribution to the dark matter in the form of massive neutrinos and also a variable $`w`$. The variable parameters common to all four models are listed in Table 1, together with the top-hat priors assumed for each. The values and priors assumed for the other parameters are shown in Table 2 for each of the four models under consideration.

The basic parameters described above completely define the models considered. It is of interest, however, also to consider the following derived parameters: $`\mathrm{\Omega }_\mathrm{b}`$, $`\mathrm{\Omega }_{\mathrm{dm}}`$, $`\mathrm{\Omega }_\mathrm{m}=\mathrm{\Omega }_\mathrm{b}+\mathrm{\Omega }_{\mathrm{dm}}`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, the age of the universe, the present day rms fluctuation in 8 $`h^1\mathrm{Mpc}`$ spheres as predicted by linear theory $`\sigma _8`$, and the optical depth to the surface of the last scattering $`\tau `$. In addition to the priors listed in Table 1 and 2, we also impose the additional constraints $`\mathrm{\Omega }_\mathrm{\Lambda }0`$ and a top-hat prior on the age of the Universe lying between 10 and 20 Gyr.

### 2.2 Bayesian analysis using MCMC samping

Our approach to Bayesian parameter estimation and model selection makes use of Markov-Chain Monte Carlo (MCMC) sampling to explore the posterior distribution of the cosmological parameters (see e.g. Paper IV). For any given model $`M`$, samples are drawn from the (unnormalised) posterior distribution given by the product of the likelihood and prior, i.e. $`\mathrm{Pr}(\text{d}|𝜽,M)\mathrm{Pr}(𝜽|M)`$, where d denotes the data under analysis and $`𝜽`$ denotes the parameters defining the model. The likelihood function is evaluated in the same manner as discussed in Paper IV, and the joint prior is simply the product of the individual priors discussed in the previous section.

The particular implementation of the MCMC approach used is a slightly adapted version of the Cosmo-mc software package (Lewis & Bridle 2002). This sampler uses CAMB (Lewis, Challinor & Lasenby 2000) as its underlying power spectrum engine, and is specifically tailored to the analysis of CMB data. It achieves high acceptance rates and rapid convergence by using a proposal function that exploits the difference between ‘fast’ and ‘slow’ parameters in CAMB.

Given an ensemble of samples from the posterior probability distributions, we can easily obtain one-dimensional marginalised posterior distributions for any desired parameter, which thus determine the constraints implied by the data. Moreover, MCMC techniques provide a natural way to perform model selection through the evaluation of the Bayesian evidence (see e.g. Hobson, Bridle & Lahav 2002). For a given model $`M`$, the evidence is given by

$$E\mathrm{Pr}(\text{d}|M)=\mathrm{Pr}(\text{d}|𝜽,M)\mathrm{Pr}(𝜽|M)𝑑𝜽,$$

which is simply the average of the likelihood over the prior. The value of the evidence naturally incorporates the spirit of Ockham’s razor: a simpler theory, having a more compact parameter space, will generally have a larger evidence than a more complicated theory, unless the latter is significantly better at explaining the data. Thus, the problem of model selection is answered simply by identifying the model with the largest evidence. We have extended the Cosmo-mc program to include the calculation of the evidence by thermodynamic integration (see e.g. Hobson & McLachlan 2002). This is implemented by numerically evaluating the appropriate integral during the ‘burn-in’ of the Markov chain. ’Burn-in’ samples are discarded for the purpose of the parameter estimation.

The resulting software was run on a 24-node Linux cluster, with each node propagating an independent chain. After burn-in, typically 5000 accepted samples were drawn from each chain, yielding a total of 120,000 accepted samples. Since successive samples from a Markov chain are, by nature, correlated, the accepted samples were thinned by a factor of 15 before analysing them further; this resulted in 8,000 independent samples. We have found that our final results are very insensitive to the thinning factor. Since each chain is run independently, we also obtained 24 separate estimates of the evidence, which are used to obtain an average value for $`E`$ and an estimate of the associated error.

For calculation of confidence limits we use the $`0.165`$, $`0.5`$ and $`0.835`$ points of the cumulative probability distribution. Thus, our parameter estimate is the median of the marginalised posterior pdf and the confidence interval encompasses 67% of the probability. Using the median instead of the maximum of the posterior has two advantages. Firstly, it gives consistent results under monotonic parameter transformation and combination (e.g. the median $`\omega _b`$ is equal to the median $`\mathrm{\Omega }_b`$ multiplied by the square of the median $`h`$). This is especially advantageous in the case of cosmological parameter estimation when it is not always clear which parameters should be considered basic (for more information see Jaynes 2002, p.621). Secondly, the method is very robust when dealing with MCMC chains: the samples may be simply sorted and searched with no need for binning and smoothing. Of course, quoting best estimates and confidence limits is simply a means of characterising the full one-dimensional marginalised distributions for each parameter. Finally, for Models A and B, the parameter estimation presented here has been repeated independently using the standard grid-based method and the results are consistent. We note that the grid based approach would be computationally prohibitive for more complex Models C and D.

### 2.3 Parameter constraints from priors alone

The priors on the individual basic parameters are summarised in Tables 1 and 2. As noted earlier, we also adopt the constraints (on derived parameters) $`\mathrm{\Omega }_\mathrm{\Lambda }>0`$ and the top-hat prior (10,20) Gyr on the age of the Universe. It is of interest to determine the effect of these combined priors alone on the one-dimensional marginalised posteriors for each parameter. Using the MCMC sampler, this is easily performed by simply setting the likelihood to a constant value, i.e. the analysis is performed using no data. The resulting marginalised distributions for Model A and Model B are shown in Fig. 1.

These results are worthy of some discussion. Consider, for example, the distribution for $`h`$. This has a clear peak around $`h0.7`$, which may be understood as follows. Models with low $`h`$ are relatively disfavoured because they would require $`\mathrm{\Omega }_\mathrm{\Lambda }<0`$ even for modest values of $`\omega _\mathrm{m}=\omega _{\mathrm{dm}}+\omega _\mathrm{b}`$, while models with high $`h`$ are disfavoured because they tend to have ages below 10 Gyr. Similar arguments may be applied to explain the shapes of the other distributions. In particular, we note that, although broad, the distributions for $`\mathrm{\Omega }_\mathrm{m}`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ peak at 0.3 and 0.7 respectively. The non-uniformity of the distributions in Fig. 1 represents the effect of the sensible, and seemingly innocuous, initial constraints on the chosen parameters. We have also repeated the analysis for Model D, and find that the resulting marginalised distributions for the extra parameters $`R`$, $`n_t`$, $`f_\nu `$ and $`w`$ accurately follow the initial flat prior. Having made explicit the impact of our priors, we can now modulate the above distributions by including increasing amounts of data from various cosmological observations via the likelihood function.

## 3 Results

We consider three different sets of data in our analysis. Firstly, we use only VSA data (from both extended and compact configurations using the main binning from Paper V), together with the COBE power spectrum points (Smoot et al. 1992) to constrain the low-$`\mathrm{}`$ normalisation. For this data set, we investigated just Models A and B. Secondly, we also include data from several recent CMB experiments, namely Maxima (Hanany et al. 2002), Boomerang (Netterfield et al. 2002), DASI (Halverson et al. 2002), CBI (Sievers et al. 2002) and Archeops (Benoit et al. 2002). For these data, we consider Models A, B and C. Finally, we include a wide range of additional independent cosmological probes to obtain the tightest possible constraints on the model parameters. In addition to CMB observations, we include constraints from: the HST Key Project (Freedman et al. 2001); the first 147,000 redshifts measurements from the 2dF Survey (Colles et al. 2001) on scales $`0.02<k/(h\text{Mpc}^1)<0.15`$ (Lewis and Bridle 2002; Percival et al. 2002); nucleosynthesis (Burles, Nollett & Turner 2001); and type IA supernovae (Perlmutter et al. 1999). We also include the constraint from the gas fraction in relaxed clusters (Allen, Schmidt & Fabian 2002; Allen et al. 2002). We do not use the local X-ray luminosity function constraints from the cited papers due to possible systematic uncertainties. Following Allen et al., we marginalise over the cluster bias parameter. Throughout the analysis we marginalise over applicable calibration and beam uncertainty (see Paper IV).

The resulting one-dimensional marginalised posterior distributions, for each of the three data sets, are shown in Figs 24 for a selection of interesting parameters. The corresponding parameter estimates, confidence limits and evidences are given in Table 3. Note that the evidence values have been scaled so that, for each data set, Model A has an evidence $`E`$ of unity.

## 4 Discussion and Conclusions

By comparing the effective priors on the cosmological parameters plotted in Fig. 1 with the posterior distributions inferred from the VSA and COBE data alone (Fig. 2), we see that the main strength of the VSA lies in its large $`\mathrm{}`$-range, which allows one to constrain a wide variety of cosmological parameters. In particular, we observe that the VSA data significantly improves the constraints on $`\omega _\mathrm{b}`$, $`\omega _{\mathrm{dm}}`$, $`\mathrm{\Omega }_\mathrm{k}`$, $`\sigma _8`$ and $`n_s`$. We note that, as a direct consequence of the pronounced third peak in the VSA power spectrum, the preferred value for $`\omega _\mathrm{b}`$ is somewhat larger than that from the nucleosynthesis constraint (Burles et al. 2001) or that from the combined CMB data. This excess is statistically significant only at the $`1.6\sigma `$ level. VSA and COBE data alone do not, however, provide significant new constraints on either $`h`$ or $`\mathrm{\Omega }_\mathrm{\Lambda }`$, especially for Model B.

On comparing Figs 2 and 3, we see that, as one would expect, the constraints on all parameters are tightened by the inclusion of all the CMB data. Indeed, for Model A, reasonable constraints are obtained on all the parameters under consideration, except for $`z_{\mathrm{re}}`$, although there is some indication that low values of $`z_{\mathrm{re}}`$ are preferred. For many of the parameters, the constraints are not significantly broader when one allows the spatial curvature to vary in Model B. Nevertheless, for this model, some parameters do become relatively unconstrained, in particular $`h`$, $`\mathrm{\Omega }_\mathrm{m}`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$, and the limits on the age of the Universe are widened considerably.

Using all recent available cosmological data sets allows one to place tight constraints on nearly all parameters, even for Model B. Indeed, only $`z_{\mathrm{re}}`$ remains relatively unconstrained, but the earlier indication that lower values are preferred appears to be reinforced. We note that, aside from the parameters $`h`$ and the age of the Universe, the constraints obtained are very similar for Model A and Model B. In fact, from Table 3, we see that, even for Models C and D, the constraints on the parameters plotted in Fig. 4 are not significantly broadened. In Fig. 5, we also plot the constraints obtained on the additional parameters $`R`$, $`n_t`$, $`f_\nu `$ and $`w`$ for Model D. We see that the data favour low values of $`R`$ and high values of $`n_t`$. More significant constraints can be placed on $`f_\nu `$ and $`w`$. We note from the figure that the marginalised distribution for each parameter possesses a single peak, although these are not pronounced features, especially for $`f_\nu `$. Of particular interest is that the preferred value of $`w`$ corresponds to a cosmological constant and that a firm upper limit is obtained on the fraction of dark matter in the form of massive neutrinos.

To distinguish between the different cosmological models considered here, we have calculated the Bayesian evidence for each model-data set combination. It must be remembered that evidence values can only be compared between different models using the same data set. From Table 3, we see that, for each data set, Model A (i.e. the simple 6-parameter $`\mathrm{\Lambda }`$CDM cosmology) is preferred, although the evidence ratio between Model A and Model B is of order unity in each case. For Models C and D, however, the evidence ratio with Model A is considerably larger and shows clearly that such general models are not necessary to explain current cosmological observations.

Since Model A is found to have the largest evidence, it is of interest to consider more closely the impressively tight parameter constraints that can be achieved in this simple case using all the available cosmological data. In particular, let us focus on $`h`$ and $`\omega _\mathrm{b}`$, which are parameters directly probed by the HST Key Project (Freedman et al. 2001) and nucleosynthesis data (Burles et al. 2001). In Fig. 6, we plot the posterior distributions of these parameters for Model A, with (thick line) and without (thin line) the corresponding direct constraint (dashed line). In each case, we see that the direct constraint is overwhelmed by the combination of the other data sets.

It is clear from the above analyses that the constraints on the parameters in the simplest models become impressively small when a wide range of cosmological probes are used. It remains to be seen, however, whether this level of precision is matched by a corresponding level of accuracy in the values determined. The model comparison exercise demonstrated here does not address the possibility that different experiments have different systematic errors. Nevertheless, one would hope that these might be reduced when many experiments are combined. In future work, we will investigate other models, with different combinations of free and fixed parameters and perform a hyperparameter analysis (Hobson et al. 2002) to reveal any discrepancies between data sets.

## ACKNOWLEDGEMENTS

We thank Antony Lewis and Sarah Bridle for providing their Cosmo-mc code and many useful discussions regarding MCMC parameter estimation. We thank Carolina Ödman for useful discussions. We also thank the staff of the Mullard Radio Astronomy Observatory, the Jodrell Bank Observatory and the Teide Observatory for invaluable assistance in the commissioning and operation of the VSA. The VSA is supported by PPARC and the IAC. RS, YH, KL, and PM acknowledge support by PPARC studentships. KC acknowledges a Marie Curie Fellowship. YH is supported by the Space Research Institute of KACST. AS acknowledges the support of St. Johns College, Cambridge. We thank Professor Jasper Wall for assistance and advice throughout the project.


text:
# 1 Introduction

## 1 Introduction

Let us consider on an hypercubic lattice of dimension $`D`$ all closed paths of length $`2k`$, starting from a given vertex. An interesting and still unanswered question concerns the counting of such paths according to the area they enclose. Besides being a challenging combinatorial problem on its own right, the question has some relevance also from a physical point of view: considering for instance the hopping of a charged particle around a two dimensional lattice in a uniform magnetic field , its energy spectrum displays an intricate hierarchic structure , whose properties are closely related to the area versus length distribution of closed paths . Such distribution finds applications also in polymer physics , and in explaining anomalous magnetoconductance . A different but somehow related problem is the one of counting random walks according to the number of distinct sites they visit : relatively compact loops are expected to enclose a small area and to visit relatively few sites, and viceversa. This counting problem is relevant when investigating random walks on lattices with static traps.

Coming back to the initial question of counting loops according to their enclosed area, in the case of a two-dimensional lattice the answer was given only in the limit of asymptotically large $`k`$ in , and the first subleading correction was provided in . In higher dimension, the problem is in general even more difficult. Here we want to address the question in the limit of infinite dimensionality of the lattice, a very peculiar situation, in which however some simplifications occur, and an exact answer can be given. The latter is to be interpreted as a “mean-field” approximation to the exact counting on a lattice of finite dimensionality $`D`$. Such a formulation of the problem was first considered by Parisi et al. . They were investigating spin models with frustration but without any quenched disorder, in order to test the conjecture that such deterministic models could behave at low temperature as some suitably chosen spin-glass model with quenched disorder. They considered the frustrated Spherical and XY spin models in the limit of large dimensionality $`D`$ of the lattice, where the saddle point approximation becomes exact. In their analysis of these models, they showed how the high temperature expansion (i.e. loop expansion) can be nicely rewritten by using the $`q`$-oscillator algebra , where $`q`$ measures the frustration per plaquette. Similar $`q`$-deformed algebraic relations have also appeared in the Hofstadter problem of quantum particles hopping on a two-dimensional lattice in a uniform magnetic field . This problem is closely related to the frustrated spin models, which can be considered as simplified models of hopping in the large $`D`$ and classical limits. The use of somehow related non-commutative geometry techniques appeared as a crucial ingredient also in Ref. .

The problem we address and solve in this paper may be stated as follows: given an infinite dimensional hypercubic lattice, what is the number $`G_{k,l}`$ of loops of length $`2k`$, starting from a given vertex, enclosing a minimal area of exactly $`l`$ plaquettesWhen facing this problem of area versus length distribution of closed random walks, the area is most often considered in its algebraic sense, i.e. as oriented; in our case, however, an average over orientations is intrinsecally performed in the approach (see Section 2), and we are therefore referring to the absolute value of the area.?

The paper is organized as follows. In next Section we summarize some results of  which are then used in Section 3 to obtain an exact expression for the counting numbers $`G_{k,l}`$. In Section 4 a generating function for these counting numbers is presented. In Section 5 a probability distribution is naturally associated to the $`G_{k,l}`$’s, and a systematic procedure is built to compute its moments of arbitrary order. The asymptotic behaviour of this probability distribution for large $`k`$ is evaluated in Section 6, while Section 7 is devoted to the presentation of some new, nontrivial identities. The last Section contains our conclusions, while some technical details are relegated into the Appendix.

## 2 Frustrated Lattice and $`q`$-oscillator Algebra

In Refs. , Parisi et al. considered the frustrated Gaussian, Spherical and XY models on a $`D`$-dimensional hypercubic lattice, in the large $`D`$ limit. Among other things, they unveiled a remarkable connection between the coefficients of the high temperature expansion (i.e., loop expansion) and the $`q`$-oscillator algebra , where $`q`$ measures the frustration per plaquette and varies continuously on the real interval $`[1,+1]`$, between the fully-frustrated case ($`q=1`$, fermionic algebra) and the ferromagnetic case ($`q=1`$, bosonic algebra). For our purposes, the relevant part of the Hamiltonians describing the models considered in Refs.  is:

$$H=\frac{1}{\sqrt{2D}}\underset{jk}{}\varphi _j^{}U_{jk}\varphi _k+\mathrm{h}.\mathrm{c}..$$
(1)
The complex field $`\varphi _j𝐂`$ is defined on the sites (labelled by $`j`$) of a $`D`$-dimensional hyper-cubic lattice; frustration is induced on the lattice through the nearest-neighbour couplings $`U_{jk}`$, which are complex numbers of modulus one and satisfy the relation $`U_{jk}=U_{kj}^{}`$; they are the link variables of a background Abelian lattice gauge field, chosen in such a way as to produce a static and constant magnetic field, suitably oriented to give the same magnetic flux $`\pm B`$ for any plaquette of the lattice (the product of the four $`U`$’s around the plaquette being $`e^{\pm iB}`$). Such a magnetic field having the same projection over all the axes, up to the sign, in Refs.  these signs were chosen randomly in order to avoid the selection of a preferred direction. The usual unfrustrated ferromagnetic spin interaction is obtained for $`B=0`$. Non-vanishing values of $`B`$ induce a frustration around each plaquette, which is maximal for $`B=\pi `$, the fully-frustrated case.

In the framework of the high temperature expansion, the free energy of such models is expressed as a sum over the contribution of loops of increasing length $`2k`$:

$$\beta F=\underset{k=0}{\overset{\mathrm{}}{}}\frac{\beta ^{2k}}{2k}G_k.$$
(2)
Each loop encloses a number of plaquettes; in the case of the models considered in , for each loop the magnetic field yields a weight, proportional to $`\mathrm{exp}(iBA)`$, where $`A`$ is the sum of plaquettes with signs depending on the orientations. The total contribution of all loops of length $`2k`$ is given by $`G_k`$. Due to the average over orientations and loops, the quantity $`G_k`$ is a polynomial in the variable $`q=\mathrm{cos}B`$, the coefficient $`G_{k,l}`$ of $`q^l`$ being given by the number of loops of length $`2k`$ and area $`l`$ The area of a loop is defined as the minimal area of a surface of lattice plaquettes having that loop as boundary.. The infinite dimensionality of the lattice ensures that for any given loop, no two plaquettes of the subtended minimal surface will lie on the same plane; this in turn accounts for such a simple averaging over orientations, performed indipendently for each plane. The order of polynomial $`G_k(q)`$ is given by $`k(k1)/2`$, the maximal area encloseable by a loop of length $`2k`$. The coefficients $`G_{k,l}`$ can also be interpreted as the number of Feynman diagrams with $`2k`$ external points, which are joined pairwise by lines (propagators) intersecting $`l`$ times. These diagrams also occur in the topological (large $`N`$) expansion of Matrix Models , where the planar limit corresponds to no intersections, i.e. to the $`q=0`$ case. Equivalently, in simple graphical terms, they just can be seen as the number of way of connecting pairwise $`2k`$ points on a circle with $`k`$ segments intersecting exactly $`l`$ times. In the following we shall refer to this last picture§§§This formulation of the problem was considered e.g. by Linus Pauling  when trying to simplify the calculation of matrix elements involved in Slater’s treatment of the electronic structure of molecules.. In Ref. , the enumeration of such diagrams was investigated. In particular, a recursion relation was found for the coefficients of the polynomial $`G_k(q)`$ – a sort of Wick theorem – which can be nicely expressed in terms of the algebra of the $`q`$-oscillators $`a_q,a_q^{}`$:

$$a_qa_q^{}qa_q^{}a_q=1.$$
(3)
These operators  act on the Hilbert space spanned by the vectors: $`|m`$, $`m=0,1,\mathrm{}`$, as follows,

$`a_q|m=`$ $`\sqrt{[m]_q}|m1,`$ $`a_q|0=0,`$
$`a_q^{}|m=`$ $`\sqrt{[m+1]_q}|m+1,`$ $`[m]_q={\displaystyle \frac{1q^m}{1q}}.`$ (4)
Using the recursion relation, the weighted multiplicities of the diagrams of eq. (2) were neatly written as an expectation value over the ground state of the $`q`$-oscillators :

$$G_k(q)=0|(a_q^{}+a_q)^{2k}|0.$$
(5)
The authors of Refs.  did not exploit the consequences of this result, preferring to turn themselves to numerical investigation of the models they were interested in. As a consequence, till now only the two limiting cases $`G_k(0)`$ and $`G_k(1)`$ were explicitly known; when $`q=1`$ it is only a matter of counting the way of connecting $`2k`$ points on a circle, with no restrictions, and this is simply the number of pairings of $`2k`$ objects: $`(2k1)!!`$ . When $`q=0`$ we are in fact evaluating the planar limit of the zero-dimensional $`2k`$-point Green function of a Matrix Model in the limit of vanishing interaction . In simple graphical terms this correspond to the number of way of joining pairwise $`2k`$ points on a circle with no intersection. In other words, this is just one of the many possible definition of Catalan numbers (see, e.g. ), given as:

$$G_k(0)=\frac{(2k)!}{k!(k+1)!}.$$
(6)
But the number $`G_{k,l}`$ of way of connecting pairwise $`2k`$ points on a circle, with exactly $`l`$ intersections (i.e. the coefficient of $`q^l`$ in $`G_k(q)`$) is not so easily accessible. In Refs. , they were found by direct enumeration of the graphs on a computer. In , a generating function for coefficients $`G_{k,l}`$ was obtained,proposed, but then the explicit evaluation of such coefficients relied upon heavy symbolic manipulations.

## 3 Exact Solution of the Enumeration Problem

We shall now present a simple, easy to evaluate, formula for a generic coefficient $`G_{k,l}`$. To this purpose let us first introduce the $`x_q`$ coordinate representation, $`x_q=a_q^{}+a_q`$, $`x_q|x=x|x`$, which is given by the so-called continuous $`q`$-Hermite polynomials . These are defined by:

$$H_n(x)=x|n𝒞_n,𝒞_n=\left(\left[n\right]_q!\right)^{1/2}𝒞_0,$$
(7)
where the normalization constant $`𝒞_0`$ is fixed by $`H_0(x)=1`$ and the $`q`$-factorial is

$$\left[n\right]_q!=[n]_q[n1]_q\mathrm{}[1]_q,[1]_q=[0]_q=1.$$
(8)
These polynomials satisfy, of course, a three-term recursion relation in the index $`n`$:

$$xH_n(x)=H_{n+1}(x)+[n]_qH_{n1}(x),n1.$$
(9)
$`x`$ ranges over the interval $`x[2/\sqrt{1q},2/\sqrt{1q}]`$, and a convenient convenient parametrisation is

$$x=\frac{2}{\sqrt{1q}}\mathrm{cos}\theta ,\theta [0,\pi ].$$
(10)
More properties of these $`q`$-Hermite polynomials can be found in Ref. , where they are defined as $`_n(\mathrm{cos}\theta )=(1q)^{n/2}H_n(x)`$. The most important property for us is the orthogonalizing measure $`\nu _q(x)`$ :

$$_{2/\sqrt{1q}}^{2/\sqrt{1q}}\nu _q(x)𝑑xH_n(x)H_m(x)=\delta _{n,m}[n]_q!,$$
(11)
$`\nu _q(x)`$ $`=`$ $`{\displaystyle \frac{\sqrt{1q}}{2\pi }}q^{1/8}\mathrm{\Theta }_1({\displaystyle \frac{\theta }{\pi }},\tau )`$ (12)
$`=`$ $`{\displaystyle \frac{\sqrt{1q}}{\pi }}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}(1)^nq^{n(n+1)/2}\mathrm{sin}[(2n+1)\theta ],`$
where $`\mathrm{\Theta }_1(z,\tau )`$ is the first Jacobi theta function :

$$\mathrm{\Theta }_1(z,\tau )=\underset{n𝐙}{}\mathrm{e}^{i\pi \tau (n+\frac{1}{2})^2+2\pi i(n+\frac{1}{2})(z\frac{1}{2})},q=\mathrm{e}^{2\pi i\tau }.$$
(13)
It is worth emphasizing that, interestingly enough, measure (12) interpolates continuously between Wigner semi-circle law for Gaussian Matrix Models, and Gaussian distribution, as $`q`$ varies from $`0`$ to $`1^{}`$; this reminds us of the behaviour of the density distribution of eigenvalues for Gaussian Ensemble of Random Matrices, as the dimension of the matrices, $`N`$, varies from $`\mathrm{}`$ to $`1`$.

Polynomial $`G_k(q)`$ may now be rewritten as

$$G_k(q)=_{2/\sqrt{1q}}^{2/\sqrt{1q}}\nu _q(x)x^{2k}𝑑x.$$
(14)
Indeed, using explicit form (12) for the integration measure, performing the integral and playing a little bit with indices, we may write

$$G_k(q)=\left(\frac{1}{1q}\right)^k\underset{l=0}{\overset{k}{}}(1)^l\left(\genfrac{}{}{0pt}{}{2k+1}{kl}\right)\frac{2l+1}{2k+1}q^{l(l+1)/2}.$$
(15)
We now only need to perform explicitly the division of the polynomial of degree $`k(k+1)/2`$ defined by the sum in the previous formula. To this purpose, let us state the following

TheoremWe are unfortunately unable to quote any reference. We are of course convinced this theorem has been enounced long time ago, indeed it is just a corollary of Ruffini’s rule. The only proof we have been able to give, absolutely inelegant and unworthy of being presented here, is by direct verification.: Let $`P(q)=_{l=0}^Np_lq^l`$ be an integer coefficient polynomial of degree $`N`$ in $`q`$, exactly divisible by $`(1q)^k`$. Then $`A(q)=P(q)/(1q)^k`$ is an integer coefficient polynomial of degree $`Nk`$ whose coefficients are simply expressed in terms of the $`p_l`$’s as follows:

$$A(q)=\underset{l=0}{\overset{Nk}{}}a_lq^l,a_l=\underset{i=0}{\overset{l}{}}\left(\genfrac{}{}{0pt}{}{k+l1i}{k1}\right)p_i.$$
(16)
We therefore readily get a simple closed expression for $`G_{k,l}`$:

$$G_{k,l}=\underset{i=0}{\overset{i_{max}}{}}(1)^i\left(\genfrac{}{}{0pt}{}{k+l1i(i+1)/2}{k1}\right)\left(\genfrac{}{}{0pt}{}{2k+1}{ki}\right)\frac{2i+1}{2k+1},l\frac{k(k1)}{2},$$
(17)
where $`i_{max}=\left[(1+\sqrt{1+8l})/2\right]_{i.p.}`$ is the largest integer $`i`$ satisfying $`i(i+1)/2l`$. Because of their definition as counting numbers, all integers $`G_{k,l}`$ should be positive. This fact is indeed not apparent from (17), and we are not able to present a rigorous proof of this statement, nevertheless we are convinced of its validity. Sensibleness arguments rely upon the intrinsic positivity of (5) and its derivatives with respect to $`q`$, for any value of $`q`$ in the interval $`[0,1)`$.

We have computed explicitly with Mathematica , through formula (17), which is very efficient, all values $`G_{k,l}`$ for $`k`$ ranging from $`1`$ to $`9`$, and verified that they indeed match the results of Ref.  which were found by direct enumeration of the graphs on a computer.

## 4 Generating Function

Let us define a generating function for the combinatorial numbers $`G_{k,l}`$:

$$R(z,q)\underset{k=0}{\overset{\mathrm{}}{}}G_k(q)z^k=\underset{k,l=0}{\overset{\mathrm{}}{}}G_{k,l}q^lz^k.$$
(18)
The running of index $`l`$ may or may not be restricted to the range $`0lk(k1)/2`$, without loss of generality: coefficients $`G_{k,l}`$ automatically vanish for $`l>k(k1)/2`$. From expression (14), we may readily write a generating functions for coefficients $`G_k(q)`$:

$$R(z,q)=_{2/\sqrt{1q}}^{2/\sqrt{1q}}\nu _q(x)\frac{1}{1zx^2}𝑑x.$$
(19)
Inserting for the measure its explicit expression (12), and performing the integral, we obtain:

$$R(z,q)=\sqrt{\frac{1q}{z}}\underset{n=0}{\overset{\mathrm{}}{}}(1)^nq^{n(n+1)/2}\left[K\left(\frac{z}{1q}\right)\right]^{2n+1},$$
(20)
where we introduced

$$K(\alpha )=\frac{1\sqrt{14\alpha }}{2\sqrt{\alpha }}.$$
(21)
For $`q=0`$ function (20) reduces to the standard generating function of Catalan numbers. The singularities of $`R(z,q)`$ in the complex $`z`$-plane are completely determined by those of the simpler function $`K(z/(1q))`$. Actually the series is very well convergent for $`|q|<1`$ as any Jacobi theta function. Moreover, for $`|q|=1`$, it is a geometric series, which is still convergent because $`|K(\alpha )|<1`$ for $`|z|<z_c(1q)/4`$.

Generating function (20), as shown in Ref.  in a different context, can be interpreted physically as the internal energy at temperature $`T=1/\sqrt{z}`$ of a Gaussian spin model defined on the sites of the infinite-dimensional frustrated lattice introduced in Ref. .

## 5 Moments of the Distribution

The interpretation of coefficients $`G_{k,l}`$ naturally suggests to introduce, for each positive integer value of index $`k`$, the following normalized (discrete) probability distribution:

$$P_k(l)=\frac{G_{k,l}}{_{l=0}^{k(k1)/2}G_{k,l}}=\frac{G_{k,l}}{G_k(1)}=\frac{G_{k,l}}{(2k1)!!},\{\begin{array}{c}k=0,1,\mathrm{}\hfill \\ l=0,1,\mathrm{},k(k1)/2,\hfill \end{array}$$
(22)
whose interpretation is obvious: given an arbitrary diagram with $`2k`$ legs, the probability for it to have exactly $`l`$ crossings is simply given by $`P_k(l)`$; equivalently, $`P_k(l)`$ measures the probability, for a randomly chosen closed path of length $`2k`$ on our infinite dimensional hypercubic lattice, to have area equal to $`l`$. The moments of the previous probability distribution are defined as

$$M_s^{(k)}l^s_k=\underset{l=0}{\overset{k(k1)/2}{}}l^sP_k(l),$$
(23)
and can be computed from the following generating function

$$\stackrel{~}{P}_k(q)q^l_k=\frac{G_k(q)}{G_k(1)}.$$
(24)
Indeed,

$$M_s^{(k)}=\left[\left(q\frac{}{q}\right)^s\stackrel{~}{P}_k(q)\right]_{q=1}=\left[\left(\frac{}{\mathrm{log}q}\right)^s\stackrel{~}{P}_k(q)\right]_{q=1}.$$
(25)
The moments of the deviations with respect to the mean $`M_1^{(k)}`$ (central moments) are given by a similar expression:

$$m_s^{(k)}(lM_1^{(k)})^s_k=\left[\left(\frac{}{\mathrm{log}q}\right)^s\stackrel{~}{P}_k(q)q^{M_1^{(k)}}\right]_{q=1}.$$
(26)
The computation of such moments, besides being interesting on its own right, is (at least for the first two) a necessary step towards the evaluation of the asymptotic form of probability distribution (22) in the limit $`k\mathrm{}`$. We therefore need a systematic procedure to calculate $`G_k(q)`$, together with its derivatives, at $`q=1`$. Indeed, from expression (14), we essentially need to evaluate the derivatives at $`q=1`$ of the integration measure (12). Unfortunately such a limit, due to the presence of the Jacobi Theta Function, is highly singular: for fixed values of $`x`$, $`\nu _q(x)`$ develops an essential singularity at $`q=1`$. A prescription should be given, in our situation the natural one being to chose to approach $`q=1`$ from real, smaller than $`1`$, values of $`q`$; in a more convenient, equivalent notation, we want to evaluate the behaviour of $`\nu _q(x)`$ in the limit $`ϵ0^+`$, $`ϵ`$ real, defined by $`q\mathrm{e}^ϵ`$. With such a choice, as shown in the Appendix, we can write any given moment $`M_s^{(k)}`$ as a simple combination of trivial Gaussian moments:

$$M_s^{(k)}=\frac{(1)^s}{(2k1)!!}\frac{1}{\sqrt{2\pi }}_{\mathrm{}}^+\mathrm{}𝑑xx^{2k}c_s(x)\mathrm{e}^{\frac{x^2}{2}},$$
(27)
where $`c_s(x)`$ is an even polynomial of degree $`4s`$ in $`x`$. The only practical difficulty is the explicit computation of the $`s`$’th coefficient $`c_s(x)`$ of the Taylor expansion around $`ϵ=0^+`$, of the function

$$\sqrt{1\mathrm{e}^ϵ}\mathrm{e}^{\frac{ϵ}{8}}\frac{1}{\sqrt{ϵ}}\mathrm{e}^{\frac{1}{2}\left[\frac{4}{ϵ}\mathrm{arcsin}^2\left(\frac{x}{2}\sqrt{1\mathrm{e}^ϵ}\right)x^2\right]}.$$
(28)
The expansion can be done by hand up to the first few orders in $`ϵ`$, but soon becomes intractable. It is however a trivial task if we use some symbolic manipulation program such as Mathematica. We shall give here explicit results for the first few moments:

$`M_0^{(k)}`$ $`=`$ $`1`$
$`M_1^{(k)}`$ $`=`$ $`k(k1)/6`$
$`M_2^{(k)}`$ $`=`$ $`k(k1)(5k^2k+12)/180`$ (29)
$`M_3^{(k)}`$ $`=`$ $`k(k1)(35k^4+14k^3+235k^2188k+24)/7560`$
$`M_4^{(k)}`$ $`=`$ $`k(k1)(175k^6+315k^5+2341k^41959k^3+2056k^25664k2160)/226800`$
Analogously, the evaluation of the central moments $`m_s^{(k)}`$, is reduced to the explicit computation of the $`s`$’th coefficient of the Taylor expansion around $`ϵ=0^+`$, of the function

$$\sqrt{1\mathrm{e}^ϵ}\mathrm{e}^{\frac{ϵ}{8}}\frac{1}{\sqrt{ϵ}}\mathrm{e}^{\frac{1}{2}\left[\frac{4}{ϵ}\mathrm{arcsin}^2\left(\frac{x}{2}\sqrt{1\mathrm{e}^ϵ}\right)x^2\right]}\mathrm{e}^{ϵk(k1)/6},$$
(30)
with the following results:

$`m_0^{(k)}`$ $`=`$ $`1`$
$`m_1^{(k)}`$ $`=`$ $`0`$
$`m_2^{(k)}`$ $`=`$ $`(k+3)k(k1)/45`$ (31)
$`m_3^{(k)}`$ $`=`$ $`(2k+3)(2k+1)k(k1)/945`$
$`m_4^{(k)}`$ $`=`$ $`k(k1)(7k^4+37k^3+7k^2108k45)/4725`$
By direct inspection of the formulae, it is possible to extract the large $`k`$ behaviour of generic moments. Since $`c_s(x)`$ is a polynomial of degree $`4s`$ in $`x`$, $`M_s^{(k)}(2k+4s1)!!/(2k1)!!=O(k^{2s})`$ as $`k\mathrm{}`$More precisely, it is easy to show that $`M_s^{(k)}=(k^2/6)^s[1+O(1/k)]`$ for large $`k`$.. Following a similar line of thought, and taking into account some cancellations which may occur more or less severely, according to the parity of $`s`$, it is also possible to guess that $`m_s^{(k)}=O(k^{[\frac{3}{2}s]_{i.p.}})`$ for asymptotically large $`k`$.

Let us further notice that, on the infinite dimensional hypercubic lattice, from $`M_1^{(k)}k^2/6`$, the area of random loops increases on average as the square of their length. This should be contrasted to the two-dimensional lattice , where the average area increases of course as the square of the average linear extension of the loop, that is as the first power of the length itself.

We shall conclude this Section by considering coefficients $`G_{k,l}`$ for fixed values of $`l`$. We then have for each integer value $`l`$ an infinite succession of integer numbers. From direct inspection of formula (17), for asymptotically large values of $`k`$ we can write:

$$G_{k,l}=G_{k,0}\frac{k^l}{l!}\left[1+O\left(\frac{1}{k}\right)\right],l\mathrm{fixed}.$$
(32)
## 6 Asymptotic Behaviour of Distribution $`P_k(l)`$

Let us turn back to probability distribution $`P_k(l)`$, eq. (22). We want to investigate its asymptotic behaviour as $`k\mathrm{}`$. The standard procedure is to rescale the discrete index $`l`$ to a continuous variable $`t`$ using the average and the standard deviation as computed in Section 5 for generic $`k`$:

$$l\mu _k+\sigma _kt,\genfrac{}{}{0pt}{}{\mu _kM_1^{(k)}=\frac{k(k1)}{6};}{\sigma _k^2m_2^{(k)}=\frac{(k+3)k(k1)}{45}.}$$
(33)
It is now possible to define the asymptotic probability distribution (22) as:

$$P(t)=\underset{k\mathrm{}}{lim}\sigma _kP_k(\mu _k+\sigma _kt).$$
(34)
Unfortunately the coefficients $`G_{k,l}`$, are defined in terms of alternate sums of very large numbers (about $`\sqrt{2l}`$ terms, each one of order $`O(k^k)`$, summing up to a much smaller number \**\**\**Just as an example, let us quote that $`G_{100,1000}`$ is expressed through formula (17) as an alternate sum of 45 numbers of order up to $`10^{200}`$, being itself of order $`10^{179}`$ i.e. 21 orders of magnitude smaller.), so that the evaluation of the asymptotic distribution $`P(t)`$ is not completely straightforward. Indeed, we are interested into the behaviour of the $`G_{k,l}`$’s for large $`k`$ and large $`l`$, $`lk^2/6`$, and the simple estimate of eq. (32) is now completely useless.

In principle the asymptotic behaviour of coefficients $`G_k(q)`$ for large value of $`k`$ is encoded in the behaviour of generating function $`R(z,q)`$, eq. (20), in the vicinity of its closest (with respect to the origin) singularity in the variable $`z`$. The latter is situated at $`z=z_c(1q)/4`$. But since we are in fact interested into the behaviour of $`G_{k,l}`$ where $`l`$ is also asymptotically large, but in a controlled way, see (33), we must inspect the generating function in the vicinity of its closest (again, with respect to the origin) singularity in variable $`q`$, that is near $`q=1`$. But in this way $`z_c`$ is collapsing towards the origin, conspiring to the construction of an essential singularity in $`z=z_c=0`$ as $`q1^{}`$. To tackle this difficulty we shall not investigate the generating function $`R(z,q)`$, but we shall rather resort to the Laplace transform method, applied directly to probability distribution $`P_k(l)`$; let us first introduce the Laplace transform

$$\varphi _k(s)=\underset{l=0}{\overset{k(k1)/2}{}}P_k(l)\mathrm{e}^{s\left(\frac{l\mu _k}{\sigma _k}\right)}$$
(35)
of probability distribution $`P_k(l)`$, suitably shifted and rescaled according to eq. (33). In the large $`k`$ limit, $`\varphi _k(s)`$ is expected to converge pointwise to the Laplace transform of distribution $`P(t)`$:

$$_{\mathrm{}}^+\mathrm{}𝑑tP(t)\mathrm{e}^{st}.$$
(36)
Since we may write

$$\varphi _k(s)=\mathrm{e}^{s\frac{\mu _k}{\sigma _k}}\mathrm{e}^{s\frac{l}{\sigma _k}}=\mathrm{e}^{s\frac{\mu _k}{\sigma _k}}\frac{G_k(\mathrm{e}^{\frac{s}{\sigma _k}})}{G_k(1)},$$
(37)
and $`\sigma _kk^{3/2}`$ for large $`k`$, the asymptotic behaviour of $`P_k(l)`$ is ruled by the behaviour of $`G_k(q)`$ in the neighbourhood of unity. Indeed, after identification of $`ϵ`$ with $`s/\sigma _k`$, we can make use of the machinery developed in Section 5 and Appendix A. We start by using the fact that, up to terms vanishing exponentially for small $`ϵ`$, ($`q=\mathrm{e}^ϵ`$)

$$G_k(q)=_{\mathrm{}}^+\mathrm{}𝑑xx^{2k}(x,ϵ),$$
(38)
with $`(x,ϵ)`$ given by equation (52). For $`ϵ=0`$ ($`q=1`$), the evaluation of the integral reduces to a straightforward application of the saddle-point method: we have to sum up the contribution of the neighbourhoods of two stationary points situated in $`x=x_\pm =\pm \sqrt{2k}`$, respectively. We readily get

$$G_k(1)=\sqrt{2}(2k)^k\mathrm{e}^k\left[1+O\left(\frac{1}{k}\right)\right],$$
(39)
which is indeed the leading asymptotics for $`(2k1)!!`$, as it should. Now, for $`ϵ`$ small but positive (we are therefore specializing to positive real values of $`s`$), for each of the two stationary point $`x=x_\pm `$, we just shift the integration variable $`x`$ to $`x_\pm +y`$ and expand the integrand, including $`(x,ϵ)`$, around $`k=\mathrm{}`$. Considering e.g. the contribution coming from the positive stationary point $`x_+`$, and taking into account only terms up to $`O(\mathrm{log}k)`$, we get

$$I_+\frac{1}{\sqrt{2\pi }}𝑑y\mathrm{e}^{\left[k\mathrm{log}2kky^2\frac{s}{2}\sqrt{5k}\sqrt{10}sy2s^2\right]},(s>0).$$
(40)
Of course the situation in which $`s<0`$ has to be computed, too, but the result turns out to be independent of the sign of $`s`$. Integrating, and adding the analogous contribution coming from the region in the neighbourhood of $`x_{}`$, we finally obtain

$$G_k(\mathrm{e}^{\frac{s}{\sigma _k}})=\sqrt{2}(2k)^k\mathrm{e}^k\mathrm{e}^{\frac{s^2}{2}\frac{s}{2}\sqrt{5k}}\left[1+O(\frac{1}{\sqrt{k}})\right],$$
(41)
and thus, for asymptotically large values of $`k`$

$$\varphi _k(s)=\mathrm{e}^{\frac{s^2}{2}}\left[1+O(\frac{1}{\sqrt{k}})\right].$$
(42)
Indeed, by comparison with eq. (36), we have thus proven that $`P(t)`$ as defined in equation (34) is a standard Gaussian distribution:

$$P(t)=\frac{1}{\sqrt{2\pi }}\mathrm{e}^{\frac{t^2}{2}}.$$
(43)
This result gives us further information about the asymptotic behaviour of the central moments $`m_s^{(k)}`$ of distribution $`P_k(l)`$: indeed, rescaling each central moment $`m_s^{(k)}`$ by a factor $`\sigma _k^sk^{\frac{3}{2}s}`$ (correspondingly to the rescaling of the discrete variable $`l`$ to the continuous one $`t`$), and performing the $`k\mathrm{}`$ limit, we must find the moments of a standard (unit variance) Gaussian distribution, namely $`0`$ or $`(s1)!!`$, for odd or even values of $`s`$, respectively. This on one hand confirms rigorously the estimate $`m_s^{(k)}=\sigma _k^s(s1)!![1+O(1/k)]`$ ($`s`$ even) and on the other hand matches the guessed behaviour $`m_s^{(k)}k^{[\frac{3}{2}s]_{i.p.}}`$ ($`s`$ odd).

As for the area versus length distribution, the result we have found, eq. (43), does not agree with the behaviour proposed in Refs. ; however the role of the $`D\mathrm{}`$ limit, inherent to our treatment, must be taken into proper account. In this limit, all loops with more than one plaquette of the subtended surface lying over the same plane are neglected; on the other hand, on the two-dimensional lattice considered in Refs.  there is only one possible plane. Moreover, having averaged over orientations, we restrict ourselves to the absolute value of area, while in , algebraic (oriented) area is considered.

When considering very compact loops, their distribution is ruled by the left tail of eq. (43). In Ref.  the different but somehow related problem of the distribution $`p(s,t)`$ of random walks of $`t`$ steps visiting exactly $`s`$ different sites was addressed on a cubic lattice. In particular, an asymptotic formula in the regime of large $`t`$ and relatively small $`s`$ was proposed, which describes the distribution of compact walks. However the two distributions are not comparable: indeed, in our situation of infinite dimensionality, the most compact loops of length $`2k`$ have zero enclosed area, but still visit a relatively large number of distinct sites, at least $`k+1`$; more precisely, when the dimensionality of the lattice is very large, the number of walks which take more than one step in the same direction is of order $`1/D`$ with respect to those stepping in $`k`$ different directions, and therefore even the most compact walk, the one coming back to the origin every two steps, at leading order in $`D`$ visits $`k`$ different sites, in addition to the origin. In other words, the connection between the two notions of compactness, as small enclosed area, or as small number of visited sites, is lost in the limit of infinite dimensionality.

## 7 An Infinite Hierarchy of Nontrivial Identities

We would like to conclude with a simple consideration, which in our opinion is however rather striking and deserve further investigation. The exact evaluation of moments $`M_s^{(k)}`$ gives as a byproduct an infinite set of nontrivial and somehow intriguing identities, the simplest of them beeing:

$$\underset{l=0}{\overset{k(k1)/2}{}}G_{k,l}=(2k1)!!k=0,1,2,\mathrm{}$$
(44)
and in general

$$\underset{l=0}{\overset{k(k1)/2}{}}G_{k,l}l^s=M_s^{(k)}(2k1)!!k=0,1,2,\mathrm{}$$
(45)
with $`G_{k,l}`$ and $`M_s^{(k)}`$ explicitely given by formulae (17) and (27-28), respectively. These identities are of course encoded in the Jacobi Theta function, essentially in its modular invariance properties, eq. (48), but they come nevertheless rather unexpected.

Other sets of identities can be obtained if we are somehow able to compute the $`G_{k,l}`$ in a indipendent and simple way. This is indeed very easy e.g. for $`G_{k,k(k1)/2}`$, which counts the number of diagrams with $`2k`$ external points and maximal crossing number, and is thus obviously equal to $`1`$ for any $`k`$ (just connect the $`j^{th}`$ point with the $`(j+k)^{th}`$ one, $`j=1,\mathrm{},k`$). From eq. (17) we may therefore write:

$$\underset{i=0}{\overset{k1}{}}(1)^i\left(\genfrac{}{}{0pt}{}{k(k+1)/21i(i+1)/2}{k1}\right)\left(\genfrac{}{}{0pt}{}{2k+1}{ki}\right)\frac{2i+1}{2k+1}=1k=0,1,2,\mathrm{}$$
(46)
Identities of this second sort are however more obvious: eq. (46) is just the simplest example of a whole bunch of identities encoded in an obvious simmetry of polynomial $`G_k(q)`$: by applying the Theorem of Section 3 (eq. (16)) to the new polynomial $`\stackrel{~}{G}_k(q)=q^{k(k1)/2}G_k(1/q)`$, with $`G_k(q)`$ expressed in terms of eq. (15) one obtains new coefficients $`\stackrel{~}{G}_{k,l}`$, which obviously satisfy $`\stackrel{~}{G}_{k,l}=G_{k,(k(k1)/2)l}`$; setting in particular $`l=0`$, we readily get eq. (46).

## 8 Conclusions

The problem of counting closed paths of given length on a lattice, according to the area they enclose, is a difficult one. In the present paper we have tackled the question on an hypercubic lattice, in the limit of infinite dimensionality. Moreover an average over orientations is inherent to the method we have used. But these conditions have allowed us to give a closed answer, which is moreover exact under the restriction we have assumed. The main flaw of the whole method resides in our opinion into the fact that the computation of even only the first corrections in $`1/D`$ to coefficients $`G_{k,l}`$ is out of reach. On the other hand, we have presented an explicit solution to a combinatorial problem which could hardly be addressed directly; in this respect, the connection between the abstract combinatorial puzzle and the physical model proposed in  is crucial. The exact solution proposed for coefficients $`G_{k,l}`$ has been investigated in detail; their asymptotic behaviour has been evaluated; finally, as a byproduct, an infinite set of new nontrivial identities has been obtained.

Acknowledgements

We would like to thank Andrea Cappelli and Stefano Ruffo for fruitful discussions, Vladimir Rittenberg for suggesting us the use of the Laplace Transform method (Section 6), and Giuseppe Santoro for informing us of Ref.  and its connection with our work.

Appendix A

We want to build up a procedure to compute in a systematic way derivatives of arbitrary order of $`G_k(q)`$ with respect to $`\mathrm{log}q`$, at $`q=1`$. From the explicit expression (14) of $`G_k(q)`$, we have to take into account both the $`q`$ dependence of the integration interval and of the integration measure $`\nu _q(x)`$, given by eq. (12). As for the latter, the task is indeed non completely straightforward, since for fixed values of $`x`$, $`\nu _q(x)`$ develops an essential singularity in the limit $`q1`$, and should somehow be regularized. In presence of an essential singularity, the result of the limiting procedure depends from the direction chosen to approach the singularity, and a regularization is a prescription on how the singularity is approached. Indeed, if we restrict $`q`$ to the real axis and approach $`q=1`$ from smaller values, it is possible to write the integration measure $`\nu _q(x)`$ ($`x=\frac{2}{\sqrt{1q}}\mathrm{cos}\theta `$) as a sum of a completely regular contribution, Taylor expandable in the real variable $`ϵ`$, ($`q=\mathrm{e}^ϵ`$) plus a singular contribution vanishing with all its derivatives for $`ϵ0^+`$. Indeed we shall show that:

$$\nu _q(x)=(x,ϵ)\left[1+O(\mathrm{e}^{1/ϵ})\right],$$
(47)
$`(x,ϵ)`$ being an analytic function of $`ϵ`$.

Using modular invariance, we can write

$$\mathrm{\Theta }_1\left(\frac{\theta }{\pi }|\tau \right)=\frac{i}{\sqrt{i\tau }}\mathrm{e}^{i\theta ^2/\pi \tau }\mathrm{\Theta }_1\left(\frac{\theta }{\pi \tau }|\frac{1}{\tau }\right);$$
(48)
on the other hand, from definition (13), letting $`\tau =\frac{iϵ}{2\pi }`$,

$$\mathrm{\Theta }_1\left(\frac{\theta }{\pi \tau }|\frac{1}{\tau }\right)=\frac{1}{i}\mathrm{e}^{\frac{\pi ^2}{2ϵ}}\mathrm{e}^{\frac{2\pi \theta }{ϵ}}\underset{n𝐙}{}(1)^n\mathrm{e}^{\frac{2\pi ^2}{ϵ}(n^2+n2n\frac{\theta }{\pi })}.$$
(49)
For arbitrary small positive $`ϵ`$, and for $`\frac{\theta }{\pi }`$ in the open interval $`(0,1)`$, the leading contribution is given by the $`n=0`$ term in the infinite sum, all other terms being exponentially small. We may therefore write, for $`\tau =\frac{iϵ}{2\pi }`$, and $`0<ϵ1`$:

$$\mathrm{\Theta }_1\left(\frac{\theta }{\pi \tau }|\frac{1}{\tau }\right)\frac{1}{i}\mathrm{e}^{\frac{\pi ^2}{2ϵ}}\mathrm{e}^{\frac{2\pi \theta }{ϵ}}\left[1+O(\mathrm{e}^{\frac{1}{ϵ}})\right],$$
(50)
and for the measure $`\nu _q(x)`$ we readily get, in the limit $`q1^{}`$, or $`q=\mathrm{e}^ϵ`$, $`ϵ0^+`$, the approximate behaviour (47), with

$`(x,ϵ)`$ $`=`$ $`\sqrt{1\mathrm{e}^ϵ}\mathrm{e}^{\frac{ϵ}{8}}{\displaystyle \frac{1}{\sqrt{2\pi ϵ}}}\mathrm{exp}\left[{\displaystyle \frac{2}{ϵ}}\left(\theta {\displaystyle \frac{\pi }{2}}\right)^2\right]`$ (51)
$`=`$ $`\sqrt{1\mathrm{e}^ϵ}\mathrm{e}^{\frac{ϵ}{8}}{\displaystyle \frac{1}{\sqrt{2\pi ϵ}}}\mathrm{exp}\left[{\displaystyle \frac{2}{ϵ}}\mathrm{arcsin}^2\left({\displaystyle \frac{x}{2}}\sqrt{1\mathrm{e}^ϵ}\right)\right]`$ (52)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}\mathrm{e}^{\frac{x^2}{2}}{\displaystyle \underset{r=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{c_r(x)}{r!}}ϵ^r,`$ (53)
where the coefficients $`c_r(x)`$ are given by the Taylor expansion in the neighbourhood of the origin $`ϵ=0`$ of the analytical function:

$$\sqrt{1\mathrm{e}^ϵ}\mathrm{e}^{\frac{ϵ}{8}}\frac{1}{\sqrt{ϵ}}\mathrm{e}^{\frac{1}{2}\left[\frac{4}{ϵ}\mathrm{arcsin}^2\left(\frac{x}{2}\sqrt{1\mathrm{e}^ϵ}\right)x^2\right]};$$
(54)
$`c_0(x)=1`$ trivially while higher coefficients $`c_r(x)`$ can be shown to be even polynomials of degree $`4r`$ in $`x`$. In the neighbourhood of $`q=1`$ we may therefore write, up to terms vanishing exponentially for small $`ϵ`$:

$$G_k(q)=\frac{1}{\sqrt{2\pi }}_{\frac{2}{\sqrt{1q}}}^{\frac{2}{\sqrt{1q}}}𝑑xx^{2k}\mathrm{e}^{\frac{x^2}{2}}\underset{r=0}{\overset{\mathrm{}}{}}\frac{c_r(x)}{r!}ϵ^r.$$
(55)
Let us now observe that the integration limits diverge as $`1/\sqrt{ϵ}`$, while the integrand falls down as a quadratic exponential. Therefore, up to terms which vanish exponentially for small $`ϵ`$, the integration interval can readily be extended to the whole real axis, and the computation of moments (25) is finally reduced to the Taylor expansion of (54) and the evaluation of a sum of standard Gaussian integrals:

$$M_s^{(k)}=\frac{(1)^s}{(2k1)!!}\frac{1}{\sqrt{2\pi }}_{\mathrm{}}^+\mathrm{}𝑑x\mathrm{e}^{\frac{x^2}{2}}x^{2k}c_s(x).$$
(56)
As an example, let us compute explicitly the first moment, $`M_1^{(k)}`$. The first coefficient of the Taylor expansion of (54) is

$$c_1(x)=\frac{1}{8}+\frac{1}{4}x^2\frac{1}{24}x^4,$$
(57)
and therefore, denoting with $`\mathrm{}`$ the Gaussian average (with variance equal to $`1`$), and exploiting the standard relation $`x^{2k}=(2k1)!!`$

$`M_1(k)`$ $`=`$ $`{\displaystyle \frac{1}{x^{2k}}}\left[{\displaystyle \frac{1}{8}}x^{2k}{\displaystyle \frac{1}{4}}x^{2k+2}+{\displaystyle \frac{1}{24}}x^{2k+4}\right]`$ (58)
$`=`$ $`{\displaystyle \frac{1}{8}}{\displaystyle \frac{1}{4}}(2k+1)+{\displaystyle \frac{1}{24}}(2k+3)(2k+1)`$ (59)
$`=`$ $`{\displaystyle \frac{k(k1)}{6}}.`$ (60)