text:
# Scaling exponents and clustering coefficients of a growing random network

## I Introduction

To study the statistical property of a complex system composed of many interacting individual components, it is often helpful to map the system into a network of nodes and links (edges). Each node in this network represents one component of the real system and the interaction, if there is any, between two components is denoted by an edge, either directed or undirected, between the two corresponding nodes in the network. One quantity of interest is the node-degree profile of the formed network: How does $`n(k)`$, the total number of nodes with a given number $`k`$ of links attached (the node degree), scales with $`k`$? Empirical observations revealed that many social and biological networks have the scale-free property Barabási and Albert (1999); Albert and Barabási (2002), that is,

$$n(k)k^\nu $$
(1)
as $`k`$ becomes large enough. The scaling exponent $`\nu `$ is typically in the range of $`2<\nu <3`$; but there are evidences that some networks have scaling exponents in the range of $`1`$-$`2`$ note1  while a few other networks have scaling exponents larger than $`3`$ (for a collection of experimental data, please refer to Table II of Albert and Barabási (2002)).

To explain the scale-free characteristic, one appealing mechanism is to assume that the network (1) keeps growing and, (2) during this growth process, new edges are generated and are attached preferentially to those nodes which have already been attached by a large number of edges Barabási and Albert (1999). Based on this mechanism several models have been suggested Barabási and Albert (1999); Dorogovtsev et al. (2000); Krapivsky et al. (2000, 2001), but they predicted the scaling exponent $`\nu `$ to be greater than $`2`$, thus failed to explain the behavior of those networks with smaller $`\nu `$. In Ref. Adamic and Huberman (2000) this $`\mathrm{`}\mathrm{`}`$preferential attachment” mechanism was questioned partly because of this apparent discrepancy between theory and empirical data. In Ref. Dorogovtsev and Mendes (2001) it was shown that if one assumes a network is growing accelerately, it is possible to generate scaling exponent in the range between $`3/2`$ and $`2`$. However, it is still not clear whether or not this condition is absolutely necessary to explain experimental observations.

An emerging property of almost all the so-far studied scale-free networks is that they can at the same time be classified as small-world networks Watts and Strogatz (1998). That is, (i) the diameter of the network scales as $`\mathrm{ln}(N)`$, where $`N`$ is the total number of nodes in the system, and (ii) the clustering coefficient $`C`$ is independent of $`N`$ and is thus much greater than that of a random network ($`k/N`$, where $`k`$ denotes the average node-degree of the network). On the theoretical side, it was confirmed that growing networks generated by the mechanism of preferential attachment will typically have diameters scales as $`\mathrm{ln}(N)`$ (see, for example, Ref. Rodgers and Darby-Dowman (2001)). However, the original Barabási-Albert model Barabási and Albert (1999); Albert and Barabási (2002) predicted a very small clustering coefficient $`CN^{0.75}`$. The clustering coefficients for other models Dorogovtsev et al. (2000); Krapivsky et al. (2000, 2001) were not reported.

To improve our understanding on scale-free networks, in this work two questions are addressed: Will it be possible to generate a scale-free network with scaling exponent $`\nu <2`$ based on the preferential attachment mechanism? Will it be possible for a scale-free network generated this way to have relatively constant clustering coefficients? We answer these questions confirmatively by studying a revised Krapivsky-Rodgers-Redner model Krapivsky et al. (2001) with Monte Carlo simulation approach.

## II The growing network model

The Krapivsky-Rodgers-Redner model Krapivsky et al. (2001) is a generalized version of the original Barabási-Albert model on scale-free networks Barabási and Albert (1999). It has the following key elements: (i) edges are directed; (ii) new nodes are added into the network and are attached preferentially to existing nodes with larger in-degrees; (iii) creation of edges between $`\mathrm{`}\mathrm{`}`$old” nodes are possible and a newly created edge also prefer to attach to nodes with larger degrees. This model is general in the sense that it takes into account directional interactions of the real network systems, and that the growth of the network is not solely caused by the inclusion of new nodes but also as a result of the increased interactions among the existing nodes of the system. It can be corresponded to the real systems including the World-Wide Web (WWW), the Internet, the food web, the transportation network, the e-mail network, etc.. Because of its generality, we re-examine this model in the present work to illustrate the property of growing scale-free network.

We noticed in the original Krapivsky-Rodgers-Redner model Krapivsky et al. (2001) the growth process permits the following two possibilities illustrated in Fig. 1: (1) a directed edge can originate from and end into the same node (self-connection) and (2) there can be more than one edges of the same direction between two nodes (multiconnection). Allowing these two possibilities makes analytical calculations possible and these authors found that in the growing network, both the in-degree and the out-degree distribution follow the power-law Krapivsky et al. (2001):

$$n_{\mathrm{in}}(k)k^{\nu _{\mathrm{in}}},\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}2}<\nu _{\mathrm{in}}=2+p\lambda <\mathrm{},$$
(2)
$$n_{\mathrm{out}}(k)k^{\nu _{\mathrm{out}}},\mathrm{\hspace{0.33em}\hspace{0.33em}2}<\nu _{\mathrm{out}}=1+1/(1p)+\mu p/(1p)<\mathrm{},$$
(3)
($`p`$, $`\lambda `$ and $`\mu `$ are the three adjustable model parameters). Although the permission of self-connections and multiconnections might be reasonable for some networks (such as WWW pages, in which a page can have several links to another page and it can have a link bring the reader from one portion to another portion of the same page), it may fail for other kinds of networks (in a co-authorship network Newman (2001), it is meaningless for a author to be the co-author of himself/herself; and in a food web Montoya and Solé (2002), there is at most one edge of the same direction between two species). Because of the preferential attachment mechanism, if a node already has a large number of incoming and outgoing edges, it has a good possibility to form self-connections and by doing so its dominance is further amplified; allowing multiconnection has similar effects. As a result, edges may concentrate on just few nodes, making the scaling behavior steeper. In the present work, the two kind of edges listed in Fig. 1 are discounted. Because of the reasoning outlined above, scaling exponents $`1<\nu <2`$ might occur in a growing network without self-connection and multiconnection. This point will be checked by Monte Carlo simulation.

## III Monte Carlo Set-up

The revised Krapivsky-Rodgers-Redner model is studied by MC simulation. Started with a single node, at each step:

With probability $`p`$, a new node is created and a directed edge from it to an existing target node $`\beta `$ is setup. Of all the $`N`$ existing nodes, $`\beta `$ will be selected with probability Krapivsky et al. (2001)

$$P_{\mathrm{attach}}(\beta )=\frac{k_{\mathrm{in}}(\beta )+\lambda }{E+\lambda N}.$$
(4)
In the above equation, $`k_{\mathrm{in}}(\beta )`$, the in-degree, is the total number of incoming edges of node $`\beta `$; $`\lambda `$ is a constant signifying the $`\mathrm{`}\mathrm{`}`$initial attractiveness” of a node Dorogovtsev et al. (2000); and $`E`$ is the total number of edges in the system before this new edge is created.

With probability $`q=1p`$, a new edge pointing from one node $`\alpha `$ to another node $`\beta `$ is created, provided that (1) $`E<N(N1)`$, (2) $`\alpha `$ and $`\beta `$ are not identical, and (3) there is no any preexisting directed edge from $`\alpha `$ to $`\beta `$. The probability that $`\alpha `$ and $`\beta `$ will be selected is governed by the probability

$$P_{\mathrm{connect}}(\alpha ,\beta )=\frac{[k_{\mathrm{out}}(\alpha )+\mu ][k_{\mathrm{in}}(\beta )+\lambda ]}{\underset{\gamma }{}\underset{\delta }{}^{}[k_{\mathrm{out}}(\gamma )+\mu ][k_{\mathrm{in}}(\delta )+\lambda ]},$$
(5)
where $`k_{\mathrm{out}}(\alpha )`$ denotes the out-degree of node $`\alpha `$; $`\mu `$ is another constant with similar physical meaning as $`\lambda `$; $`\underset{\delta }{}^{}`$ denotes the summation over all the nodes $`\delta `$ which is not yet approached by a directed edge from node $`\gamma `$.

For large system size $`N`$ it turns out to be quite inefficient and complicated when performing procedure (ii) based on a direct application of Eq. (5). This is partly because of the fact that, after a new edge has been created between node $`\alpha `$ and $`\beta `$, one must update the value of the summation in Eq. (5) by $`O(N)`$ iterations. To speed up procedure (ii), the selection of two nodes and the connection of an edge between them is finished actually through the following way:

Select a outgoing node $`\alpha `$ with probability $`[k_{\mathrm{out}}(\alpha )+\mu ]/(E+\mu N)`$;

Select an in-coming node $`\beta `$ with probability $`[k_{\mathrm{in}}(\beta )+\lambda ]/(E+\lambda N)`$;

If $`\alpha `$ and $`\beta `$ are identical, or if there is already a directed edge from $`\alpha `$ to $`\beta `$, repeat steps (1) and (2); if else, accept $`\alpha `$ and $`\beta `$ and update the system.

It is not difficult to prove that by this method the probability for nodes $`\alpha `$ and $`\beta `$ be chosen is identical to Eq. (5).

The algorithm code is written in C++ language Stroustrup (1997), with some of its standard containers (including map and set) being exploited.

To estimate the scaling exponents from the simulated data, we use two methods. One can directly fit the data with Eq. 1. Alternatively, one can define the cumulative degree distribution by

$$P(k)=\underset{k^{}k}{}n(k^{})k^{(\nu 1)},$$
(6)
and from the cumulative distribution data an estimation of the value of $`\nu `$ could be obtained.

## IV  Degree distribution of the growing network

In the growing network model, there are three adjustable parameters, namely $`p`$, $`\lambda `$, and $`\mu `$. Figure 2 shows the relations between the average number of nodes and the in-degree and out-degree for both the original and the revised Krapivsky-Rodgers-Redner model. Figure 3 shows the corresponding cumulative degree distributions for the two models. The network is the result of $`N=10^6`$ growing steps.

At $`p=0.133334`$ (a new node will be included on average every $`7.5`$ steps), $`\lambda =0.75`$ and $`\mu =3.55`$ the original model Krapivsky et al. (2001) predicts $`\nu _{\mathrm{in}}=2.1`$ and $`\nu _{\mathrm{out}}=2.7`$. From the MC data we obtain that $`\nu _{\mathrm{in}}=2.066\pm 0.014`$ and $`\nu _{\mathrm{out}}=2.626\pm 0.036`$, in close agreement with the analytical values. At these same parameters, the revised model has $`\nu _{\mathrm{in}}=1.925\pm 0.007`$ and $`\nu _{\mathrm{out}}=2.269`$. Thus, excluding of self-connection and multiconnection leads to decreased values for the scaling exponents. Other quantitative differences are:

In the revised model there is a cut-off in the in-degree distribution: no node has in-degree $`k>9\times 10^3`$. While in the original model, there are nodes with in-degree as large as $`k=4\times 10^5`$.

In the revised model the out-degree distribution has an exponential cut-off around $`k=250`$; while in the original model, there are nodes with out-degree as large as $`k=1.510^4`$.

The value of $`n(k)`$ is much larger in the revised model than in the original model for a given $`k`$ (less than the cut-off value). This holds both for the in-degree distribution and for the out-degree distribution.

These observations lead to the following picture: By prohibiting self-connection and multiconnection, edges which originally belong to just few $`\mathrm{`}\mathrm{`}`$suppernodes” are now redistributed to those nodes of small or intermediate in- and out-degrees. Consequently, the number of nodes with small and intermediate node-degrees increases considerably, resulting in a smaller scaling exponent in the power-law decrease of the distribution and a cut-off in the tail of this distribution.

In Fig. 4 we demonstrate the simulation result when the probability of node addition is changed to $`p=0.05`$ (a new node will be included on average every $`20`$ steps) while the other two parameters are kept the same values as in Fig. 2. In these parameters, the original model predicts $`\nu _{\mathrm{in}}=2.04`$ (theory) and $`2.018\pm 0.015`$ (MC) and $`\nu _{\mathrm{out}}=2.24`$ (theory) and $`2.190\pm 0.014`$ (MC); while the revised model has an in-degree exponent $`\nu _{\mathrm{in}}=1.672\pm 0.003`$ and an out-degree exponent $`\nu _{\mathrm{out}}=1.764`$, both of which are markedly smaller than $`2`$.

Therefore, the exclusion of self-connection and multi-connection can change the scaling exponent of the scale-free network dramatically when each node has a relatively large average node-degrees. It can be anticipated that similar behavior will be observed when the initial attractiveness parameters, $`\lambda `$ and $`\mu `$, of each node are varied. It is therefore possible for the present model to explain networks with scaling exponent $`\nu <2`$.

However, could the observation of scaling exponents $`1<\nu <2`$ be an artifact caused by finite-size effects? In Fig. 5 the calculated in-degree scaling exponent is plotted as a function of the total growth steps. This figure strongly indicates that even for an infinite system the scaling exponents will still be less than $`2`$.

Adamic and Huberman Adamic and Huberman (2000) studied the WWW by mapping each web domain (rather than each web-page) as a single node, and they reported an in-degree scaling exponent of $`\nu _{\mathrm{in}}=1.94`$. Mossa and co-authors, upon their reinterpretation of the WWW data of Barabási et al. Barabási and Albert (1999), reported an exponent of $`\nu =1.25`$ Mossa et al. (2002). The e-mail network studied in Ref. Ebel et al.  has a scaling exponent in the range $`1.47<\nu <1.82`$. And even smaller scaling exponents are reported in several other networks Albert and Barabási (2002). To attain quantitative agreement to these empirical data, one need nevertheless more information to help fixing the values of the adjustable parameters.

A persistent property of the revised network model is that there is a rapid decay in the out-degree distribution (which occurs at $`k_{\mathrm{out}}400`$ in Fig. 4B). Such a rapid decay was not observed in the original model. This feature is also absent in the in-degree distribution of both the revised and the original model, although the in-degree distribution of the revised model does have a cut-off for large $`k`$. Experimentally, it was reported that both the WWW network Huberman and Adamic (1999); Mossa et al. (2002) and the e-mail network Ebel et al.  show exponential cut-off in the node-degree distribution.

## V  Cluster coefficients of the growing network

As was mentioned in the introduction, many real scale-free network at the same time show small-world behavior Albert and Barabási (2002): having small diameters and being highly clustered. For the original Krapivsky-Rodgers-Redner model, it has been reported in Ref. Rodgers and Darby-Dowman (2001) that the average minimum path scales as $`\mathrm{ln}(N)`$. We anticipate this to be hold also for the revised model. Here we focus on the clustering characteristics of the revised model system, with the parameters setting equal to those of Fig. 2, i.e., $`p=0.133334`$, $`\lambda =0.75`$ and $`\mu =3.55`$, and a total $`10^6`$ growing steps.

### V.1 Mutual-connection coefficient $`C_{\mathrm{mutual}}`$

Denote $`G_{\mathrm{down}}(\alpha )`$ as the complete set of nodes which are the $`\mathrm{`}\mathrm{`}`$downstream” neighbors of node $`\alpha `$, namely there exists a directed edge from $`\alpha `$ to each nodes in $`G_{\mathrm{down}}(\alpha )`$; similarly we define $`G_{\mathrm{up}}(\alpha )`$ as the complete set of nodes which are the $`\mathrm{`}\mathrm{`}`$upstream” neighbors of node $`\alpha `$. $`|G_{\mathrm{down}}(\alpha )|`$ means the size of set $`G_{\mathrm{down}}(\alpha )`$.

The mutual-connection coefficient is defined as

$$C_{\mathrm{mutual}}=\frac{\underset{\alpha }{}|G_{\mathrm{down}}(\alpha )G_{\mathrm{up}}(\alpha )|/|G_{\mathrm{down}}(\alpha )|}{N},$$
(7)
which signifies to what extent the $`\mathrm{`}\mathrm{`}`$downstream neighbors” of one node intersect with its $`\mathrm{`}\mathrm{`}`$upstream neighbors”. The value of $`C_{\mathrm{mutual}}`$ averaged over $`20`$ realization of the growing network is $`0.0010`$. This indicates that the interaction between one node and its $`\mathrm{`}\mathrm{`}`$neighbors” in the network is usually not bi-directional. However, this value is still much larger than the value for a random network of the same size ($`N=133271`$) and the same average degree of out-going edges ($`k=7.50`$), for which $`C_{\mathrm{mutual}}=5.63\times 10^5`$.

### V.2 Incoming clustering coefficient $`C_{\mathrm{in}}`$

Suppose a given node $`\alpha `$ has in-degree $`k_{\mathrm{in}}(\alpha )`$. The maximal number of edges existing between the nodes in $`G_{\mathrm{up}}(\alpha )`$ is $`k_{\mathrm{in}}(\alpha )(k_{\mathrm{in}}(\alpha )1)`$. Denote $`i_{\mathrm{actual}}(\alpha )`$ as the actual number of edges existing between these edges. We define the incoming clustering coefficient as

$$C_{\mathrm{in}}=\frac{\underset{\alpha }{}^{}i_{\mathrm{actual}}(\alpha )/[k_{\mathrm{in}}(\alpha )(k_{\mathrm{in}}(\alpha )1)]}{N^{}}$$
(8)
where $`\underset{\alpha }{}^{}`$ indicates summation over all the nodes whose incoming edges is larger than $`1`$, and $`N^{}`$ is the total number of nodes with this property. We find that $`C_{\mathrm{in}}=0.0044`$. This value indicates that the degree of cliqueness of the upstream neighbors of a given node is usually very small. For a completely random graph, $`C_{\mathrm{in}}=5.63\times 10^5`$.

### V.3  Outgoing clustering coefficient $`C_{\mathrm{out}}`$

The definition of the out-going clustering coefficient $`C_{\mathrm{out}}`$ is similar to that of $`C_{\mathrm{in}}`$. Suppose a particular node $`\alpha `$ has $`k_{\mathrm{out}}(\alpha )`$ out-going edges, and $`i_{\mathrm{actual}}(\alpha )`$ is the total number of edges between the nodes in $`G_{\mathrm{down}}(\alpha )`$, then

$$C_{\mathrm{out}}=\frac{\underset{\alpha }{}^{}i_{\mathrm{actual}}(\alpha )/[k_{\mathrm{out}}(\alpha )(k_{\mathrm{out}}(\alpha )1)]}{N^{\prime \prime }},$$
(9)
where $`N^{\prime \prime }`$ is the total number of nodes whose out-degree is larger than $`1`$. We find $`C_{\mathrm{out}}=0.229`$ for the present growing network. Compared with the small values of $`C_{\mathrm{mutual}}`$ and $`C_{\mathrm{in}}`$, such a large value of $`C_{\mathrm{out}}`$ is surprising. It suggests the average interaction between the nodes which are in the downstream group of a given node is considerably strong. How to understand this kind of asymmetry, namely $`1C_{\mathrm{out}}C_{\mathrm{in}}`$? We suggest the following possibility: In the network, there are some nodes which are so popular that a large population of the whole nodes will have an edge pointing to them (see Fig. 2A). Consequently, these nodes will have great possibility to belong to the downstream group of any particular node, and they will also have great possibility to be pointed to by other members of this group, making $`C_{\mathrm{out}}`$ to be proportional to unity. However, the number of nodes decays quickly when the out-degree increases to about $`250`$ (see Fig. 2B). Therefore, in the network there is no node which are so “generous” that it points to a large population of the whole network. This may make the value of $`C_{\mathrm{in}}`$ small. In other words, it might be the existence of an steep cut-off in the out-degree profile that accounts for the difference in the clustering coefficients $`C_{\mathrm{out}}`$ and $`C_{\mathrm{in}}`$.

### V.4 Triangle coefficient $`C_{\mathrm{triangle}}`$

For a particular node $`\alpha `$, suppose node $`\beta G_{\mathrm{down}}(\alpha )`$. Then $`i_{\mathrm{}}(\alpha ,\beta )=|G_{\mathrm{down}}(\alpha )G_{\mathrm{down}}(\beta )|`$ is the total number of nodes that are pointed to by both $`\alpha `$ and $`\beta `$. We define

$$C_{\mathrm{triangle}}=\frac{\underset{\alpha }{}\underset{\beta G_{\mathrm{down}}(\alpha )}{}i_{\mathrm{}}(\alpha ,\beta )/k_{\mathrm{out}}(\beta )}{N}.$$
(10)
The triangle coefficient $`C_{\mathrm{triangle}}`$ signifies to what extent, if there is a directed edge from node $`\alpha `$ to node $`\beta `$ and there is a directed edge from $`\beta `$ to node $`\gamma `$, there will also be a directed node from $`\alpha `$ to $`\gamma `$. Calculation revealed that $`C_{\mathrm{triangle}}0.011`$.

## VI  Conclusion and Discussion

In this work we have used a revised Krapivsky-Rodgers-Redner model to investigate the degree distribution of growing random network and to investigate whether such a kind of growing network could be regarded as a small-world network. After excluding the possibility of self-connection (Fig. 1A) and requiring that there is at most one directed edge from one node to any another node (Fig. 1B), the Monte Carlo simulation demonstrated that scale-free network with degree distribution coefficient $`\nu `$ less than $`2`$ can be generated. And it is also revealed that the average interactions between the nodes which are belong to the downstream group of a particular node is very strong, suggesting the growing network at the same time forms a $`\mathrm{`}\mathrm{`}`$small world”. The strong interaction in the downstream group of a particular node was suggested to be closely related to the existence of several $`\mathrm{`}\mathrm{`}`$popular” nodes which are pointed to by a large fraction of the total population in the node system. Previous efforts often predicted that the scaling exponent $`\nu `$ should be larger than $`2`$, and there is still not many efforts to understand why many scale-free networks are at the same time small-world networks. It is hoped that the present work will help to improve our understanding of the occur of scale-free networks with $`\nu <2`$ and to improve our understanding of the close relationship between scale-free and small-world networks.

The present work suggests that, by excluding the possibility of self-connection and multiconnection, those many edges, which were associated with several nodes of extremely large in- or out-degrees in the original Krapivsky-Rodgers-Redner network, are now redistributed to nodes of small or intermediate degrees. This may explain why a dramatic decrease in scaling exponents could be observed.

## VII Acknowledgement

This work was financially supported by the Alexander von Humboldt foundation. The author is grateful for the hospitality of Professor Reinhard Lipowsky.


format: text

text:
# Statistical mechanical aspects of joint source-channel coding

## Abstract

An MN-Gallager Code over Galois fields, $`q`$, based on the Dynamical Block Posterior probabilities (DBP) for messages with a given set of autocorrelations is presented with the following main results: (a) for a binary symmetric channel the threshold, $`f_c`$, is extrapolated for infinite messages using the scaling relation for the median convergence time, $`t_{med}1/(f_cf)`$; (b) a degradation in the threshold is observed as the correlations are enhanced; (c) for a given set of autocorrelations the performance is enhanced as $`q`$ is increased; (d) the efficiency of the DBP joint source-channel coding is slightly better than the standard gzip compression method; (e) for a given entropy, the performance of the DBP algorithm is a function of the decay of the correlation function over large distances.

With the rapid growth of information content in today’s wire and wireless communication, there is an increasing demand for efficient transmission systems. A significant gain in the transmission performance can be achieved by the application of the joint source-channel coding technique, which has attracted much attention during the recent past, see for instance . Roughly speaking, source coding is mainly a data compression process that aims at removing as much redundancy as possible from the source signal, whereas channel coding is the process of intelligent redundancy insertion so as to be robust against channel noise. These two processes, source coding and channel coding, seem to act in opposition, where the first/second process shrinks/expands the transmitted data. For illustration, assume that our compression shrinks the size of the source signal by a factor $`2`$ and in order to be robust against channel noise we have to expend our compressed file by a factor $`4`$. Hence, the length of the transmitted sequence is only twice the length of the uncompressed source.

The source-channel coding theorem of Shannon indicates that if the minimal achievable source coding rate of a given source is below the capacity of the channel, then the source can be reliably transmitted through the channel, assuming an infinite source sequence. This theorem implies that source coding and channel coding can be treated separately without any loss of overall performance, hence they are fundamentally separable. Practically, the source can be first efficiently compressed and then an efficient error correction method can be used.

The objective of joint source-channel coding is to combine both source (compression) and channel (error correction) coding into one mechanism in order to reduce the overall complexity of the communication while maintaining satisfactory performance. Another possible advantage of the joint source-channel coding is the reduction of the sensitivity to a bit error in a compressed message.

In a recent paper a particular scheme based on a statistical mechanical approach for the implementation of the joint source-channel coding was presented and the main steps are briefly summarized below. The original boolean source is first mapped to a binary source $`\left\{x_i\pm 1\right\}i=1,\mathrm{},L`$, and is characterized by a finite set, $`k_0`$, of autocorrelations

$$C_k=\frac{1}{L}\underset{i=1}{\overset{L}{}}x_ix_{\left(i+k\right)\mathrm{𝐦𝐨𝐝}L}$$
(1)
where $`kk_0`$ is the highest autocorrelation taken. The number of sequences oeying these $`k_0`$ constraints is given by

$$\mathrm{\Omega }=Tr\underset{k=1}{\overset{k_0}{}}\delta \left(\underset{j}{}x_jx_{j+k}C_kL\right)$$
(2)
Introducing the Fourier representation of the delta functions and then using the standard transfer matrix (of size $`2^{k_0}\times 2^{k_0}`$) method, one finds $`\mathrm{\Omega }=𝑑y_k\mathrm{exp}\{L[y_kC_k\mathrm{ln}\lambda _{max}(\{y_k\})]\}`$, where $`\lambda _{max}`$ is the maximal eigenvalue of the corresponding transfer matrix. For large $`L`$, using the saddle point method the entropy, $`H_2(\{C_k\})`$, is given in the leading order by

$$H_2\left(\{C_k\}\right)=\frac{\frac{1}{k_0}\mathrm{ln}\lambda _{max}\left(\{y_k\}\right)_{i=1}^{k_0}y_kC_k}{\mathrm{ln}2}$$
(3)
where $`\{y_k\}`$ are determined from the saddle point equations of $`\mathrm{\Omega }`$. Assuming Binary Symmetric Channel (BSC) and using Shannon’s lower bound, the channel capacity of sequences with a given set of $`k_0`$ autocorrelations is given by

$$C=\frac{1H_2\left(f\right)}{H_2(\{C_k\})H_2\left(P_b\right)}$$
(4)
where $`f`$ is the channel bit error rate and $`p_b`$ is a bit error rate. The saddle point solutions derived from eq. 3 indicate that the equilibrium properties of the one dimensional Ising spin system

$$H=\underset{i}{}\underset{k=1}{\overset{k_0}{}}\frac{y_k}{\beta }x_ix_{i+k}$$
(5)
obey in the leading order the autocorrelation constraints of eq. 2. Note that in the typical scenario of statistical mechanics, one of the main goals is to calculate the partition function and the equilibrium properties of a given Hamiltonian. Here we present a prescription of how to solve the reverse question. Given the desired macroscopic properties, the set of the autocorrelations, the goal is to find the appropriate Hamiltonian obeying these macroscopic constraints. This property of the effective Hamiltonian, eq. 5, is used in simulations below to generate an ensemble of signals (source messages) with the desired set of autocorrelations.

The decoding of symbols of $`k_0`$ successive bits is based on the standard message passing introduced for the MN decoder over Galois fields with $`q=2^{k_0}`$ and with the following modification. The horizontal pass is left unchanged, but a dynamical set of probabilities assigned for each block are used in the vertical pass. The Dynamical Block Probabilities (DBP), $`\{P_n^c\}`$, are determined following the current belief regarding the neighboring blocks and are given by

$`\gamma _n^c`$ $`=`$ $`S_I\left(c\right)\left({\displaystyle \underset{l=1}{\overset{q}{}}}q_L^lS_L(l,c)\right)\left({\displaystyle \underset{r=1}{\overset{q}{}}}q_R^rS_R(c,r)\right)`$ (6)
$`P_n^c`$ $`=`$ $`{\displaystyle \frac{\gamma _n^c}{_{j=1}^q\gamma _n^j}}`$ (7)
where $`l/r/c`$ denotes the state of the left/right/center ($`n1/n+1/n`$) block respectively and $`q_L^l/q_R^r`$ are their posterior probabilities. $`S_I(c)=e^{\beta H_I}`$, where $`H_I`$ is the inner energy of a block of $`k_0`$ spins at a state $`c`$, see eq. 5. Similarly $`S_L(l,c)`$ ($`S_R(c,r)`$) stands for the Gibbs factor of consecutive Left/Center (Center/Right) blocks at a state $`l,c`$ $`(c,r)`$.

Note that the complexity of the calculation of the block prior probabilities is $`O(Lq^2/\mathrm{log}q)`$ where $`L/\mathrm{log}q`$ is the number of blocks. The decoder complexity per iteration of the MN codes over a finite field $`q`$ can be reduced to order $`O(Lqu)`$, where $`u`$ stands for the average number of checks per block. Hence the total complexity of the DBP decoder is of the order of $`O(Lqu+Lq^2/\mathrm{log}q)`$.

In this Letter we examine the efficiency of the DBP-MN decoder as a function of the maximal correlation length taken, $`k_0`$, the strength of the correlations, the size of the finite fields $`q`$, and we compare this efficiency with the standard gzip compression procedure. A direct answer to the above questions is to implement exhaustive simulations on increasing message length, various finite fields $`q`$, and sets of autocorrelations, which result in the bit error probability versus the flip rate $`f`$. Besides the enormous computational time required , the conclusions would be controversial since it is unclear how to compare, for instance, the performance as a function of $`q`$; with the same number of transmitted blocks or with the same number of transmitted bits.

In order to overcome these difficulties, for a given MN-DBP code over GF(q) and a set of autocorrelations, the threshold $`f_c`$ is estimated from the scaling argument of the convergence time, which was previously observed for $`q=2`$. The median number of message passing steps, $`t_{med}`$, necessary for the convergence of the MN-DBP algorithm is assumed to diverge as the level of noise approaches $`f_c`$ from below. More precisely, we found that the scaling for the divergence of $`t_{med}`$ is independent of $`q`$ and is consistent with

$$t_{med}=\frac{A}{f_cf}$$
(8)
where for a given set of autocorrelations and $`q`$, $`A`$ is a constant. Moreover, for a given set of autocorrelations and a finite field $`q`$, the extrapolated threshold $`f_c`$ is independent of $`L`$, as demonstrated in Fig. 1. This observation is essential to determine the threshold of a code based on the above scaling behavior. Note that the estimation of $`t_{med}`$ is a simple computational task in comparison with the estimation of low bit error probabilities for large $`L`$, especially close to the threshold. We also note that the analysis is based on $`t_{med}`$ instead of the average convergence time, $`t_{av}`$, since we wish to prevent the dramatic effect of a small fraction of finite samples with slow convergence or no convergence.

All simulation results presented below are derived for rate $`1/3`$ and the construction of the matrices $`A`$ and $`B`$ of the MN code are taken from . In all examined sets of autocorrelations, $`10^3L5\times 10^4`$ and $`4q64`$, the scaling for the median convergence time was indeed recovered. At this stage an efficient tool to estimate the threshold of an MN-DBP decoder exists and we are ready to examine the efficiency of the DBP decoder as a function of $`\{C_k\}`$ and $`q`$.

Results of simulations for $`q=4,8,16`$ and $`32`$ and selected sets of autocorrelations are summarized in Table I, and the definition of the symbols is: $`\{C_k\}`$ stand for the imposed values of autocorrelations as defined in eqs. 1-2; $`\{y_k\}`$ are the interaction strengths, eqs. 3 and 5; $`H`$ represents the entropy of sequences with the given set of autocorrelations, eq. 2; $`f_c`$ is the estimated threshold of the MN-DBP decoder derived from the scaling behavior of $`t_{med}`$; $`f_{Sh}`$ is the Shannon’s lower bound, eq. 4; Ratio is the efficiency of our code $`f_c/f_{Sh}`$; $`Z_R`$ indicates the gzip compression rate averaged over files of the sizes $`10^6`$ bits with the desired set of autocorrelations. We assume that the compression rate with $`L=10^6`$ achieves its asymptotic ratio, as was indeed confirmed in the compression of files with different $`L`$; $`1/R^{}`$ indicates the ideal (minimal) ratio between the transmitted message and the source signal after implementing the following two steps: compression of the file using gzip and then using an ideal optimal encoder/decoder, for a given BSC with $`f_c`$. A number greater than (less than) $`3`$ in this column indicates that the MN-DBP joint source-channel coding algorithm is more efficient (less efficient) in comparison to the channel separation method using the standard gzip compression. The last four columns of Table I are devoted for the comparison of our DBP algorithm to advanced compression methods. $`PPM_R`$ and $`AC_R`$ stand for the compression rate of files of the size $`10^6`$ bits with the desired autocorrelations using the Prediction by Partial Match and for the Arithmetic Coder, respectively. Similarly to the gzip case, $`1/R_{PPM}`$ and $`1/R_{AC}`$ stand for the optimal (minimal) rate required for the separation process (first a compression and then an ideal optimal encoder/decoder) assuming a BSC with $`f_c`$.

Table I indicates the following main results: (a) For $`q=4`$ (the upper part of Table I) a degradation in the performance is observed as the correlations are enhanced and as a result the entropy decreases. The degradation seems to be significant as the entropy is below $`0.3`$ (or for the test case $`R=1/3`$, $`f_c0.3`$). A similar degradation was observed as the entropy decreases for larger values of $`q`$. (b) The efficiency of our joint source-channel coding technique is superior to the alternative standard gzip compression in the source channel separation technique. For high entropy the gain of the MN-DBP is about $`510\%`$. This gain disappears as the entropy and the performance of the DBP algorithm are decreased. (c) In comparison to the standard gzip, the compression rate is improved by $`25\%`$ using the AC method. A further improvement of a few percents is achieved by the PPM compression. This later improvement seems to be significant in the event of low entropy. (d) Our DBP joint source-channel coding seems to be superior (by $`3\%`$)to the separation method based on the PPM compression for high entropy. However for ensemble of sequences characterized by low entropy this gain disappears. (e) With respect to the computational time of the source channel coding, our limited experience indicates that the DBP joint-source channel coding is faster than the AC separation method and the PPM separation method is substantially slower.

For a given set of autocorrelations where $`C_{k_0}`$ is the maximal one taken, the MN-DBP algorithm can be implemented with any field $`q2^{k_0}`$. If one wishes to optimize the complexity of the decoder it is clear that one has to work with the minimal allowed field, $`q=2^{k_0}`$. However, when the goal is to optimize the performance of the code and to maximize the threshold, the selection of the optimal field, $`q`$, is in question. In order to answer this question we present in Fig. 2 results for $`k_0=2`$ ($`C_1=C_2=0.86`$) and $`q=4,16,64`$. It is clear that the threshold, $`f_c`$, increases as a function of $`q`$ as was previously found for the case of unbiased signals. More precisely, the estimated thresholds for $`q=4,16,64`$ are $`0.293,0.3,0.309`$, respectively, and the corresponding Ratios ($`f_c/f_{Sh}`$) are $`0.913,0.934,0.962`$. where $`f_{Sh}=0.321`$. Note that the extrapolation of $`f_c`$ for large $`q`$ seems asymptotically to be consistent with $`f_c(q)0.3160.18/q`$.

For a given $`q`$, there are many sets of autocorrelations (or a finite fraction of $`\{C_k\}`$ in $`k_01`$ dimensions) obeying the same entropy. An interesting question is whether the performance of our DBP algorithm measured by the Ratio $`(f_c/f_{Sh})`$ is a function of the entropy only. Our numerical simulations indicate that the entropy is not the only parameter which controls the performance of the DBP algorithm. For the same entropy and $`q`$ the Ratio can fluctuate widely among different sets of correlations. For illustration, in Table II results for two sets of autocorrelations with the same entropy are summarized for each $`q=4,8,16`$ and $`32`$. It is clear that as the Ratio $`(f_c/f_{Sh})`$ is much degradated the gzip performance is superior (the second example with $`q=8`$ and $`32`$ in Table II where the Ratio is $`0.8`$ and $`0.72`$, respectively). The crucial question is to find the criterion to classify the performance of the DBP algorithm among all sets of autocorrelations obeying the same entropy. Our generic criterion is the decay of the correlation function over distances beyond two successive blocks. However, before the examination of this criterion, we would like to turn back to some aspects of statistical physics.

The entropy of sequences with given first $`k_0`$ correlations are determined via the effective Hamiltonian consisting of $`k_0`$ interactions, eqs. 2-3. As a result the entropy of these sequences is the same as the entropy of the effective Hamiltonian, $`H\{y_k\}`$, at the inverse temperature $`\beta =1`$, eq. 5. As for the usual scenario of the transfer matrix method, the leading order of quantities such as the free energy and the entropy are a function of the largest eigenvalue of the transfer matrix only. On the other hand the decay of the correlation function is a function of the whole spectrum of the $`2^{k_0}`$ eigenvalues. Asymptotically, the decay of the correlation function is determined from the ratio between the second largest eigenvalue and the largest eigenvalue, $`\lambda _2/\lambda _{max}`$. From the statistical mechanical point of view one may wonder, why the first $`k_0`$ correlations can be determined using the information of $`\lambda _{max}`$ only. The answer to this question is that once the transfer matrix is defined as a function of $`\{y_k\}`$, eqs. 3-5, all eigenvalues are determined as well as $`\lambda _{max}`$. There is no way to determine $`\lambda _{max}`$ independently of all other eigenvalues.

In Table II results of the DBP-MN algorithm for $`q=4,8,16,32`$ are presented. For each $`q`$, two different sets of autocorrelations characterized by the same entropy and threshold $`f_{Sh}`$ are examined. The practical method we used to generate different sets of autocorrelations with the same entropy was a simple Monte Carlo over the space of $`\{C_k\}`$. The additional column in Table II (in comparison with Table I) is the ratio between $`\lambda _2/\lambda _{max}`$, which characterizes the decay of the correlation function over large distances. Independent of $`q`$, it is clear that for a given entropy as $`\lambda _2/\lambda _{max}`$ increases/decreases the performance of the DBP algorithm measured by the Ratio $`f_c/f_{Sh}`$ is degradated/enhanced. The new criterion to classify the performance of the DBP algorithm among all sets of autocorrelations obeying the same entropy is the decay of the correlation function. This criterion is consistent with the tendency that as the first $`k_0`$ autocorrelations are increased/decreased a degradation/enhancement in the performance is observed (see Table I). The physical intuition is that as the correlation length increases, the relaxation time increases and flips on larger scales than nearest neighbor blocks are required.

Note that the decay of the correlation function in the intermediate region of a small number of blocks is a function of all the $`2^{k_0}`$ eigenvalues. Hence, in order to enhance the effect of the fast decay of the correlation function in the case of small $`\lambda _2/\lambda _{max}`$, we also try to enforce in our Monte Carlo search that all other $`2^{k_0}2`$ eigenvalues be less than $`A\lambda _{max}`$ with the minimal possible constant $`A`$. This additional constraint was easily fulfilled for $`q=4`$ with $`A=0.1`$, but for $`q=32`$ the minimal $`A`$ was around $`0.5`$.

Finally we raise the question of whether for a given entropy a degradation in the typical performance of the DBP algorithm is expected as $`q`$ increases. This is crucial since the superiority, if any, of the DBP joint source-channel coding method over advanced compression methods is in question. As explained above, our Monte Carlo simulations indicate that for a given entropy the suppression of the correlation function is more difficult as $`q`$ increases. This is a strong indication that as $`q`$ increases a degradation in the typical performance of the DBP decoder is expected, but its nature and significance have still to be examined in further research.

We thank Shlomo Shamai and David MacKay for valuable discussions.


text:
# General self-flattening surfaces

## Abstract

Recently Jeong and Kim \[Phys. Rev. E 66, 051605 (2002)\] investigated the scaling properties of equilibrium self-flattening surfaces subject to a restricted curvature constraint. In one dimension (1D), they found numerically that the stationary roughness exponent $`\alpha 0.561`$ and the window exponent $`\delta 0.423`$. We present an analytic argument for general self-flattening surfaces in $`D`$ dimensions, leading to $`\alpha =D\alpha _0/(D+\alpha _0)`$ and $`\delta =D/(D+\alpha _0)`$ where $`\alpha _0`$ is the roughness exponent for equilibrium surfaces without the self-flattening mechanism. In case of surfaces subject to a restricted curvature constraint, it is known exactly that $`\alpha _0=3/2`$ in 1D, which leads to $`\alpha =3/5`$ and $`\delta =2/5`$. Small discrepancies between our analytic values and their numerical values may be attributed to finite size effects.

Fluctuation properties of equilibrium surfaces have been studied extensively for the last few decades . Surface roughness is well documented and classified into a few universality classes. The Edwards-Wilkinson (EW) class is generic and robust for equilibrium surfaces with local surface tension . The EW surfaces can be described by the continuum Langevin-type equation

$$\frac{h(\stackrel{}{r},t)}{t}=\nu ^2h(\stackrel{}{r},t)+\eta (\stackrel{}{r},t),$$
(1)
where $`h(\stackrel{}{r},t)`$ is the height at site $`\stackrel{}{r}`$ and time $`t`$, $`\eta (\stackrel{}{r},t)`$ is an uncorrelated Gaussian noise, and $`\nu `$ represents the strength of local surface tension.

The surface fluctuation width $`W(L,t)`$, defined as the standard deviation of the surface height $`h(\stackrel{}{r},t)`$ starting from a flat surface of lateral size $`L`$, satisfies the dynamic scaling relation

$$W(L,t)=L^\alpha f\left(t/L^z\right),$$
(2)
where the scaling function $`f(x)\mathrm{const}.`$ for $`x1`$ and $`f(x)x^\beta `$ $`(\beta =\alpha /z)`$ for $`x1`$.

The EW universality class is characterized by the dynamic exponent $`z=2`$ and the roughness exponent $`\alpha =(2D)/2`$ ($`D`$: substrate dimension). In the absence of local surface tension $`(\nu =0)`$, higher-order local suppression terms like $`^{2m}h`$ $`(m=2,3,\mathrm{})`$ become relevant to determine the scaling properties of surface roughness. In this case, the scaling exponents become $`z=2m`$ and $`\alpha =(2mD)/2`$.

Recently Kim, Yoon, and Park  introduced a global-type suppression (self-flattening) mechanism which reduces growing (eroding) probability only at the globally highest (lowest) point on the surface. They found that this global-type suppression changes the scaling properties of the EW surfaces: $`z3/2`$ and $`\alpha =1/3`$ in 1D and $`z5/2`$ and $`\alpha =0(\mathrm{log})`$ in 2D. The 1D roughness exponent $`\alpha =1/3`$ characterizing the stationary surface fluctuations could be derived exactly by mapping the surface evolution model onto the static self-attracting random walk model .

More recently Jeong and Kim (JK)  investigated the effect of the self-flattening mechanism on the so-called restricted curvature (RC) model . The RC model is known to have no local surface tension term $`(\nu =0)`$ and its dominant suppression term is of the fourth order $`(^4h)`$. Accordingly, the ordinary RC model have the scaling exponents $`z=4`$ and $`\alpha =(4D)/2`$. Using Monte Carlo simulations, JK found for the 1D self-flattening RC model that

$$z=1.69(5),\alpha =0.561(5),\text{and}\beta =0.332(5).$$
(3)
Again, the self-flattening mechanism changes the scaling properties of the RC surfaces.

JK also studied the height-height correlation function $`G(r)`$, defined as the average of the square of height differences at two sites separated by a distance $`r`$. They found an extra length scale $`\xi `$ (window size) where the correlation function starts to saturate. The correlation function in the steady state satisfies the scaling relation

$$G(r)=L^{2\alpha }g\left(r/\xi \right)\text{with}\xi L^\delta ,$$
(4)
where the scaling function $`g(x)\mathrm{const}.`$ for $`x1`$ and $`f(x)x^{2\alpha ^{}\kappa }`$ $`(\alpha ^{}=\alpha /\delta )`$ for $`x1`$. This type of a crossover scaling, due to the existence of a smaller length scale $`(\delta <1)`$ compared to the system size, has been previously identified in the so-called even-visiting random walk (EVRW) models (see Eq. (29) in reference ). In fact, the EVRW model is intimately related to the EW-type surface model with the self-flattening mechanism . JK’s numerical estimates for the exponents are

$$\delta 0.423,\kappa 0.868,\text{and}\alpha ^{}=\alpha /\delta 1.33.$$
(5)
In this Comment, we present an analytic argument that predicts the values of the scaling exponents, $`\alpha `$ and $`\delta `$, associated with the stationary properties of surface fluctuations. We consider a general equilibrium surface growth model, of which the stationary roughness exponent is known exactly as $`\alpha _0`$. The partition function for equilibrium self-flattening surfaces of lateral size $`L`$ can be written as

$$Z_L(K)=\underset{𝒞}{}e^{KH(𝒞)},$$
(6)
where the summation is over all possible height configurations $`𝒞`$ subject to a given constraint, $`K`$ is a temperaturelike parameter, and $`H(𝒞)`$ is the height excursion width (the globally maximum height minus the globally minimum height) for a given configuration $`𝒞`$. Global suppression for self-flattening dynamics is simply Metropolis type evolution algorithm with this partition function to reach the equilibrium . For the EW surfaces, one can take the surface height configurations subject to the restricted solid-on-solid (RSOS) constraint, where the step heights are allowed to take finite values. In the case of the RC surfaces, the local curvature $`^2h`$ is restricted to be finite.

One can decompose the configurational space into sectors with a constant excursion width $`H`$. Then, the partition function can be rewritten as

$$Z_L(K)=_0^{\mathrm{}}𝑑H\omega _L(H)e^{KH},$$
(7)
where $`\omega _L(H)dH`$ is the number of configurations with the height excursion width between $`H`$ and $`H+dH`$. One can define $`\mathrm{\Omega }_L(H)`$ as the number of configurations with the height excursion width less than $`H`$ as

$$\mathrm{\Omega }_L(H)=_0^H𝑑H^{}\omega (H^{}),$$
(8)
and it is clear that $`\mathrm{\Omega }_L(\mathrm{})=Z_L(K=0)`$.

We estimate $`\mathrm{\Omega }_L(H)`$ as follows. Consider a flat surface between two walls separated by a distance $`H`$ and parallel to the surface. The surface starts to fluctuate, following its ordinary evolution dynamics. Whenever the surface hits either of the two walls, the number of possible configurations inside the walls will be reduced by a constant factor, compared to the no-wall case ($`H=\mathrm{}`$). This entropic reduction can be roughly translated as

$$\mathrm{\Omega }_L(H)\mathrm{\Omega }_L(\mathrm{})\mathrm{exp}[aN_c],$$
(9)
where $`N_c`$ is the typical number of contacts between the walls and the surface in the stationary state and $`a`$ is a positive constant of $`O(1)`$. Typically, there will be $`O(1)`$ contacts over a block of lateral size $`\mathrm{}`$, within which the stationary surface width $`W_0(\mathrm{})\mathrm{}^{\alpha _0}`$ is of the same order of magnitude as the wall distance $`H`$, i.e.,

$$N_c\left(L/\mathrm{}\right)^DL^D/H^{D/\alpha _0}.$$
(10)
With the above estimates, the partition function becomes for nonzero $`K`$

$`Z_L(K)`$ $`=`$ $`Z_L(0){\displaystyle _0^{\mathrm{}}}𝑑H\left({\displaystyle \frac{}{H}}e^{aL^D/H^{D/\alpha _0}}\right)e^{KH},`$ (11)
$`=`$ $`Z_L(0)K{\displaystyle _0^{\mathrm{}}}𝑑He^{KHaL^D/H^{D/\alpha _0}}.`$ (12)
This integral can be evaluated by the saddle point method in the limit of large $`L`$. The maximum contribution comes from

$$H^{}\left(\frac{L^D}{K}\right)^{\alpha _0/(D+\alpha _0)}.$$
(13)
The stationary surface fluctuation width should be proportional to this typical height excursion width, so we predict that

$$WK^{\alpha _0/(D+\alpha _0)}L^\alpha ,$$
(14)
with $`\alpha =D\alpha _0/(D+\alpha _0)`$.

The length scale $`\mathrm{}`$ arising naturally in our argument should be proportional to the window size $`\xi `$ in the correlation function. Inside of this length scale, the surface does not feel the global self-flattening suppression. As our length scale $`\mathrm{}`$ is explicitly related to the wall distance in Eq. (10), we can also predict that

$$\mathrm{}^{}\xi L^\delta ,$$
(15)
with $`\delta =D/(D+\alpha _0)`$.

In summary, for the general self-flattening equilibrium surfaces, we predict that

$$\alpha =\frac{D\alpha _0}{D+\alpha _0},\delta =\frac{D}{D+\alpha _0},\text{and}\alpha ^{}=\frac{\alpha }{\delta }=\alpha _0,$$
(16)
where $`\alpha _0`$ is the stationary roughness exponent for the equilibrium surfaces without the self-flattening mechanism.

For the EW surfaces, $`\alpha _0=1/2`$ in 1D, which leads to $`\alpha =1/3`$ and $`\delta =2/3`$. These results agree with those by the other analytic (healing time) arguments  and those by numerical simulations . For the RC surfaces, $`\alpha _0=3/2`$ in 1D, which leads to $`\alpha =3/5`$ and $`\delta =2/5`$. The short range behavior governing the value of $`\kappa `$ should be identical to the ordinary RC model, so we also expect that $`\kappa =1`$. These results have small discrepancies from the JK’s numerical results in Eqs. (3) and (5). We believe that this may be due to rather small sizes used in their numerical simulations. The RC surfaces have $`\alpha _0=(4D)/2`$ in $`D`$ dimensions, so $`\alpha =D(4D)/(4+D)`$ and $`\delta =2D/(4+D)`$. It may be interesting to check these predictions for the RC self-flattening surfaces in $`D=2`$ and 3.

At the upper critical dimensions ($`D=2`$ for the EW and $`D=4`$ for the RC surfaces), the surface roughness becomes logarithmic ($`\alpha _0=0`$) and the self-flattening mechanism induces only corrections to scaling in the surface fluctuations. The dominant correction seems to be independent of system size $`L`$ , but its functional dependence on $`K`$ is not fully explored. In higher dimensions, the EW $`(D>2)`$and the RC $`(D>4)`$ surfaces are asymptotically flat ($`\alpha _0<0`$). The self-flattening mechanism induces power-law-type corrections to scaling as given in Eq. (16).

We could not present any analytic explanation for the dynamic exponents $`z`$ and $`\beta `$. This may be due to the lack of a continuum-type equation to govern the self-flattening mechanism. It would be very interesting to find such an equation which should contain a global coupling term in space.

We thank J. M. Kim and H.-C. Jeong for bringing us their work  in our attention. This work was supported by Grant No. 2000-2-11200-002-3 from the Basic Research Program of KOSEF.


text:
# Moving a perturbation across quantum dots: tuning of Fano resonances

## I acknowledgments

M. Mendoza would like to acknowledge the Brazilian agency CAPES for financial support, while P.A.S. is grateful to the continuous support provided by FAPESP.


text:
# A connection between the morphology of the X-ray emission and rotation for solar-type stars in open clusters1footnote 11footnote 1This is paper 14 in the WIYN Open Cluster Study (WOCS).

## 1 Introduction

The soft ($`<2.5`$ keV) X-ray emission from late-type cluster stars is currently believed to delineate three kinds of dependence on stellar rotation or more precisely, Rossby number. This was originally noted in Patten & Simon (1996), was reviewed in Randich (2000), and has been reiterated recently by Feigelson et al. (2003). Other studies have noted similar results for cluster and field stars (eg. James et al. 2000). In the regime of slow rotation, a clear rotation-activity relationship is observed, with X-ray emission declining steadily with increasing rotation period. For faster rotators, there is evidence of X-ray “saturation,” a phenomenon describing a relative constancy of the X-ray emission at an elevated level regardless of increasing rotation speed. Finally, for the very fastest rotators, there appears to be a moderate decline in the X-ray emission, and this phenomenon is sometimes called “super-saturation.” The physical basis for these effects is presently unclear.

It has recently been proposed that the morphology of the rotational observations of cluster and field stars may be understood in terms of a magnetic paradigm (Barnes 2003). According to this paradigm, the fastest rotators in open clusters lie on a convective sequence, where a convective magnetic field dominates, and the outer convection zone of the stars shears past the interior radiative zone. It is proposed that this shear generates an interface dynamo which attempts to couple the two zones, driving the star through a gap (in rotation) region, towards the interface sequence, one occupied by stars with well-developed dynamos, and including the present-day Sun.

If this paradigm is valid, it should also be reflected in the X-ray properties of stars, since these are believed to originate in stellar magnetism. Here, I attempt to show that in the case of specific open clusters where high-quality rotation periods and X-ray data are both available, the X-ray emission does, in fact, conform to the new paradigm.

## 2 Observations

A meaningful test of the paradigm requires clusters that have not only been well-studied in X-rays but where the cluster stars can be placed on the underlying rotational sequences with absolute unambiguity. This latter requirement implies the use of rotation periods.

Four such groups of stars are included here - IC 2391:  Patten & Simon (1996); IC 2602:  Barnes et al. (1999); Alpha Per:  Prosser & Grankin (1997); and a group of field M dwarfs from James et al. (2000). Using the $`BV`$ color and rotation period, it is possible to classify each star in this set in terms of the paradigm advanced in Barnes (2003) as lying on the convective sequence (C), interface sequence (I), or in the rotational gap (g) between these sequences. This classification is displayed graphically in the left-hand panels of Figure 1. It is based on the rotation period of the star in question, taking account of its (dereddened) color and age, except in the case of the M dwarfs whose heterogeneous ages do not permit an unambiguous identification with the relevant sequence. An ambiguity discussed thoroughly in Barnes (2003) with respect to very late-type stars suggests a more cautious labelling of the three slowest rotators in this set merely by the label “S.”

Having classified each star as C, g, or I/S, we consider the X-ray properties. For each of these stars, the soft X-ray emission values are available from the following studies - IC 2391:  Patten & Simon (1996); IC 2602:  Randich et al. (1995); Alpha Per:  Randich et al. (1996) and Prosser et al. (1996); and for field M stars: James et al. (2000). In the right-hand panels of Figure 1, we correspondingly plot the X-ray emission from each of these stars222James et al. (2000) suggest that some caution is warranted in interpreting data where the bolometric corrections for late-type stars are calculated using $`BV`$ colors. against the Rossby Number, $`N_R=P/\tau _c`$, where $`P`$ is the rotation period of the star and $`\tau _c`$ is the convective turnover timescale, estimated using the relationship from Kim and Demarque (1996). In these plots, as before, we mark each star by the symbol C, g, or I/S, based on the classification described earlier.

One observes that the stars marked C lie in the super-saturated region of X-ray activity, those marked g in the saturated region, and those marked I in the region of validity of the rotation-activity relationship.

This picture is further reinforced by Figure 1 in Randich (2000), which shows that all the solar-type field stars, expected to lie on the Interface sequence (I) in the new paradigm, obey the standard rotation-activity relationship.

## 3 Interpretation

The foregoing observations encourage us to think that the paradigm suggested to account for the rotational evolution of stars might also be useful in understanding the X-ray emission from stars.

Purely observationally, we know that the X-ray emission $`X=X(N_R)=X(P/\tau _c)`$. Now we know from Barnes (2003) that $`P`$ is not a smoothly varying function of stellar properties but rather that

$$P=\{\begin{array}{cc}P_C\hfill & \text{for stars on the convective sequence;}\hfill \\ P_g\hfill & \text{for stars in the gap;}\hfill \\ P_I\hfill & \text{for stars on the interface sequence.}\hfill \end{array}$$
(1)
where $`P_C`$ and $`P_I`$ have separate functional dependences on stellar age and color, and $`P_g`$ merely represents stars with rotation periods in the gap between the Convective and Interface sequences. Very briefly, $`P_I\sqrt{t/I_s}`$ but $`P_Ce^{t/(KI_c)}`$ where $`t`$ is the age, $`I_s`$ and $`I_c`$ are, respectively, the moments of inertia of the entire star or only the convection zone, and $`K`$ is a constant. The origin of these sequences in terms of stellar properties has been explained at length in Barnes (2003).

Thus, the fact that the X-ray emission has the particular three-segmented morphology that has been noted by various studies can be traced to the fact that, for the three sections, we have

$$X=\{\begin{array}{cc}X(P_C/\tau _c)\hfill & \text{for stars on the convective sequence;}\hfill \\ X(P_g/\tau _c)\hfill & \text{for stars in the gap;}\hfill \\ X(P_I/\tau _c)\hfill & \text{for stars on the interface sequence.}\hfill \end{array}$$
(2)
Although this scheme suggests the immediate reason for the differing dependences on Rossby number, the underlying reason must have to do with the differing magnetic fields on these sequences that are manifested in differing X-ray emission.

If the X-ray emission scales simply with the total closed topology magnetic field, then we have a scheme here which roughly parallels the suggestion in Barnes (2003) that $`CgI`$ represents an evolutionary progression for solar-type stars. In that case, the right-hand panels of Figure 1 imply, in order, the creation/strengthening of a magnetic field on the C sequence (super-saturation in X-ray parlance), the attainment of a maximal field in the gap g (saturation), and the steady decay of this field for stars on the I sequence (rotation-activity paradigm).

The paradigm suggested for the rotational evolution of stars is strikingly parallel. It was suggested that stars on the convective sequence are initially equipped with only a convective or small-scale magnetic field which cannot couple to the radiative interior. Thus, only the external convection zone is braked. This results in a large shear at the convective/radiative interface, which creates an interface-type dynamo and magnetic field. The field itself then couples the convective and radiative zones of the star, and drives the star through the gap to the interface sequence, where the entire star is braked subsequently. In this picture, the X-ray observations would therefore be interpreted to indicate the creation and strengthening of the interface dynamo among the C sequence stars, the attainment of a maximal field for the g stars, and the decay of the interface field among the I sequence stars.

## 4 Chromospheric Ca II H & K emission

The Chromospheric Ca II H & K emission from solar-type stars in the field has been known to decline steadily with age (eg. Wilson 1963; Skumanich 1972). In fact, this decline is well-behaved enough that it is currently used to derive stellar ages for field stars (Soderblom, Duncan & Johnson 1991; Donahue 1998). These stars might reasonably be expected to be similar to those that obey the rotation-Xray activity paradigm, and to be located on the interface sequence since they are mostly older than several hundred megayears.

In young clusters, however, there appear to be complex patterns of Ca II H & K emission that cannot easily be mapped onto the standard age-rotation-activity paradigm (C.M. Anderson, 2002, personal communication). If the new rotational paradigm is correct, it ought to be able to explain these variations in a manner similar to the X-ray variations. This task is currently underway.

## 5 Conclusion

The X-ray observations of cluster and field stars have a three-segmented morphology largely because there are three kinds of rotating stars - those with only a convective magnetic field, here identified with the super-saturation segment, those with a well-developed interface magnetic field, here identified with those stars obeying the rotation-activity relationship, and those in the process of conversion from the former to the latter state, identified with the saturated activity segment.

Because rotating stars evolve systematically in the direction: $`CgI`$, the X-ray morphology of cluster and field stars might merely reflect the build-up of large scale magnetic fields on the convective sequence during this high-shear era, the attainment of a maximal field (saturation) for the gap stars, and the gradual decline of X-rays with magnetic field and rotation on the interface sequence.

SAB would like to acknowledge support from the NSF through AST-9986962.


text:
# Logic and thermodynamics: the heat-engine axiomatics of the second law

## Abstract

We challenge the statement that the principle of Thomson and the principle of Clausius are equivalent. A logical mistake in the supposed textbook proof of their equivalency is indicated. On this account we refine the heat-engine axiomatics. We consider the energy exchange in the configuration comprised of two heat reservoirs and one mechanical device and show explicitly the domains banned by the laws of thermodynamics.

In the preface to the textbook Kubo  on thermodynamics Ryogo Kubo passed the following observation. ”As in contrast to the atomic theory, thermodynamics does not find a support in our intuition. This is one of the reasons why students consider thermodynamics difficult for mastering and can not apply it to concrete problems.” As a confirmation of his words we see Dmitriyev  a confusion in the very foundations of the theory. There is an error in the axiomatics of phenomenological thermodynamics or, may be, imperfection, that moves over from one textbook to another for a century (see e.g. Ref. Huang, ; Munster, ; Goldstein, ). This is the well-known proof by contradiction of that the principle of Thomson is equivalent to the principle of Clausius and vice versa. This theorem is mostly decorative and so does not prevent us from the proper conclusions. However, it may be just the reason of the difficulties in comprehending the whole structure of thermodynamics. We use it as a guide to make the formulation of the second law more coherent and unequivocal. The role of formal logic in constructing deductive theories is thus emphasized.

We consider the exchange of the energy in the ternary configuration comprised of the two thermostats $`h`$ and $`c`$ and one mechanical system $`m`$. The thermomechanical device that implements this exchange in a cycle is known as the heat engine. The first law restricts the values $`\mathrm{\Delta }E`$ of the increments of the energy by the plane

$$\mathrm{\Delta }E_h+\mathrm{\Delta }E_m+\mathrm{\Delta }E_c=0.$$
(1)
The second law puts further restraints on the vector $`(\mathrm{\Delta }E_h,\mathrm{\Delta }E_m,\mathrm{\Delta }E_c)`$. The principle of Thomson forbids the region in the plane (1) where

$$\mathrm{\Delta }E_h0,\mathrm{\Delta }E_m>0,\mathrm{\Delta }E_c0.$$
(2)
The principle of Clausius forbids in the plane (1) the ray

$$(\mathrm{\Delta }E_h,0,\mathrm{\Delta }E_c):\mathrm{\Delta }E_h>0$$
(3)
and thus ranks the bodies into cold $`c`$ and hot $`h`$. The principles (2) and (3) do not specify the whole region banned by the second law. But they play a crucial role in constructing this region. As we see, these principles belong to different domains. So, they can not be equivalent.

Using the both principles we may prove Huang  the following lemma:

$$\text{if}\mathrm{\Delta }E_m>0\text{then}\mathrm{\Delta }E_c>0.$$
(4)
Indeed, let $`\mathrm{\Delta }E_c=0`$. Then we have from Eq. (1) $`\mathrm{\Delta }E_h<0`$, that contradicts the Thomson’s principle (2). Let $`\mathrm{\Delta }E_c<0`$. The Thomson’s principle (2) does not forbid the process

$$\mathrm{\Delta }E_h>0,\mathrm{\Delta }E_m>0,\mathrm{\Delta }E_c<0.$$
(5)
It also allows a process

$$\mathrm{\Delta }^{}E_h>0,\mathrm{\Delta }^{}E_m<0,\mathrm{\Delta }^{}E_c=0.$$
(6)
Summing up (5) and (6) for $`\mathrm{\Delta }^{}E_m=\mathrm{\Delta }E_m`$ we get

$$\mathrm{\Delta }E_h+\mathrm{\Delta }^{}E_h>0,\mathrm{\Delta }E_m+\mathrm{\Delta }^{}E_m=0,\mathrm{\Delta }E_c+\mathrm{\Delta }^{}E_c<0.$$
(7)
The process (7) contradicts the Clausius’ principle (3). That proves lemma (4).

To find the region allowed by thermodynamics we accept the postulate of the existence of reversible processes ($`r`$). It states that for each $`\mathrm{\Delta }E_h^{(r)}`$ there exists $`\mathrm{\Delta }E_c^{(r)}`$ such that we can find $`\mathrm{\Delta }^{}E_c^{(r)}`$ and $`\mathrm{\Delta }^{}E_h^{(r)}`$ complying with

$`\mathrm{\Delta }E_h^{(r)}+\mathrm{\Delta }^{}E_h^{(r)}`$ $`=`$ $`0`$ (8)
$`\mathrm{\Delta }E_c^{(r)}+\mathrm{\Delta }^{}E_c^{(r)}`$ $`=`$ $`0.`$ (9)
The reverse process is unique that is easily proved using the Thomson’s principle (2).

Proceeding from (4) and (8), (9) we may prove the Carno’s theorem. It says that the efficiency of a reversible heat engine is greater or equal to the efficiency of any other heat engine. Below we give a more accurate proof of the Carno’s theorem than it is usually done Kubo .

Consider the energy exchange in the ternary configuration according to Eq. (1). We want to compare (1) with the reversible process

$$\mathrm{\Delta }E_h^{(r)}+\mathrm{\Delta }E_m^{(r)}+\mathrm{\Delta }E_c^{(r)}=0$$
(10)
under the condition

$$\mathrm{\Delta }E_h=\mathrm{\Delta }E_h^{(r)}.$$
(11)
To this end, we will carry out firstly (1) and then the process

$$\mathrm{\Delta }^{}E_h^{(r)}+\mathrm{\Delta }^{}E_m^{(r)}+\mathrm{\Delta }^{}E_c^{(r)}=0$$
(12)
that is reverse to (10). Summing up Eq. (1) and Eq. (12) and using in it (11) and (8) we get

$$\mathrm{\Delta }E_m+\mathrm{\Delta }^{}E_m^{(r)}+\mathrm{\Delta }E_c+\mathrm{\Delta }^{}E_c^{(r)}=0.$$
(13)
By virtue of the Thomson’s principle (2) we have from Eq. (13)

$$\mathrm{\Delta }E_c+\mathrm{\Delta }^{}E_c^{(r)}0.$$
(14)
Excluding $`\mathrm{\Delta }^{}E_c^{(r)}`$ from Eq. (14) and (9) gives

$$\mathrm{\Delta }E_c\mathrm{\Delta }E_c^{(r)}0.$$
(15)
Let $`\mathrm{\Delta }E_m^{(r)}>0`$. Then by the lemma (4) $`\mathrm{\Delta }E_c^{(r)}>0`$. Inequality (15) says that with the condition (11) the work made by the reversible heat engine is greater than the work done by the irreversible heat engine. Dividing (15) by $`\mathrm{\Delta }E_c^{(r)}`$ and using in it tautologically (11) we get

$$\frac{\mathrm{\Delta }E_h}{\mathrm{\Delta }E_h^{(r)}}+\frac{\mathrm{\Delta }E_c}{\mathrm{\Delta }E_c^{(r)}}0.$$
(16)
Let $`\mathrm{\Delta }E_m^{(r)}<0`$. Then by the lemma (4) and (9) we have for reversible processes $`\mathrm{\Delta }E_c^{(r)}<0`$. In this case inequality (16) should be replaced by

$$\frac{\mathrm{\Delta }E_h}{\mathrm{\Delta }E_h^{(r)}}\frac{\mathrm{\Delta }E_c}{\mathrm{\Delta }E_c^{(r)}}0.$$
(17)
At last, we may unite (16) and (17) into the inequality

$$\frac{\mathrm{\Delta }E_h}{|\mathrm{\Delta }E_h^{(r)}|}+\frac{\mathrm{\Delta }E_c}{|\mathrm{\Delta }E_c^{(r)}|}0.$$
(18)
Combine $`n`$ irreversible and $`n^{(r)}`$ reversible heat engines such that

$$\frac{\mathrm{\Delta }E_h}{\mathrm{\Delta }E_h^{(r)}}=\frac{n^{(r)}}{n}.$$
(19)
Then we may deduce Huang  (18) for a general case that does not obey the restriction (11). As you see, the proof of the Carno’s theorem makes use of the Thomson’s principle (2). It can also be proved using the Clausius’ principle (3).

Let us take instead of (1) a reversible process $`(r_0)`$. Then, reversing the processes, we may turn Kubo  (18) into equality. This equality can be written as

$$\frac{|\mathrm{\Delta }E_h^{(0)}|}{|\mathrm{\Delta }E_h^{(r)}|}=\frac{|\mathrm{\Delta }E_c^{(0)}|}{|\mathrm{\Delta }E_c^{(r)}|}.$$
(20)
Now we can define Kubo ; Huang  an absolute temperature scale by the relations

$$|\mathrm{\Delta }E_h^{(0)}|=\kappa \theta _h,|\mathrm{\Delta }E_c^{(0)}|=\kappa \theta _c$$
(21)
where $`\kappa >0`$ is a constant. Substituting (20) and (21) into (18) gives

$$\frac{\mathrm{\Delta }E_h}{\theta _h}+\frac{\mathrm{\Delta }E_c}{\theta _c}0.$$
(22)
Taken together with Eq. (1), inequality (22) specifies the domain of the processes allowed by thermodynamics. This is a half-plane in the energy space $`(\mathrm{\Delta }E_h,\mathrm{\Delta }E_m,\mathrm{\Delta }E_c)`$.

As we see from the above, the principle of Thomson and the principle of Clausius are not equivalent. However, there is a popular textbook proof Huang ; Munster ; Goldstein  of their equivalence. In order to show explicitly the error in the supposed proof we must give both the verbal and logical formulation of these principles. The Thomson’s principle states: One can not take heat from a body and totally convert it to work. The Clausius’ principles states: ”Heat cannot, of itself, pass from a colder to a hotter body.” The equivalence is proved by contradiction (i.e. from the contrary). The two theorems are needed.

A: If the Thomson’s principle is wrong then the Clausius’ principle is wrong as well.

B: If the Clausius’ principle is wrong then the Thomson’s principle is wrong as well.

It is sufficient to analyze the supposed proof of the first theorem. It runs as follows.

1. Let the Thomson’s principle be wrong.

2. Then we may take a heat from the cold body and convert it wholly to a work.

3. Next we may convert this work into the internal energy of the hot body.

4. The pure result of these operations is that heat was transferred from the cold body to the hot body.

That contradicts the Clausius’ principle and hence the statement A is valid.

Here the error is in the step 3. To explicate it we formulate the Thomson’s principle in logical terms. Let us define on the line $`\mathrm{\Delta }E+\mathrm{\Delta }E_m=0`$ the Boolean function that describes the energy exchange between a termostat and a mechanical system:

$`I(\mathrm{\Delta }E_m)`$ $`=`$ $`true\text{for}\mathrm{\Delta }E_m0`$ (23)
$`I(\mathrm{\Delta }E_m)`$ $`=`$ $`false\text{for}\mathrm{\Delta }E_m>0.`$ (24)
This function just expresses the Thomson’s principle. Negating it we get

$`{}_{}{}^{\neg }I(\mathrm{\Delta }E_m)`$ $`=`$ $`false\text{for}\mathrm{\Delta }E_m0`$ (25)
$`{}_{}{}^{\neg }I(\mathrm{\Delta }E_m)`$ $`=`$ $`true\text{for}\mathrm{\Delta }E_m>0.`$ (26)
This allows $`\mathrm{\Delta }E_m>0`$ and forbids $`\mathrm{\Delta }E_m0`$. We made use (26) in step 2 and violated (25) in step 3. The origin of this mistake is in that the verbal formulation of the Thomson’s principle says explicitly (24) and does not speak out (23). However, the latter is implicit as in an ellipsis or enthymeme. The literal understanding of the words gave rise to the error in the logical inference. This illustrates the significance of formal logic in constructing deductive theories.


text:
# A clock transition for a future optical frequency standard with trapped atoms

## Abstract

We report the first direct excitation of the strongly forbidden $`5s^2{}_{}{}^{1}S_{0}^{}5s5p{}_{}{}^{3}P_{0}^{}`$ transition in 87Sr. Its frequency is 429 228 004 235 (20) kHz. A resonant laser creates a small leak in a magneto-optical trap (MOT): atoms build up to the metastable $`{}_{}{}^{3}P_{0}^{}`$ state and escape the trapping process, leading to a detectable decrease in the MOT fluorescence. This line has a natural width of $`10^3`$ Hz and can be used for a new generation of optical frequency standards using atoms trapped in a light shift free dipole trap.

In addition to being at the heart of the international system of units, extreme frequency metrology finds a wide range of applications throughout physics. As an example, tests of the stability of fundamental constants are performed by comparing the frequencies of different atomic transitions Marion03 ; Bize03 . These laboratory tests are already competitive with those performed at the cosmological scale Damour96 ; Webb01  and improve with the accuracy of the atomic transitions frequency measurements. Microwave frequency standards such as atomic fountains now have relative accuracies better than $`10^{15}`$ with potential improvements down to $`10^{16}`$ Bize02 . Going further, however, seems extremely difficult due to the limited quality factor of the atomic resonance in these devices ($`10^{10}`$). Optical frequency standards exhibit much higher line-$`Q`$, up to $`10^{14}`$ as recently demonstrated in Ref. Rafac00 . These devices reach fractional frequency instabilities below $`10^{14}`$ over 1 s Diddams01 , almost one order of magnitude better than fountains Bize02 . Their current accuracy lies near 1 part in $`10^{14}`$ with anticipated large room for improvement Niering00 ; Bize03 ; Riehle02 ; Curtis02 ; StengerOL01 .

Two different approaches are commonly used to develop optical frequency standards StAndrews . The first one is based on the spectroscopy of a single trapped ion, the second on the spectroscopy of a large ensemble of free falling neutral atoms. It is a commonly shared opinion that the ion approach may lead to a better ultimate frequency accuracy due to the ”perfect” control of the ion motion, while the atomic approach should lead to a better frequency stability thanks to the numerous quantum references contributing to the signal. Recently H. Katori proposed a scheme which would combine the advantages of both approaches Katori02 . The idea is to trap neutral atoms in the Lamb-Dicke regime in an optical lattice operating at a wavelength where the light-shift of the clock transition vanishes.

H. Katori proposed to use 87Sr probed on the $`5s^2{}_{}{}^{1}S_{0}^{}5s5p{}_{}{}^{3}P_{0}^{}`$ line at 698 nm, which indeed seems an ideal system for the realization of this scheme. This $`J=0J=0`$ resonance is only slightly allowed by hyperfine coupling and its metrological properties are exquisite. A natural linewidth of $`1`$ mHz can be derived from the hyperfine data found in Ref. Kluge74 . It has a high insensitivity to external electromagnetic fields. The light-shift cancellation is expected to occur near 800 nm, a wavelength far from any atomic resonance and for which powerful and practical laser sources are readily available. Apart from a residual effect due to hyperfine structure, the light shift cancellation is independent on polarization. Finally, higher order effects due to the trapping field are expected to be extremely small on this particular line Pal02 . In this Letter we report the first direct observation and frequency measurement of this transition.

The experiment is performed on a sample of cold atoms collected in a magneto-optical trap (MOT). An atomic beam is decelerated in a Zeeman slower and captured at the crossing-point of three retro-reflected beams tuned 40 MHz to the red of the $`{}_{}{}^{1}S_{0}^{}^1P_1(F=9/2)`$ transition (Fig. 1). The setup is described in details in Ref. Courtillot03 . For this particular experiment $`3\times 10^6`$ 87Sr atoms are trapped at steady state at a temperature of 2 mK. The atomic cloud has a $`1/e^2`$ diameter of $`2`$ mm. To probe the 698 nm line, a 14 mW laser beam passes four times through the atomic cloud in a standing wave configuration. The probe beam waist radius is 1.3 mm. At first sight, a direct detection of the very forbidden line with our system seems a desperate task: the resonance is expected to be Doppler broadened to 1.5 MHz (FWHM) more than three orders of magnitude higher than the Rabi frequency (0.8 kHz). Thus only $`10^3`$ atoms are expected to be resonant with the laser at a time. In addition, an efficient fluorescence detection of atoms in the $`{}_{}{}^{3}P_{0}^{}`$ state is made difficult by the absence of cycling transition from this state. Finally the energy of the $`{}_{}{}^{3}P_{0}^{}`$ state is known from databases to within a few hundred MHz only.

To overcome these problems we have first determined the resonance frequency with an accuracy of 110 kHz by measuring the frequency of the $`{}_{}{}^{1}S_{0}^{}^3P_1`$ line at 689 nm and of the frequency difference between the $`{}_{}{}^{3}P_{1}^{}^3S_1`$ and $`{}_{}{}^{3}P_{0}^{}^3S_1`$ transitions at 688 nm and 679 nm respectively (see Fig 1) Note . Second, to amplify the fraction of cold atoms transferred to the $`{}_{}{}^{3}P_{0}^{}`$ state by direct excitation, we have taken advantage of the fact that the lifetime of the MOT is 40 times larger than the Rabi oscillation period. A laser tuned to resonance induces a leak in the MOT leading to a detectable decrease of the number of trapped atoms of 1%.

The experimental set up used for all frequency measurements is described in Fig. 2. An extended cavity laser diode (ECLD1) is locked to a high finesse cavity using the Pound-Drever-Hall method Drever83  with performances reported in Ref. Quessada03 . Its frequency is continuously measured vs a hydrogen maser with a scheme based on a self referenced femtosecond Ti:Sapph laser Holzwarth00 ; Jones00 . The relative resolution of this measurement is typically $`3\times 10^{13}`$ for a one second averaging time. A second laser (ECLD2) is offset-phase locked to ECLD1. The beat-note between both lasers is mixed with the output of a radio frequency synthesizer to generate the offset phase lock error signal. The bandwidth of the servo control is 2 MHz. With this scheme the light of ECLD2, which is sent to the atoms, can be tuned to any frequency between two modes of the cavity of free spectral range 1.5 GHz by actuating the RF synthesizer. Two sets of lasers are used. One can be tuned from 675 to 685 nm, the other one from 685 to 698 nm.

The frequency of the $`{}_{}{}^{1}S_{0}^{}^3P_1(F=9/2)`$ transition is measured with an atomic beam independent from the cold atom setup. The obtained frequency is 434 829 342 950 (100) kHz, with an uncertainty mainly dominated by statistical noise and imperfect knowledge of the magnetic field.

The determination of the fine structure splitting between $`{}_{}{}^{3}P_{0}^{}`$ and $`{}_{}{}^{3}P_{1}^{}`$ is performed with the cold atoms. While trapped and cycling on the $`{}_{}{}^{1}S_{0}^{}^1P_1`$ transition, atoms eventually emit a spontaneous photon which brings them to the $`{}_{}{}^{1}D_{2}^{}`$ state (Fig. 1). This state has two main decay channels: to the $`{}_{}{}^{3}P_{1}^{}`$ and $`{}_{}{}^{3}P_{2}^{}`$ states. Atoms in the $`{}_{}{}^{3}P_{1}^{}`$ state decay back to the ground state and are kept in the trap while atoms in the metastable $`{}_{}{}^{3}P_{2}^{}`$ state are lost. This process limits the lifetime of the MOT to some 50 ms. We make the fine structure measurement by modifying this escape process. If a laser resonant to one of the hyperfine components of the $`{}_{}{}^{3}P_{1}^{}^3S_1`$ transition is added to the trap, atoms in the corresponding $`{}_{}{}^{3}P_{1}^{}`$ state are pumped to the $`{}_{}{}^{3}P_{2}^{}`$ and $`{}_{}{}^{3}P_{0}^{}`$ metastable states. They escape the trap instead of decaying back to the ground state. This decreases the trapped atom number. Fig. 3(a) shows the fluorescence of the MOT as a function of the 688 nm laser detuning from the $`{}_{}{}^{3}P_{1}^{},F=9/2^3S_1,F=11/2`$ resonance.

A direct detection of the $`{}_{}{}^{3}P_{0}^{}^3S_1`$ transition by the same technique is not possible because $`{}_{}{}^{3}P_{0}^{}`$ is not populated in the MOT. Instead, a first signal is obtained by detecting the light shift of the $`{}_{}{}^{3}P_{1}^{}^3S_1`$ transition induced by a 679 nm laser close to the $`{}_{}{}^{3}P_{0}^{}^3S_1`$ resonance. With an intensity of 1.8 mW/mm2, a shift of 100 kHz is observed for a detuning of 100 MHz. For a 679 nm detuning smaller than the width of the 688 nm resonance, a dip appears in the resonance profile due to coherent population trapping (CPT): when the frequency difference between both lasers matches the atomic fine structure, there exists a coherent superposition of $`{}_{}{}^{3}P_{1}^{}`$ and $`{}_{}{}^{3}P_{0}^{}`$ states which is not coupled to $`{}_{}{}^{3}S_{1}^{}`$ Alzetta76 . Atoms in this dark state are not pumped to $`{}_{}{}^{3}P_{2}^{}`$ but decay back to the ground state in a few 10 $`\mu `$s due to the $`{}_{}{}^{3}P_{1}^{}`$ instability. They are kept in the MOT. When the CPT dip is centered on the 688 nm resonance, the 679 nm laser is tuned to resonance. The observed signal in this configuration is shown in Fig. 3(b). The fine structure measurement is performed light shift free with both lasers locked on resonance. We measured a frequency of 5 601 338 650 (50) kHz with an uncertainty mainly due to the magnetic field gradient of the MOT. According to these measurements, the $`{}_{}{}^{1}S_{0}^{}^3P_0`$ transition is expected to have a frequency of 429 228 004 300 (110) kHz.

For the direct observation, we induce in the MOT a leak to the $`{}_{}{}^{3}P_{0}^{}`$ state with a laser tuned to resonance. The lifetime of the trap is two orders of magnitude longer than the duration of a $`\pi `$ pulse on the forbidden transition with the laser parameters given above. This leads to a build-up by the same factor of the fraction of atoms escaping the MOT if the transfer rate to $`{}_{}{}^{3}P_{0}^{}`$ is constant and if the atoms, once in the $`{}_{}{}^{3}P_{0}^{}`$ state, actually escape the trapping process. One way to fulfil the first condition could consist in using the MOT to rethermalize atoms to fill the dip in the velocity distribution created by excitation to $`{}_{}{}^{3}P_{0}^{}`$. For the escape condition, one could think of optical pumping to $`{}_{}{}^{3}P_{2}^{}`$ with a laser tuned to $`{}_{}{}^{3}P_{0}^{}^3S_1`$. In our experiment however, the Doppler effect induced by gravity is sufficient to fulfil both conditions. It amounts to several times the Rabi frequency per millisecond with the $`45^{}`$ angle formed by the 698 nm probe beam and vertical. The dip in the velocity distribution induced by the excitation laser is then permanently refilled. On the other hand, atoms in the $`{}_{}{}^{3}P_{0}^{}`$ state are rapidly detuned from the excitation laser and do not emit a stimulated photon back to the ground state. The trapped atom number should then decrease by several %.

In practice, the experiment is operated sequentially: by means of acousto-optic modulators we alternate a capture and cooling phase with the blue lasers and a probe phase with the 698 nm laser. The main motivation for this pulsed operation is to avoid any light shift of the forbidden transition. Another reason is related to the escape process of the $`{}_{}{}^{3}P_{0}^{}`$ atoms. It would be less efficient with the MOT beams on due to the instability of the $`{}_{}{}^{1}S_{0}^{}`$ state confered by the coupling to $`{}_{}{}^{1}P_{1}^{}`$: the reduced lifetime of the atomic coherences due to spontaneous emission of blue photons would lead to an effective broadening of the 698 nm excitation much larger than the Doppler effect due to gravity. The optimization of the time sequence results from a trade-off between the capture efficiency of the MOT, the ballistic expansion of the atomic cloud during probe phases, and the efficiency of excitation to $`{}_{}{}^{3}P_{0}^{}`$. With phases of capture of 3 ms duration and probe phases of 1 ms duration, we have $`1\times 10^6`$ atoms at steady state in the MOT and the contrast of the resonance is 1%. In figure 4 is shown the fluorescence of the trapped atoms versus the 698 nm laser detuning from resonance. A narrow sub-Doppler structure is expected at the center of the resonance due to the standing wave configuration. With the present signal to noise ratio, we are not able to detect it.

We have locked the laser to the $`{}_{}{}^{1}S_{0}^{}^3P_0`$ resonance. Its frequency is 429 228 004 235 kHz with a standard deviation of 20 kHz for an averaging time of 2 hours. At that level, the systematic effects are negligible. The first order Doppler effect vanishes due to the standing wave configuration. The Zeeman effect induced by the MOT field gradient is less than 1 kHz, since the Landé factors of the $`{}_{}{}^{1}S_{0}^{}`$ and $`{}_{}{}^{3}P_{0}^{}`$ levels are both equal to $`1.3\times 10^4`$. The light shift which would result from imperfect extinction of the blue light is also estimated below 1 kHz. Tab. 1 summarizes the indirect and direct measurements and confirms the good agreement between both methods.

The technique used for the direct detection of this line can be extended to other atomic species. As an example, Yb has the same level scheme as Sr with several isotopes having non-zero nuclear spin. With Yb atoms trapped on the $`{}_{}{}^{1}S_{0}^{}^3P_1`$ transition of natural width 200 kHz, the lifetime of the MOT can be much longer than the blue Sr MOT, up to several seconds in Ref. Kuwamoto99 . The contrast of the forbidden line would then be close to 100%.

The observation of the $`{}_{}{}^{1}S_{0}^{}^3P_0`$ transition is a first step towards the realization of an optical frequency standard using trapped neutral atoms. Clearly the next step is the measurement of the wavelength of the dipole trap beam where light-shift cancellation occurs. This measurement requires better frequency resolution than achieved here, i.e. higher laser intensity at 698 nm and/or colder atoms.

In this new type of optical frequency standards a line-$`Q`$ as high as $`10^{15}`$ is achievable if the spontaneous emission rate in the optical trap is less than one per second. With a laser intensity of $`2\times 10^7`$ W/m2, the trap oscillation frequency can be 50 kHz and the spontaneous emission rate 0.6 s-1. With a line-$`Q`$ of $`10^{15}`$ and a reasonable trapped atom number of the order of $`10^6`$, ultimate performances are orders of magnitude better than existing devices. The frequency noise of the laser used to probe the atoms will then be of decisive importance Quessada03 . State of the art ultra stable lasers use macroscopic resonators as a reference and exhibit $`1/f`$ frequency noise at Fourier frequencies below a few Hz Young99 . We propose to circumvent this problem with a first stage servo control to the atomic transition with a separate setup fully optimized for frequency stability: large number of atoms, favorable duty cycle…At the price of an increase of the experimental width of the atomic resonance the response time of the servo-control can be as fast as desired. With this additional degree of freedom, the optimization of ultra-stable lasers should lead to a large improvement of their performances. One could then approach the demanding requirements of the optical standard using trapped neutral atoms.

We thank Ouali Acef, André Clairon, Michel Lours and Giorgio Santarelli for helpful discussions, and the optoelectronic group of the university of Bath (UK) for providing the photonic cristal fiber of the frequency chain. A. B. acknowledges his grant from the european Research Training Network CAUAC. BNM-SYRTE is Unité Associée au CNRS (UMR 8630).