format: text

text:
# How the No-Cloning Theorem Got its Name

## Abstract

I was the referee who approved the publication of Nick Herbert’s FLASH paper, knowing perfectly well that it was wrong. I explain why my decision was the correct one, and I briefly review the progress to which it led.

The no-cloning theorem  is of fundamental importance in quantum theory. It asserts that no quantum amplifier can duplicate accurately two or more nonorthogonal quantum states. A simple proof requires only a few lines . Why wasn’t that theorem discovered fifty years earlier? What are the events that actually led to its discovery and publication?

This is the story of my own personal contribution to the no-cloning theorem, made public for the first time after more than twenty years. Early in 1981, the editor of Foundations of Physics asked me to be a referee for a manuscript by Nick Herbert, with title “FLASH—A superluminal communicator based upon a new kind of measurement.” It was obvious to me that the paper could not be correct, because it violated the special theory of relativity. However I was sure this was also obvious to the author. Anyway, nothing in the argument had any relation to relativity, so that the error had to be elsewhere.

Herbert’s apparatus was an idealized laser gain tube which would have macroscopically distinguishable outputs when the input was a single arbitrarily polarized photon. Indeed, the word LASER is an acronym for “Light Amplification by Stimulated Emission of Radiation.” However, besides the stimulated emission, there is also spontaneous emission, which results in noise. Herbert’s claim was that the noise would not prevent identifying the polarization of the incoming photon, at least statistically. For this, he used the notion of quantum compounds  and many properties of laser physics with which I was not familiar.

I recommended to the editor of Foundations of Physics that this paper be published . I wrote that it was obviously wrong, but I expected that it would elicit considerable interest and that finding the error would lead to significant progress in our understanding of physics. Soon afterwards, Wootters and Zurek  and Dieks  published, almost simultaneously, their versions of the no-cloning theorem. The tantalizing title “A single quantum cannot be cloned” was contributed by John Wheeler. How the present paper got its name is another story .

There was another referee, GianCarlo Ghirardi, who recommended to reject Herbert’s paper. His anonymous referee’s report contained an argument which was a special case of the theorem in references . Perhaps Ghirardi thought that his objections were so obvious that they did not deserve to be published in the form of an article (he did publish them the following year ). Other objections were raised by Glauber , and then by many other authors whom I am unable to cite, because of space limitations.

With some hindsight, it is now clear that the no-cloning interdiction was implicitly used by Stephen Wiesner in his seminal paper Conjugate Coding which was submitted circa 1970 to IEEE Transactions on Information Theory, and promptly rejected because it was written in a jargon incomprehensible to computer scientists (this actually was a paper about physics, but it had been submitted to a computer science journal). Wiesner’s work was finally published in its original form in 1983  in the newsletter of ACM SIGACT (Association for Computing Machinery, Special Interest Group in Algorithms and Computation Theory). Another early article, Unforgeable Subway Tokens , also tacitly assumes that exact duplication of a quantum state is impossible. As it often happens in science, these things were well known to those who know things well.

Anyway, the no-cloning theorem was no justification for rejecting Herbert’s paper. His manuscript mentioned exact duplication as a theoretical ideal, but it was clear that what was actually produced by the laser gain tube was a messy state, where many different outputs were entangled with the laser final state.

The no-cloning theorem and its practical implications for quantum cryptography led to an enormous wave of interest, as I had predicted. Every month one or more papers on this subject are put in the quant-ph electronic archive or appear in other media. Typical questions that are discussed are how to achieve optimal, if not perfect cloning , a question which led to the important notion of quantum disentanglement ; the no-broadcasting theorem , which is a generalization of no-cloning to general impure density matrices; $`MN`$ probabilistic cloning ; rigorous and stronger theorems ; and in particular the relation of optimal cloning to the impossibility of signaling  and bounds on state estimation that are imposed by the no signaling condition . The problem was finally laid to rest by a fundamental result: completely positive maps, the only ones that can be realized by quantum mechanical systems, cannot increase distinguishability .

There is no doubt that the term “superluminal” was one of the causes this subject became so attractive. Who would not be happy to beat the relativistic limit on the speed of transmission of information? Actually, superluminal group velocities have been observed in barrier tunneling in condensed matter . However, special relativity does not forbid the group velocity to exceed $`c`$. It is the front velocity of a wave packet that is the relevant criterion for signal transmission, and the front velocity never exceeds $`c`$.

Clearly, the term superluminal has no room in the present considerations. The speed of light never enters in them, nor does the theory of relativity. The issue simply is: given two observers with an entangled bipartite quantum system (or even numerous such systems: they are equivalent to a single, bigger, physical system) it is impossible, by means of local quantum operations (LQO), to transmit any information whatsoever from one observer to the other, without transmitting real material objects between them.

This is what quantum mechanics says , and a speculative question could be: what kinds of modifications of quantum mechanics would allow information transfer by means of LQOs. In the following, I’ll show that a nonlinear evolution of the density matrix $`\rho `$ would give such a result. A similar conclusion was also reached by Svetlichny , who however only considered the special case of an initial maximally entangled pure state, and made further restrictive assumptions (and also used the term “superluminal” without ever using the theory of relativity).

Consider our two familiar observers, Alice and Bob, who have a bipartite quantum system in a known state $`\rho `$. They perform LQOs that are mathematically represented by positive operator valued measures (POVMs) with elements $`A_j=\text{1l}`$ and $`B_\mu =\text{1l}`$ respectively. The probability for the joint result $`j\mu `$ is

$$P_{j\mu }=\mathrm{Tr}_A\mathrm{Tr}_B(\rho A_jB_\mu ),$$
(1)
where the double trace is taken on the indices of Alice and of Bob. If Bob is not informed that Alice got result $`j`$, the probability that he gets $`\mu `$ is

$$\underset{j}{}P_{j\mu }=\mathrm{Tr}(\rho _BB_\mu ),$$
(2)
where $`\rho _B=\mathrm{Tr}_A(\rho )`$ is Bob’s reduced density matrix. This result does not depend on Alice’s choice of a POVM. This is what quantum mechanics says, and all this is well known.

If Bob is informed that Alice got result $`j`$, it follows from (1) that the probability that Bob has result $`\mu `$ is $`\mathrm{Tr}_B(\rho _jB_\mu )`$, where

$$\rho _j=\mathrm{Tr}_A(\rho A_j)/p_j.$$
(3)
Here the denominator $`p_j=_\nu P_{j\nu }`$ is the probability that Alice gets result $`j`$, so that $`\rho _j`$ has unit trace as it should.

Therefore everything happens as if the state of Bob’s subsystem actually was $`\rho _j`$. Quantum mechanics does not claim that this is true, but also gives no way of showing that this realistic point of view is false. It’s just a matter of belief. Bernard d’Espagnat calls Bob’s reduced density matrix an “improper” mixture . I also considered such situations in , where I distinguished pure states, mixtures, and “compounds” (a compound is a mixture that has a unique decomposition into pure states if additional information is supplied). Here, I am slightly more general: Bob’s reduced density matrix can be split in a unique way into other density matrices (not only pure states) if Alice reveals to him which result she got.

The question is whether Bob can do that without Alice’s help, and thereby find out which POVM she chose to perform. For this, one has to violate quantum mechanics in some way, such as cloning, as Nick Herbert originally proposed. Here, it is essential to assume that cloning an improper mixture (a “compound”) means to clone each component separately, as if the state really was one of these components, with probability $`p_j`$, and we merely ignored which one.

Let us suppose that Bob can perform a trace-preserving nonlinear transformation

$$\rho _j\stackrel{~}{\rho }_j,$$
(4)
on the unknown “true” states of his subsystem. Then his reduced density matrix evolves as

$$\rho _B=\underset{j}{}p_j\rho _j\underset{j}{}p_j\stackrel{~}{\rho }_j.$$
(5)
Note that the result is not determined by $`\rho _B`$ alone, but the decomposition of $`\rho _B`$ into a definite set of $`\rho _j`$ is essential. In particular, the result depends on the choice made by Alice of a particular POVM, as Eq. (3) explicitly shows.

If Alice chooses a different POVM, $`A_s=\text{1l}`$, giving Bob states $`=\mathrm{Tr}_A(\rho A_s)/p_s`$ with probability $`p_s`$, Bob would obtain, after the hypothetical nonlinear evolution,

$$\rho _B=\underset{s}{}p_s\rho _s\underset{s}{}p_s\stackrel{~}{\rho }_s.$$
(6)
Let us call the right hand sides of the last two equations $`\rho ^{}`$ and $`\rho ^{\prime \prime }`$, respectively. In the generic case, they are not equal to each other and they can be distinguished statistically. If enough copies are supplied to Bob, he would be able to know which POVM Alice chose.

Needless to say, the assumption of a nonlinear evolution of $`\rho `$ violates quantum mechanics. Contrarywise, if $`\rho `$ evolves linearly, then $`\rho ^{}=\rho ^{\prime \prime }`$ and no information can be transmitted in this way. Note that all the discussion is about the density matrix $`\rho `$. It is quite possible to have a nonlinear evolution of pure states which corresponds to a linear evolution of $`\rho `$. For example, the evolution

$$\left(\genfrac{}{}{0pt}{}{\alpha }{\beta }\right)\left(\genfrac{}{}{0pt}{}{1}{0}\right)$$
(7)
is nonlinear in terms of pure states, because

$$\left(\genfrac{}{}{0pt}{}{\alpha }{\beta }\right)\alpha \left(\genfrac{}{}{0pt}{}{1}{0}\right)+\beta \left(\genfrac{}{}{0pt}{}{0}{1}\right),$$
(8)
and the right hand side would give $`(\alpha +\beta )\left(\genfrac{}{}{0pt}{}{1}{0}\right)`$ if the evolution were linear.

On the other hand, the same evolution expressed in terms of density matrices appears as

$$\left(\genfrac{}{}{0pt}{}{ab}{b^{}1a}\right)\left(\genfrac{}{}{0pt}{}{10}{00}\right).$$
(9)
This is a linear operation, generated by a pair of Kraus matrices

$$\left(\genfrac{}{}{0pt}{}{10}{00}\right)\mathrm{and}\left(\genfrac{}{}{0pt}{}{01}{00}\right).$$
(10)
In summary, Nick Herbert’s erroneous paper was a spark that generated immense progress. There also are many wrong papers that have been published in reputable journals, some of them by renowned scientists. Their bad influence may last for years. For these, I decline all responsibility. I was not the referee of these papers and I could not protect the good reputation of their authors.

I am grateful to Erika Andersson, Gilles Brassard, Dagmar Bruß, Chris Fuchs, GianCarlo Ghirardi, Nobu Imoto, Danny Terno, and Bill Wootters for many clarifying comments and for helping me to locate some of the references below. This work was supported by the Gerard Swope Fund and the Fund for Encouragement of Research.


text:
# A note on ”A Matrix Realignment Method for Recognizing Entanglement,” quant-ph/0205017 v1

In  Chen, Wu and Yang formulated a necessary separability criterion based on a realignment method for matrices. This note is to point out that this criterion is identical to the necessary cross norm criterion put forward in . The cross norm criterion in  is based on a canonical isomorphism between Hilbert-Schmidt operators acting on $`^n^n`$ and Hilbert-Schmidt operators mapping $`^{n^2}`$ into itself. The matrix realignment method proposed in  is this canonical isomorphism. Therefore the criteria formulated in  and  and subsequent results, including examples, are identical. Moreover, from the formulation of the criterion given in  the statement of Proposition 1 in  follows as an immediate corollary and it is also immediately obvious that the dual criterion formulated in Theorem 2 in  is equivalent to the original criterion.


text:
# Proof of Bose-Einstein Condensation for Interacting Gases with a One-Particle Spectral Gap

## 1 Introduction

The long-standing problem of proving the existence of a Bose-Einstein condensation (BEC) in non-ideal real Bose gases with a standard two-body interaction has recently been warmed up by the great success of observing this phenomenon in trapped gases. Notice that much earlier BEC has been observed in liquid $`{}_{}{}^{4}\mathrm{He}`$, which however remained always under discussion.

In the present note we announce our result together with a sketch of the proof about the existence of the standard or zero-mode BEC in Bose gases with realistic superstable two-body interactions and with gap in the one-particle excitation energy spectrum. To the best of our knowledge, this is a first proof of its sort for homogeneous systems. We are not using any scaling limits, (e.g. type van der Waals limits ) or truncation of particle interactions . We prove that BEC occurs by constructing a positive lower bound for the condensate density which is valid for low enough temperature and appropriate large density of particles.

We consider a gas of interacting Bosons in cubic boxes $`\mathrm{\Lambda }=L\times L\times L^3`$ with periodic boundary conditions. We look here in detail at the $`3`$ dimensional case, and comment later on the case of other dimensions. Denote by $`V=L^3`$ the volume of the box $`\mathrm{\Lambda }`$. The grand-canonical Hamiltonian of this system reads

$$H_{\mathrm{\Lambda },g}^\mathrm{\Delta }=T_\mathrm{\Lambda }^\mathrm{\Delta }\mu N_\mathrm{\Lambda }+gU_\mathrm{\Lambda },g>0,$$
(1)
where $`T_\mathrm{\Lambda }^\mathrm{\Delta }`$ is the kinetic energy with gap $`\mathrm{\Delta }>0`$ in its spectrum,

$$T_\mathrm{\Lambda }^\mathrm{\Delta }=\underset{k\mathrm{\Lambda }^{}}{}\frac{\mathrm{}^2k^2}{2m}a_k^{}a_k\mathrm{\Delta }a_0^{}a_0.$$
(2)
The sum $`k`$ runs over the set $`\mathrm{\Lambda }^{}`$, dual to $`\mathrm{\Lambda }`$, i.e.

$$\mathrm{\Lambda }^{}=\left\{k^3;k_\alpha =2\pi n_\alpha /L;n_\alpha =0,\pm 1,\mathrm{};\alpha =1,2,3\right\}.$$

The operators $`a_k^{}`$ and $`a_k`$ are the Bose creation and annihilation operators for mode $`k`$. As usual, the mode-occupation number operators are denoted by $`N_k=a_k^{}a_k`$, $`N_\mathrm{\Lambda }=_{k\mathrm{\Lambda }^{}}N_k`$ is the total number operator in $`\mathrm{\Lambda }`$.

We assume a priori the presence of a gap $`\mathrm{\Delta }`$ in the one-particle excitations spectrum, isolating the lowest (zero-mode) energy level. This can be realised by taking appropriate boundary conditions, attractive boundary conditions , or such a gap can also be realised by specific interparticle interactions effectively incorporated in general two-body interactions .

Of course $`\mu `$ is the chemical potential and the interaction between the particles is modelled by the two-body interaction term

$$U_\mathrm{\Lambda }=\frac{1}{2}_{\mathrm{\Lambda }^2}dxdya^{}(x)a^{}(y)v(xy)a(y)a(x),$$
(3)
where $`a^{}(x),a^{}(y)`$ and $`a(y),a(x)`$ are the creation and annihilation operators for the Bose particles at $`x,y^3`$. The interaction potential $`v`$ is assumed to be spherically symmetric, superstable , i.e. it satisfies the inequality

$$\underset{1i<jn}{}v(x_ix_j)\frac{A}{2V}n^2Bn$$
(4)
for some constants $`A>0`$, $`B0`$, and all $`n2`$, $`x_i\mathrm{\Lambda }`$. Consequently, the interaction term (3) satisfies

$$U_\mathrm{\Lambda }\frac{A}{2V}N_\mathrm{\Lambda }^2BN_\mathrm{\Lambda }.$$
(5)
This superstability property is, together with the spectral gap (2), the physical foundation of our proof. Intuitively, one might understand that condensation in the groundstate (i.e. $`k=0`$), which is energetically isolated by a gap $`\mathrm{\Delta }`$ can survive the switching-on of a gentle interaction, and that fluctuations must be of a macroscopical size to overcome this gap and lift particles out of the isolated groundstate. More technical ingredients of the proof are the convexity properties of thermodynamical potentials, such as the pressure, and the use of an optimal choice for the constants $`A,B`$ in the superstability criterion (4). Indeed, it was proved  that continuous $`L^1`$-functions of positive type $`v:^3`$ are superstable potentials if and only if

$$\widehat{v}(0)\widehat{v}(q)=_^3dxv(x)\mathrm{e}^{\mathrm{i}qx}0,q^3,$$
(6)
and $`\widehat{v}(0)>0`$. Moreover, Lewis, Pulè, and de Smedt  proved the existence of the optimal constants $`A=\widehat{v}(0)(1ϵ)`$ and $`B=v(0)/2`$ in (4) for this type of potentials. Here, $`ϵ>0`$ is an arbitrarily small positive constant. It is related to the size of the system, and can be put to zero after the thermodynamic limit. This optimal choice is of determining importance in our proof of the zero-mode Bose-Einstein condensation.

## 2 Sketch of Proof

The main idea of the proof is to estimate the Bose condensate of the full model (1) by the condensate of a particularly chosen reference system for which one knows the occurrence of condensation. The clever choice of this reference system is the subtle point of our proof. The reference system is a so-called mean-field Bose gas, an exactly solvable model of Bosons  (for a review see ) defined by the grand-canonical Hamiltonian as

$$H_{\mathrm{\Lambda },g,\lambda }^\mathrm{\Delta }=T_\mathrm{\Lambda }^\mathrm{\Delta }\mu N_\mathrm{\Lambda }+g\frac{\lambda }{2V}N_\mathrm{\Lambda }^2.$$
(7)
The kinetic energy operator $`T_\mathrm{\Lambda }^\mathrm{\Delta }`$ (2) is as for the general interacting system (1), but the interaction term (3) is replaced by a mean-field interaction term.

The reference system (7), a mean-field Bose gas, emerges as the van der Waals limit of the fully interacting system . In that case, the constant $`\lambda `$ equals $`\widehat{v}(0)`$ (6), which means that the van der Waals limit reduces the full interaction to the $`\widehat{v}(0)`$ contribution. In our proof we tune the constant $`\lambda `$ in order to get the best possible lower bound for the condensate density of the full interaction model.

Remark that our reference system (7) does show Bose condensation for large enough densities (i.e. for $`\mu `$ large) at any given temperature. Moreover, systems of the type of our reference system have better properties than the ideal Bose gas which is in many ways a pathological model, e.g. in the sense that there is no equivalence of ensembles  and in the sense that the chemical potential is limited by a zero upperbound in order to safeguard thermodynamical stability. These are the reasons for our strategy of using the free of those pathologies reference system (7).

Thermodynamic properties of the reference system (7). It is well known that our reference system (7) is a soluble model. In the case of vanishing gap $`\mathrm{\Delta }=0`$, the complete solution can be found at several places in the literature . In particular there is condensation for all dimensions $`D3`$ at any temperature for densities large enough. It is a student exercise to work out now the case with gap $`\mathrm{\Delta }>0`$. There is one main difference with the gapless case, namely the presence of the gap provokes a shift in the chemical potential and its treshholds, and one gets condensation in all dimensions $`D1`$. One derives straightforwardly that for the reference model (7) one gets condensation for all values of the chemical potential satisfying,

$$\mu >g\lambda \rho ^P(\beta ,\mathrm{\Delta })\mathrm{\Delta }.$$

where $`\rho ^P(\beta ,\mathrm{\Delta })`$ is the total particle density of the perfect Bose Gas (PBG) at inverse temperature $`\beta `$ and chemical potential equal to $`\mathrm{\Delta }`$. Moreover, the condensate density

$$\rho _{0,g,\lambda }^\mathrm{\Delta }(\beta ,\mu )=\underset{V\mathrm{}}{lim}\frac{1}{V}N_0_{H_{\mathrm{\Lambda },g,\lambda }^\mathrm{\Delta }}(\beta ,\mu ),$$

i.e. the particle density at the zero mode in the thermodynamic limit $`(V\mathrm{})`$ of the grand-canonical Gibbs states $`_{H_\mathrm{\Lambda }}(\beta ,\mu )`$, in volumes $`\mathrm{\Lambda }`$ for a certain choice of inverse temperature and chemical potential $`(\beta ,\mu )`$ and Hamiltonian $`H_\mathrm{\Lambda }`$, is explicitly given by

$$\rho _{0,g,\lambda }^\mathrm{\Delta }=\frac{\mu +\mathrm{\Delta }}{g\lambda }\rho ^P(\beta ,\mathrm{\Delta }).$$
(8)
One also computes the total particle density for given $`(\beta ,\mu )`$ to be

$$\rho _{g,\lambda }^\mathrm{\Delta }=\underset{V\mathrm{}}{lim}\frac{1}{V}N_\mathrm{\Lambda }_{H_{\mathrm{\Lambda },g,\lambda }^\mathrm{\Delta }}=\frac{\mu +\mathrm{\Delta }}{g\lambda }.$$
(9)
To prove Bose condensation in the full model (1), we subtract from (3) the long-range part of the interaction proportional to $`N_\mathrm{\Lambda }^2/2V`$, tune it with a factor, taking into account the optimal stability constants, and add it to the kinetic-energy term. The latter serves as our reference system (7), from which we establish a lower bound on the zero-mode condensate density $`\rho _{0,g}^\mathrm{\Delta }(\beta ,\mu )=lim_V\mathrm{}N_0/V_{H_{\mathrm{\Lambda },g}^\mathrm{\Delta }}`$ in the fully interacting system (1). In the lemma below, a lower bound on $`\rho _{0,g}^\mathrm{\Delta }(\beta ,\mu )`$ is given.

###### Lemma 1

The zero-mode condensate density $`\rho _{0,g}^\mathrm{\Delta }(\beta ,\mu )`$ in the thermodynamic limit of grand-canonical Gibbs states of interacting systems (1) with superstable two-body potentials $`v`$ (6), has the following lower bound:

$`\rho _{0,g}^\mathrm{\Delta }(\beta ,\mu )`$ $``$ $`{\displaystyle \frac{\mu }{g\widehat{v}(0)}}+{\displaystyle \frac{g\widehat{v}(0)}{2\mathrm{\Delta }}}\left(\rho ^P(\beta ,\mathrm{\Delta })\right)^2{\displaystyle \frac{\mu +\mathrm{\Delta }}{\mathrm{\Delta }}}\rho ^P(\beta ,\mathrm{\Delta })`$ (10)
$`{\displaystyle \frac{gv(0)}{2\mathrm{\Delta }}}\rho _g^{(\mathrm{\Delta }=0)}(\beta ,\mu )\rho _c^P(\beta ).`$
Here $`\rho _g^{(\mathrm{\Delta }=0)}(\beta ,\mu )`$ is the total density of the interacting gas without gap (1). $`\rho ^P(\beta ,\mathrm{\Delta })`$ refers to the total density of the PBG at inverse temperature $`\beta `$ and chemical potential $`\mu =\mathrm{\Delta }`$, $`\rho _c^P(\beta )`$ is the critical density of the PBG. The bound is valid for values $`\mu >g\widehat{v}(0)\rho _c^P(\beta )`$, and dimensions $`D3`$.

Idea of the proof. Using the Bogoliubov convexity inequality  one gets:

$$\frac{g}{V}W_\mathrm{\Lambda }^\lambda _{H_{\mathrm{\Lambda },g}^\mathrm{\Delta }}p_\mathrm{\Lambda }[H_{\mathrm{\Lambda },g,\lambda }^\mathrm{\Delta }]p_\mathrm{\Lambda }[H_{\mathrm{\Lambda },g}^\mathrm{\Delta }]\frac{g}{V}W_\mathrm{\Lambda }^\lambda _{H_{\mathrm{\Lambda },g,\lambda }^\mathrm{\Delta }}.$$
(11)
This gives upper and lower bounds on the difference of the grand-canonical pressure $`p_\mathrm{\Lambda }[H_\mathrm{\Lambda }]`$ of the mean-field reference Bose gas (7) and the full model (1). The operator $`W_\mathrm{\Lambda }^\lambda `$ is the difference between the interactions of the fully interacting and the mean-field Bose gases, i.e. $`W_\mathrm{\Lambda }^\lambda =U_\mathrm{\Lambda }\frac{\lambda }{2V}N_\mathrm{\Lambda }^2`$. The expectation values in (11) can be estimated using, for the lower bound, the superstability properties of the interaction, and, for the upper bound, the properties of the equilibrium states of the mean-field reference Bose gas. The lower bound in (11) follows from (5), and from the tuning of the interaction parameter $`\lambda `$ for the mean-field reference Bose gas (7) to the constant $`A`$ in (5),

$$\frac{g}{V}W_\mathrm{\Lambda }^A_{H_{\mathrm{\Lambda },g}^\mathrm{\Delta }}\frac{gB}{V}N_\mathrm{\Lambda }_{H_{\mathrm{\Lambda },g}^\mathrm{\Delta }}.$$
(12)
On the other hand, using the mode by mode gauge invariance of the Gibbs states of the mean-field Bose gas we arrive at the following upper bound for the pressure difference (11),

$$\frac{g}{V}W_\mathrm{\Lambda }^A_{H_{\mathrm{\Lambda },g,A}^\mathrm{\Delta }}\frac{g}{V^2}CN_\mathrm{\Lambda }^2\frac{\widehat{v}(0)}{2}N_0^2_{H_{\mathrm{\Lambda },g,A}^\mathrm{\Delta }},$$
(13)
where $`C=\widehat{v}(0)A/2`$. It follows from the properties of the mean-field Bose gas that the expectation values in the rhs of (13) in the limit $`(V\mathrm{})`$ are given by $`gC(\rho _{g,A}^\mathrm{\Delta }(\beta ,\mu ))^2g\widehat{v}(0)(\rho _{0,g,A}^\mathrm{\Delta }(\beta ,\mu ))^2/2`$.

The pressure $`p_\mathrm{\Lambda }[H_\mathrm{\Lambda }^\mathrm{\Delta }]`$ is an increasing convex function of $`\mathrm{\Delta }0`$. Since the condensate density $`\rho _{0,g}^\mathrm{\Delta }(\beta ,\mu )`$ is the derivative of the pressure with respect to $`\mathrm{\Delta }`$, by convexity we find for it a lower bound, given by

$$\frac{1}{V}N_0_{H_{\mathrm{\Lambda },g}^\mathrm{\Delta }}\frac{p_\mathrm{\Lambda }[H_{\mathrm{\Lambda },g}^\mathrm{\Delta }]p_\mathrm{\Lambda }[H_{\mathrm{\Lambda },g}^{(\mathrm{\Delta }=0)}]}{\mathrm{\Delta }}.$$
(14)
Analogously, by virtue of the same convexity property, the difference of the pressures between the mean-field Bose gas with gap and without gap, is bounded from below by the condensate density for the mean-field gas without gap, i.e.

$$\frac{p_\mathrm{\Lambda }[H_{\mathrm{\Lambda },g,A}^\mathrm{\Delta }]p_\mathrm{\Lambda }[H_{\mathrm{\Lambda },g,A}^{(\mathrm{\Delta }=0)}]}{\mathrm{\Delta }}\frac{1}{V}N_0_{H_{\mathrm{\Lambda },g,A}^{(\mathrm{\Delta }=0)}}.$$
(15)
Adding up these two inequalities and using the Bogoliubov convexity inequality (11) once at $`\mathrm{\Delta }>0`$ and once at $`\mathrm{\Delta }=0`$, together with the bounds (12) and (13), one gets the following lower bound for the condensate density $`\rho _{0,g}^\mathrm{\Delta }(\beta ,\mu )`$ in the thermodynamic limit ($`V\mathrm{}`$):

$`\rho _{0,g}^\mathrm{\Delta }(\beta ,\mu )`$ $``$ $`\rho _{0,g,A}^{(\mathrm{\Delta }=0)}(\beta ,\mu )+g{\displaystyle \frac{\widehat{v}(0)}{2\mathrm{\Delta }}}(\rho _{0,g,A}^\mathrm{\Delta }(\beta ,\mu ))^2`$ (16)
$`{\displaystyle \frac{g}{\mathrm{\Delta }}}\left(B\rho _g^{(\mathrm{\Delta }=0)}(\beta ,\mu )+C(\rho _{g,A}^\mathrm{\Delta }(\beta ,\mu ))^2\right).`$
The lower bound (10) now follows from (16) by use of the explicit expressions for total density and the condensate density of the mean-field Bose gas with gap, (9) and (8), and by the well-known expression for the condensate density in the gapless mean-field model:

$$\rho _{0,g,\lambda }^{(\mathrm{\Delta }=0)}=\underset{V\mathrm{}}{lim}\frac{1}{V}N_0_{H_{\mathrm{\Lambda },g,\lambda }^{(\mathrm{\Delta }=0)}}=\frac{\mu }{g\lambda }\rho _c^P(\beta ),$$

if $`\mu >g\lambda \rho _c^P(\beta )`$, where $`\rho _c^P(\beta )`$ is the critical density for the PBG at inverse temperature $`\beta `$. As a last step, we need the optimal superstability constants for continuous $`L^1`$ potentials of positive type , i.e. we take $`A=\widehat{v}(0)`$, and $`B=v(0)/2`$ after the thermodynamic limit, and get the expression for the lower bound (10) in the lemma. QED

Now we get our main result:

###### Theorem 2

Consider a three dimensional system of interacting Bose particles (1), with a superstable two-body potential. There exists a minimal gap $`\mathrm{\Delta }_{min}`$ such that for any finite $`\mathrm{\Delta }\mathrm{\Delta }_{min}`$, one has $`\rho _{0,g}^\mathrm{\Delta }(\beta ,\mu )>0`$ or zero-mode condensation.

Proof. Take any $`\eta >0`$, fix a temperature and a chemical potential $`(\beta ,\mu )`$ such that $`\mu >g\widehat{v}(0)(\rho _c^P(\beta )+3\eta )`$, then there exists a minimal gap $`\mathrm{\Delta }_{min}`$ such that for any finite $`\mathrm{\Delta }\mathrm{\Delta }_{min}`$,

$$\left|\frac{g\widehat{v}(0)}{2\mathrm{\Delta }}\left(\rho ^P(\beta ,\mathrm{\Delta })\right)^2\frac{gv(0)}{2\mathrm{\Delta }}\rho _g^{(\mathrm{\Delta }=0)}(\beta ,\mu )2\frac{\mu +\mathrm{\Delta }}{\mathrm{\Delta }}\rho ^P(\beta ,\mathrm{\Delta })\right|<\eta $$

For these values of the gap $`\mathrm{\Delta }`$ we have $`\rho _{0,g}^\mathrm{\Delta }(\beta ,\mu )>\eta >0`$, by virtue of the lower bound (10) in Lemma 1, and hence we proved condensation. QED

## 3 Discussion

First of all let us remark that our proofs hold without any change in all dimensions $`D3`$. For $`D=1`$ or $`2`$ a similar lower bound on the condensate density can be derived, on the basis of modified convexity arguments (14)–(15), i.e. one has to consider pressure differences of the form $`p_\mathrm{\Lambda }[H_\mathrm{\Lambda }^\mathrm{\Delta }]p_\mathrm{\Lambda }[H_\mathrm{\Lambda }^{\mathrm{\Delta }_0}]`$, with $`0<\mathrm{\Delta }_0<\mathrm{\Delta }`$, instead of with $`\mathrm{\Delta }_0=0`$ in (14)–(15). This yields the substitution of $`\rho _{0,g,A}^{\mathrm{\Delta }_0}(\beta ,\mu )`$ and $`\rho _g^{\mathrm{\Delta }_0}(\beta ,\mu )`$ for $`\rho _{0,g,A}^{(\mathrm{\Delta }=0)}(\beta ,\mu )`$ and $`\rho _g^{(\mathrm{\Delta }=0)}(\beta ,\mu )`$ in (16). Hence, also in one and two dimensional interacting Bose gases with gap (1), Bose condensation is proved, in contrast to the Bogoliubov–Hohenberg theorem  which yields the absence of BEC for $`1D`$ or $`2D`$ translation invariant continuous Bose systems without gap.

Finally we want to stress that our results are for homogeneous systems. On the other hand we have to mention here the interesting recent result about the condensation for trapped Bose gases , i.e. for inhomogeneous systems where a rigorous proof is given of Bose condensation in the so-called Gross-Pitaevskii limit. As such trapped systems always have a gap in the one-particle spectrum, we consider our result as a bridge between the homogeneous gapless systems on one hand and the trapped systems on the other hand. Moreover, we hope that our work might be inspiring to establish a proof of BEC for homogeneous systems when the gap tends to zero.

J. L. gratefully acknowledges financial support from K.U.Leuven grant FLOF-10408, and V. A. Z. acknowledges ITF K.U.Leuven for hospitality.


text:
# New nearby stars selected in a high proper motion survey by DENIS photometry

## 1 Introduction

The solar neighbourhood serves as a fundamental constrain for our understanding of the stellar physics and the Galaxy. The nearest stars provide accurate data on luminosities, temperatures, masses, which are fundamental parameters of stellar astronomy. Their kinematics, chemical composition, and mass function also hold important clues to the nature and history of the Milky Way.

However, among the nearest stars, many of the low luminosity stars such as white, red, and brown dwarfs remain undetected. Henry et al. (1997) estimated that more than 30 % of the stars within 10 pc are currently missing from the solar neighbourhood sample. The deficit is estimated to be twice within 25 pc (Henry et al., 2002) and is largest south of decination -30. As a confirmation, stars are continuously identified as nearby stars in proper motion catalogues or in Schmidt plates (Henry et al., 1997, 2002; Gizis & Reid, 1997; Scholz et al., 1999, 2001, 2002; Jahreiss et al., 2001; Phan-Bao et al., 2001; Delfosse et al., 2001).

In addition, new nearby faint stars, mainly brown dwarfs, are detected in deep sky surveys (Ruiz et al., 1997; Delfosse et al., 1997, 1999; Martín et al., 1999; Kirkpatrick et al., 1999, 2000; Reid et al., 2000; Fan et al., 2000; Gizis et al., 2000).

One way to detect nearby stars is to select the red faint objects obtained in near-infrared surveys: the Two Micron All Sky Survey (2MASS) (Kirkpatrick et al., 1999, 2000; Reid et al., 2000; Gizis et al., 2000) and the Deep Near-Infrared Survey of the Southern Sky (DENIS) (Delfosse et al., 2001; Phan-Bao et al., 2001). Most of the nearest stars found in these surveys turned out to be high proper motion objects. Given this fact, Scholz et al. (2001) searched for new stars in the solar neighbourhood by combining proper motion catalogues with the 2MASS data base.

Using the same approach, we cross-identify a proper motion catalogue with DENIS data and estimate the distances using near-infrared photometry. The proper motion catalogue is described in Sect. 2. In Sect. 3, we give the result of the cross-identification with the DENIS data base. Sect. 4 presents the distance estimation. In Sect. 5, we discuss the effect of metallicity on the distance estimation. In Sect. 6, we list the stars found to be within 50 pc and discuss the precision on our distance determinations.

## 2 High proper motion sample

One of the criteria for finding nearby stars is their large proper motion. The main source for the identification of such stars are the high-proper motion catalogues of Luyten: the Luyten Half Second proper motion catalogue (LHS) and the New Luyten Catalogue of Stars with Proper Motions Larger than Two Tenths of an Arcsecond (NLTT) based on observations with the Palomar Schmidt telescope. Dawson (1986) showed that the LHS catalogue is incomplete south of declination -33, $`R>18`$ or $`\mu >2.5^{\prime \prime }`$ yr-1. Wroblewski & Torres (1989, 1991, 1994, 1996, 1997) improved the completness in the southern sky and found new stars with $`\mu >0.15^{\prime \prime }`$ yr-1 and photographic magnitudes $`m_{pg}<20`$. Ruiz et al. (1993) found proper motion stars in the ESO Schmidt plates with $`\mu >0.1^{\prime \prime }`$ yr-1 and $`R<20`$.

In order to complete the existing proper motion catalogues at fainter magnitudes in the southern hemisphere, Scholz et al. (2000) used Automatic Plate Measuring (APM) measurements of sky survey plates taken with the UK Schmidt telescope. Plates in the passbands $`B_J`$ and $`R`$ were taken with typical epoch difference of 15 years. The resulting catalogue contains 693 stars (out of which 195 have already been published by Scholz et al. (2000)) south of declination -15, in the regions 0h-7h, 10h-14h, and 19h-23h in right ascension. The overall search area is about 3000 square degrees. The lower proper motion limit of the survey is typically 0.25 arcsec yr-1. But in fields with smaller epoch difference between the $`B_J`$ and $`R`$ plates, only objects with even larger proper motions could be detected, respectively. The distribution of the epoch differences as a function of declination was shown in Scholz et al. (2000). The largest proper motions found are about 1.3 arcsec yr-1. The practical limit for searching for high proper motion objects using two plates in different passbands is about 1 magnitude above the respective plate limits, i.e. $`R=20.0`$ and $`B_J=21.5`$. For a comparison with other high proper motion catalogues and an estimate of the completeness of the survey, see Scholz et al. (2000). About half of the objects detected are known proper motion stars, whereas the remaining objects mainly at fainter magnitudes are new high proper motion stars.

With these measurements, they were able to fill the gap in the high proper motion survey of faint stars with $`\delta <`$ -33. In addition they obtained more accurate positions and photometry for the already known high proper motion stars. The accuracies are 0.1 to 0.2 ′′ in position, 0.03 ′′ yr-1 on proper motion, and $``$ 0.2 mag on photometry.

## 3 Cross-identification with DENIS

The DENIS survey (Epchtein et al., 1997) provides a full coverage of the southern sky in the optical band $`I`$ (0.85 $`\mu `$m) and the near-infrared bands $`J`$ (1.25 $`\mu `$m) and $`K_s`$ (2.17 $`\mu `$m). The position accuracy is 0.5 ′′ and the photometric accuracy is better than 0.1 mag. About 40% of the data are up to now calibrated at the Paris Data Analysis Center (PDAC).

We look for counterparts of the given proper motion sample in the DENIS point source catalogue extracted at the PDAC with a search radius of $`3^{\prime \prime }`$. We get 301 recoveries among the 693 objects. None of them have a counterpart in the USNO A2.0 catalogue, as expected for proper motion stars. The density of stars with no USNO counterpart in a field in the Galactic plane is less than 0.01 star within the search region. Thus the possibility of false cross-identification is low, smaller than 1%. Most of the unrecovered stars are in the unprocessed part of the DENIS survey. Dots in Fig. 1 shows the ($`I`$,$`IJ`$) colour-magnitude diagram of the recovered stars. Most stars have $`IJ>`$ 1. They are M-dwarfs (Leggett, 1992). Stars with $`0.5<IJ<1`$ are G or K-dwarfs. 1 star, with $`IJ>`$ 3, is a possible brown dwarf (square). Stars with $`IJ<`$ 0.5 are either white dwarfs (circles) or distant subdwarfs (triangles), as discussed in Sect. 4.

## 4 Photometric distance estimation

Photometric distances are estimated using the $`IJ`$ colour index. Theoretical colour-magnitude relations ($`M_I`$,$`IJ`$) for stars placed at 10 pc are superimposed to the observations in Fig 1. The solid line shows the relation for M-dwarfs with solar metallicity (Baraffe et al., 1998), the dashed line shows the relation for solar metallicity G and K-dwarfs (Lejeune et al., 1997), and the dotted line shows the relation for white dwarfs with 0.6 M(Bergeron et al.,, 1995).

We compute photometric distances from the DENIS $`IJ`$ colour index and the corresponding theoretical colour-magnitude relation: Baraffe et al. (1998) relation for M-dwarfs with $`IJ>`$ 1, Lejeune et al. (1997) relation for G and K-dwarfs in the range 0.5 $`<IJ<`$ 1. Stars with $`IJ`$ $`<`$ 0.5 can be either nearby white dwarfs or distant subdwarfs. The 5 faintest blue stars (circles in Fig. 1) are white dwarfs. Indeed if they would be subdwarfs with that faint apparent magnitude, they would be at $`d>`$ 1600 pc. Their proper motion would lead to a tangential velocity $`>`$ 4000 km s-1, much larger than the Galactic escape speed (Leonard & Tremaine, 1990; Meillon et al., 1997). For these stars, we use Bergeron et al., (1995) model atmosphere for white dwarfs of mass 0.6 M, before the turnover ($`M_I`$ 14.5), assuming they are not cool white dwarfs.

The 3 brighter blue stars (triangles in Fig. 1) are more probably late G-stars of the spheroid with large tangential velocity, rather than disc white dwarfs closer than 3 pc. The theoretical relation from Lejeune et al. (1997) for the typical spheroid metallicity $`[\frac{Fe}{H}]=1.8`$ gives $`M_I`$ 5 mag at $`IJ`$ 0.5. The derived photometric distances of these 3 stars ranges between 160 and 260 pc, and their tangential velocities between 280 and 380 km s-1, in agreement with the spheroid kinematics.

$`IK_s`$ is also a good estimator of the luminosity or effective temperature of M-dwarfs. Photometric distances differ by 15% depending on which colour index is used. However, Delfosse et al. (1997) showed that the $`K_s`$ band of DENIS can differ by $``$0.1 mag compared with the standard $`K`$ band for which models are computed. We then choose to estimate the distances from the $`IJ`$ colour index and ($`M_I`$,$`IJ`$) theoretical relations.

## 5 Metallicity effect

Within the sample magnitude range and proper motion range, one expects to get disc stars, but also thick disc and spheroid stars. Using a theoretical colour-magnitude relation at solar metallicity would overestimate the distance for the stars of these old populations. We use the Besançon model of population synthesis to reproduce the stellar content in the fields of the sample. The Besançon model has been described in (Haywood et al.,, 1997). Since then, Robin et al. (2000) improved the halo density law and mass function, while Reylé & Robin (2001) constrained the thick disc population description. New constraints available on the galactic potential, the local luminosity function and the kinematics versus age, obtained with the Hipparcos mission, are also taken into account to derive the thin disk density law (Reylé et al., to be submitted).

The simulation shows that the different populations are separated in the plane ($`H,IJ`$), where $`H=I+5\text{log}\mu +5`$ is the reduced proper motion. We define three regions in this plane where stars are most probably disc stars, thick disc stars, or spheroid stars, as represented in Fig. 2.

We determine the photometric distance using a theoretical colour-magnitude relations ($`M_I`$,$`IJ`$) at solar metallicity for stars in the disc region, at metallicity $`[\frac{Fe}{H}]=0.8`$ in the thick disc region, and at metallicity $`[\frac{Fe}{H}]=1.8`$ in the spheroid region. These relations are obtained by interpolation of the available relations with $`[\frac{Fe}{H}]=0.5,1,1.5,2`$ from Lejeune et al. (1997) for G and K-stars and Baraffe et al. (1998) (private communication for low metallicities) for M-stars.

For stars in between two regions, we make two determinations with both metallicity hypotheses (disc/thick disc or thick disc/spheroid). The probability that a star with a magnitude $`I`$, a colour index $`IJ`$ and a proper motion $`\mu _{tot}`$, belongs to one population is estimated from the population distribution of the stars in the simulated sample. In this case, both distances and tangential velocities are computed.

Disc stars with tangential velocities as large as 250 km s-1 still are found with a large probability in the simulated sample. Indeed, the sample is not representative of the entire disc population, but as a high proper motion sample, it selects the stars having a high velocity.

## 6 Stars within 50 pc

107 stars are found to have photometric distances less than 50 pc. They are mainly M-stars, but 4 stars are white dwarfs. According to the colour-reduced proper motion diagram ($`H,IJ`$), 88 stars among them belong to the disc, 14 to the thick disc stars, and 1 to the spheroid. In addition, 33 stars belong either to the disc or to the thick disc. Distance estimations with both hypothesis are made for these stars and they are possibly within 50 pc. Likewise, 7 stars belong either to the spheroid or to the thick disc and may enter or not within the 50 pc limit.

Astrometric, photometric and kinematical characteristics, estimated distances and spectral types are given in Table 1. The spectral types for M-stars are derived from the $`IJ`$ colour index (Leggett, 1992). 4 late M-type stars have been previously identified by Phan-Bao et al. (2001) and 1 by Henry et al. (2002). 6 stars already have spectroscopic distance estimates (Scholz et al., 2002). The extremely red object is the already known brown dwarf LP944-20 at 5 pc (Tinney, 1996). Trigonometric distances from the Yale parallax catalogue or the Hipparcos catalogue, or photometric distances from the ARI database for nearby stars (ARICNS), available for 25 stars, are also indicated in Table 1.

Among the list of 107 nearby stars, 15 new stars are within the 25 pc limit of the Catalogue of Nearby Stars (CNS3, Gliese & Jahreiss, 1991). New discoveries within 15 pc are LHS5045, LP225-57, LP831-45, LHS1767, and WT792. LP225-57 is presumably closer than 10 pc, at $`d`$ = 9.5 pc. Its estimated spectral type is M3.5. The star LHS124 may be a M1-star of the thick disc at 7.2 pc, although the simulation of the stellar sample tends to show it is more likely a more distant star in the disc. The distance indicated in the ARICNS (21 pc) seems also to favor a disc star.

It should be noted that the distances estimated in Table 1 have large uncertainties, particularly for stars with $`IJ<`$ 1.5, for which the colour index is less sensitive to the luminosity (see Fig. 1). Moreover, Phan-Bao et al. (2001) plotted the colour-magnitude relation of M-stars with measured distances and noted an intrinsic scatter of $`\pm `$ 1 mag. They also superimposed the Baraffe et al. (1998) theoretical track which does not fit perfectly the mean observed relation around $`IJ`$ = 1.5 where the $`M_I`$ magnitude is lower by 1 magnitude compared to the magnitude obtained with a polynomial fit of the observed data (Phan-Bao et al., 2001). For this reason, our distance determination for the 4 stars already identified by Phan-Bao et al. (2001) is smaller by a factor up to 15%.

Fig. 3 gives the difference between our estimated distance and other distances found in the litterature, as indicated in Table 1 footnotes, versus our determination. Open circles are used when the other determination is photometric, filled circles when it is trigonometric, circled dots when it is spectroscopic. Our determination is not biased compared to trigonometric parallaxes, giving confidence in our distances, while a systematic deviation appears at large distances between other photometric distances and ours, due to the different theoretical colour-magnitude relations used. The mean error is of the order of 5 pc at distances closer than 25 pc but grows to more than 10 pc at larger distances.

## 7 Conclusion

Thanks to the deep and large surveys at optical and near-infrared wavelengths, many low-luminosity nearby stars are identified. Most of the nearest stars turned out to have a large proper motion. In this paper, we present the identification of 107 stars closer than 50 pc among high proper motion stars cross-identified with DENIS. Their large proper motions exclude that they are giants. They are mainly M-stars, while 4 stars are white dwarfs. 31 stars among them have already measured distances. In addition, 40 stars may enter in the 50 pc limit depending on which population they belong to. 6 stars among them have previously measured distances.

15 new stars are within the 25 pc limit of the CNS3 catalogue. 5 stars, LHS5045, LP225-57, LP831-45, LHS1767, and WT792, are closer than 15 pc,with 1 object, LP225-57, at 9.5 pc. In order to insure the spectral types and therefore distances, we plan to obtain low-resolution spectroscopy of the closest candidates.

###### Acknowledgements.

The authors thank Guy Simon who provided photometric measures in a strip not available in the PDAC database, the whole DENIS staff and all the DENIS observers who collected the data. The DENIS project is supported by the SCIENCE and the Human Capital and Mobility plans of the European Commission under grants CT920791 and CT940627 in France, by l’Institut National des Sciences de l’Univers, the Ministère de l’éducation Nationale and the Centre National de la Recherche Scientifique (CNRS) in France, by the State of Baden-Würtemberg in Germany, by the DGICYT in Spain, by the Sterrewacht Leiden in Holland, by the Consiglio Nazionale delle Ricerche (CNR) in Italy, by the Fonds zur Förderung der wissenschaftlichen Forschung and Bundesministerium für Wissenschaft und Forschung in Austria, and by the ESO C & EE grant A-04-046.


text:
# 1 Introduction

## 1 Introduction

The potential resolution of Monte Carlo (MC) computer simulations has increased substantially over the past few years. This has been due, in part, to the dramatic rise in the performance of computers, but, more importantly, to the development of more powerful data analysis and computer simulation techniques. Histogram methods allow us to extract much more information from simulation data than was previously possible. By providing the ability to continuously vary the temperature or other intensive parameters of a simulation, these techniques have greatly simplified the analysis of simulation data by traditional means and, in addition, have also played an important role in the development of new methods of analyzing simulation data. These methods are most effective when very large numbers of spin configurations have been generated, and it is the common belief that the number needed is enlarged by correlations between successive states. More recently, a new generation of algorithms to calculate the density of states accurately via a random walk in energy space have been devised for producing canonical averages of thermodynamical quantities at essentially any temperature.  Simulation techniques have also improved immensely. Fast implementations of local update (Metropolis) algorithms have been developed for a variety of models, while cluster-flipping algorithms , which can dramatically reduce the correlation time in a simulation, now exist for several classes of models.

A different approach to increasing the performance of computer simulations is to combine several different algorithms into a single, hybrid algorithm. This idea is not new; hybrid Monte Carlo, Hybrid Molecular Dynamics, Metropolis with overrelaxation and Multi-Hit Swendsen-Wang are some examples of hybrid algorithms. In these cases, however, the two algorithms that are combined perform the simulation in different ensembles, either canonical/microcanonical or canonical/fixed-cluster-distribution. The approach we consider here is to combine algorithms that work in the same ensemble, for our examples the canonical ensemble, so that each of the individual component algorithms is a self-sufficient simulation technique. This eliminates any concerns about how the mixing of ensembles could potentially affect the quality or correctness of the results. We will, however, discuss the generalization of these “proper” hybrid algorithms to include mixed-ensemble cases.

Our aim in this work is two-fold. We first discuss, in the next section, the efficiency of a general hybrid algorithm and show how it can be improved in the case when Wolff plus Metropolis is applied to the spin-1/2 two-dimensional Ising model. Second, in section III, we apply a hybrid algorithm to the spin-3/2 two-dimensional Ising model for which a correct single cluster algorithm is not immediately obvious since the simple version does not take into account transitions between states having different spin moduli (for instance, transitions between $`\pm 3/2`$ and $`\pm 1/2`$ spin values). Further discussion and some concluding remarks are given in the last section.

## 2 General Hybrid Algorithm: spin-1/2 Ising Model

Consider a MC study of some model in which $`N`$ measurements of some observable quantity $`A`$ (energy, magnetization, susceptibility, cumulants, $`etc.`$) are made, and for which there exist several different algorithms that could be used to perform the simulation. In order to compare the efficiency of the different techniques, one needs to know both the speed with which measurements are made and the degree to which successive measurements are correlated. For this section, we will define the efficiency $`e`$ for an algorithm as

$$e=\frac{\mathrm{\#}\mathrm{of}\mathrm{measurements}\mathrm{generated}\mathrm{per}\mathrm{second}}{2\tau _A+1},$$

where the integrated autocorrelation time $`\tau _A`$ is given by

$$\tau _A=\underset{t=1}{\overset{t=N}{}}\left(1\frac{t}{N}\right)\varphi _A(t),$$

with the time-displaced correlation function $`\varphi _A(t)`$ for the quantity $`A`$ calculated as

$$\varphi _A(t)=\frac{<A(0)A(t)><A>^2}{<A^2><A>^2}.$$

Note that the correlation time, and therefore the efficiency of an algorithm, can depend strongly on the particular quantity $`A`$ measured.

Now consider a hybrid simulation algorithm that combines several different component algorithms. To set up some notation, let $`a`$ represent the number of different algorithms used, $`N_i`$ the number of measurements made with simulation technique $`i`$ and $`t_i`$ the time (in seconds) required for performing the update and making a measurement for technique $`i`$. This time will, of course, depend strongly on the implementation of the algorithms and the particular computers on which they are run.

The time in seconds needed to produce a measurement using the hybrid algorithm is

$$\frac{\underset{i=1}{\overset{a}{}}N_it_i}{\underset{i=1}{\overset{a}{}}N_i},$$

so that the efficiency of the hybrid algorithm becomes

$$e=\frac{\underset{i=1}{\overset{a}{}}N_i}{\left(2\tau _A+1\right)\times \left(\underset{i=1}{\overset{a}{}}N_it_i\right)}.$$

Please note that the correlation time, and therefore the efficiency, of the algorithm will depend on its specific implementation. For example, in a hybrid algorithm consisting of two components, $`1`$ and $`2`$, the correlation time for the sequence $`122122122122\mathrm{}`$ would, most likely, be different from the sequence $`112222112222\mathrm{}..`$ Which of the two would produce the smaller correlation time would depend on the dynamics (kinetics) of the individual algorithms.

We now demonstrate the development of the above hybrid algorithms by considering a specific example, the spin -1/2 nearest-neighbor square-lattice Ising model at its critical temperature, $`T_c`$. The Ising model has traditionally been used to test new simulation algorithms and data analysis techniques because of its simplicity and the exact finite and infinite-system solutions in the two-dimensional model. Because of the large amount of work done with the Ising model, there exist several different simulation algorithms for it. These can be broken up into two major classes: 1) single-spin update algorithms, including Metropolis, heat bath and microcanonical algorithms, and 2) cluster algorithms, including the Swendsen-Wang and Wolff algorithms. We will concentrate on two of these algorithms: Metropolis and Wolff.

Each of these algorithms has its strengths and weaknesses. The Metropolis algorithm is very efficient at equilibrating short-range fluctuations in the system, and there exist highly-optimized multi-spin coding implementations of the Metropolis algorithm. Unfortunately, the Metropolis method is not efficient at decorrelating the long-range clusters that characterize the behavior of the system near the critical point. The Wolff algorithm, on the other hand, concentrates its effort on the large clusters leading to greatly reduced correlation times and a much smaller dynamic critical exponent $`z`$. However, smaller-scale structures in the system, in particular regions of disorder, are not handled efficiently by the Wolff algorithm. The speed of the Wolff algorithm, based on the number of spins updated per second, is also lower for Wolff than for multi-spin coding implementations of Metropolis. Because of the Wolff algorithm’s smaller dynamic exponent, it is clear that it will become more efficient than Metropolis for sufficiently large lattices; however, “sufficiently large” might well be larger than the range of sizes of interest in a particular study. Work by Ito and Kohring estimates that Metropolis remains more efficient than Wolff, in terms of independent measurements per second, for system sizes as large as $`L=70`$ in two dimensions and $`L=100`$ in three dimensions (running on a scalar workstation). This is, of course, strongly dependent on the type of computer and the particular implementation of the algorithms used. For example, with the programs, algorithms and computers used in this study, we estimate that Wolff becomes more efficient than Metropolis for $`L32`$ for $`d=2`$ and $`L16`$ for $`d=3`$.

Another concern with the Wolff algorithm is its sensitivity to flaws in the random number generator used in the simulation. Small but significant systematic deviations from the exactly-known answer in the $`d=2`$ Ising model and other systems have been reported and investigated- using a variety of popular random number generators.- While the results of any simulation method can be biased by subtle correlations in the random numbers, the Wolff algorithm was found to be particularly susceptible. Despite these concerns about the Wolff algorithm, the dramatic reduction in the correlation time is a very tantalizing effect. If the speed and efficiency at equilibrating small-scale structures of the Metropolis algorithm is combined with the strength in decorrelating large-scale structures of the Wolff algorithm, the resulting hybrid algorithm could, in fact, be more efficient than either Metropolis or Wolff individually.

To test this possibility, we implemented a scalar hybrid algorithm which combines the Metropolis and Wolff algorithms. Simulations were performed on IBM RISC/6000, DEC Alpha, PC Linux and SGI Power Challenge workstations. The spins were stored as bit variables, with up to 32 spin variables packed into a single computer word. Note that the Metropolis algorithm can take advantage of this packing arrangement by effectively updating many spins in parallel using multi-spin coding techniques. This will result in substantial improvement in performance for increasing system size until all 32 bits are filled (for $`L64`$ in this implementation). While the Wolff algorithm cannot make full use of the multi-spin coding, it does benefit from the smaller memory requirements of the packed-spin representation. (Smaller memory means that more of the system can be stored in the computer’s cache memory which results in much better performance.)

The random number generators used for the simulation must be chosen with great care, especially for the Wolff algorithm. After performing extensive tests of several generators, we selected the following as being the fastest random number generators that would give us the correct answer within the precision of our testing. For the Wolff algorithm, we used a combination generator by l’Ecuyer that was recommended as a “perfect” random number generator in the Numerical Recipes column in Computers in Physics. With this program, we can produce a random number in 840.2 nsec on an SGI Power Challenge workstation with a 194 MHz R10000 processor. For the Metropolis part of the simulation, we used a faster, shift-register generator, R1279, which can produce a random number in 21.4 nsec.

To see how the hybrid algorithm behaves when poor random numbers are used, we ran a series of simulations deliberately using a bad random number generator for the Wolff algorithm. We thus used the R250 shift-register generator which is known to introduce significant systematic errors for the d=2 Ising model. We performed hybrid updates consisting of one Wolff update followed by 0, 1, 2, 3 and 4 Metropolis updates. (A simulation consisting of only Metropolis was also performed for completeness.) For each hybrid, 16 independent simulations consisting of $`3\times 10^6`$ hybrid steps were performed. The results for $`L=16`$ are shown in Fig. 1 for the energy and specific heat. For the internal energy the Wolff algorithm yields the wrong answer by an amount which is more than 35 times the calculated error bar. With the inclusion of $`50\%`$ of Metropolis flips this error is reduced by a factor of 10, and with $`80\%`$ Metropolis flips, no discernible error is seen. Very much the same behavior is seen in the specific heat, although the rate of convergence to the correct answer is slightly different.

Not only are the results like to be more correct if the different flipping mechanisms are mixed, but the performance is also improved. For the magnetization the relative efficiency of the hybrid algorithm, with $`5080\%`$ Metropolis flips added, is about $`30\%`$ greater than for Wolff alone, as can be seen from Fig. 2. It is surprising that even for $`L=64`$, where the Metropolis algorithm is much less efficient than the Wolff algorithm, the hybrid is significantly more efficient. Although for pure Metropolis, the relative performance becomes markedly worse as the lattice size increases, the same is not true for the hybrid algorithm. For the internal energy the relative efficiency, also shown in Fig. 2, is much better still for the hybrid algorithm, by more than a factor of two.

## 3 Hybrid algorithm: spin-3/2 Ising model

For models with higher values of spin, not only are Monte Carlo simulations, as well as specific algorithms, less ubiquitous than for their two-state counterpart, but no exact solution is still available for their critical temperatures. Thus, the basic ideas of the last section need to be extended to more general models, e.g. the spin-$`3/2`$ Ising model, where each spin state can assume values $`\pm 3/2,\pm 1/2`$. Although some spin-$`1`$ - and spin-$`3/2`$ - models have already been studied through Metropolis technique and cluster spin-flipping , there is still a lack of a detailed analysis of the statistical and systematic errors even in the simple Ising limit. So, before starting to implement a hybrid algorithm to this model it is interesting to see first what one gets with single-spin flipping procedures.

To analyse the statistical errors of some observable thermodynamic quantity $`A`$ we first applied just the Metropolis algorithm to the spin-$`3/2`$ Ising model. We ran $`6.02\times 10^6`$ Monte Carlo steps (MCS) per spin with $`2\times 10^4`$ configurations discarded for thermalization on different lattice sizes $`L`$ ($`8L128`$) and using the “perfect” random-number generator. We measured the energy $`E`$, magnetization $`M`$, fourth-order cumulant of the magnetization $`U`$ and the quadrupole moment $`Q`$ (the mean value of the square of the spins). Typical results of the relative error $`\mathrm{\Delta }A/<A>`$ for different lattice sizes $`L`$ are shown in Fig. 3 at $`t=k_BT/J=3.29`$, a value close to the critical temperature. Using ’coarse-graining’  we have estimated $`\mathrm{\Delta }A`$ through

$$(\mathrm{\Delta }A)^2=(<A^2><A>^2)/N,$$

for large enough $`N`$, where we divided our data into $`N=MCS/n`$ bins of different lengths $`n`$ ($`n`$ ranging from $`5`$ to $`10^5`$). The relative error in the magnetization and its cumulant increases as $`L`$ increases while for the energy and the quadrupole moment it stays almost constant. In terms of different degrees of self-averaging, $`M`$ and $`U`$ are non-self-averaging while $`E`$ and $`Q`$ exhibit a lack of self-averaging (the number of ‘effectively’ independent measurements through the computation of the correlation time $`\tau `$ is certainly necessary  for a more detailed analysis of the errors. This is, however, outside of the scope of the present work). We have also noted no significant changes in the errors by using different random number generators, even taking the poorer congruential one, and the data are also depicted Fig. 3. Fig. 4 shows the results of the magnetization cumulant $`U`$ as a function of temperature for different MCS for the lattice size $`L=128`$. Here, we have used the histogram technique at $`t_o=3.29`$ in order to obtain estimates for other values of temperatures close to $`t_o`$. It is worthwhile to analyse such behavior since we will use the crossings of the cumulants $`U`$ to locate the critical temperature of the present model. Besides having large error bars one can see that the mean value of the cumulant is strongly dependent on the number of MCS used to obtain the statistics. The dependence with the number of states generated for this lattice size is more pronounced by using the ’perfect’ random generator (note that the mean values of $`U`$ with $`6\times 10^6`$ MCS and ’perfect’ generator are comparable to those with $`3\times 10^6`$ MCS and the congruential one) although both converge to the same limit as the number of MCS gets very large. Within the error bars we also notice almost no systematic error due to the use of different random number generators, in contrast to the case of the Wolff algorithm which, with a bad random number generator for the spin-1/2 model, gives wrong results for the energy and specific heat (see Fig. 1). The same qualitative behavior of Fig. 4 (large error bars and a strong dependence of the cumulant with the number of MCS) is also seen for other lattice sizes $`L`$. Even by substantially increasing the MCS per spin one still gets large errors, mainly for the magnetization and its cumulant (see also Fig. 3).

We show in Fig. 5 the reduced pseudo-critical temperature $`t_c`$ as a function of $`L^1`$ (in fact $`L^{1/\nu }`$, where $`\nu =1`$ for the two-dimensional Ising universality class) obtained from the crossings of the fourth-order cumulant of the magnetization for different values of $`L`$ using just the Metropolis algorithm. Each point in that figure represents the crossing point of the cumulant $`U_L`$ of the lattice size $`16L128`$ with the corresponding cumulant of the smallest lattice $`U_8`$. Only a poor estimate of the critical reduced temperature can be achieved in this case which can be ascribed to the large error bars obtained in computing $`U_L`$ as well as its strong dependence on the number of MCS taken in the statistics. In particular, we have $`t_c=3.288(1)`$ with ’perfect’ and $`t_c=3.287(1)`$ with congruential random number generators which are, even so, comparable to the more recent series expansion result $`t_c=3.2878(22)`$. We can see that, in general, no systematic error due to random number generator is observed for the Metropolis algorithm. Moreover, within the error bars, very similar results are also obtained by running the symmetric heat bath single spin-flip procedure.

It is clear that one way to improve the accuracy of the location of the critical temperature with Metropolis can be done by increasing the MCS in order to achieve better statistics. This will require, of course, much more computer time, mainly for large lattice sizes. We can, however, use the results of the previous section in order to construct a hybrid algorithm where, with not much extra computer time, more precise results could be obtained. The first step is thus to implement a Wolff algorithm for this model. In a straightforward way, this implementation can be done by activating bonds between parallel nearest-neighbors spins $`S_i`$ and $`S_j`$ according to the probability $`p(K_{ij})=1\mathrm{exp}(2K_{ij}S_iS_j)`$ and, when the full cluster has been activated, all its spins are reversed. This procedure has, however, two main differences regarding the spin-1/2 systems which we have to keep in mind: i) now, the probability $`p(K_{ij})`$ depends on the particular configuration of the parallel spins and can have three possibilities, depending on whether $`\{S_iS_j\}`$ are $`\{\frac{3}{2}\frac{3}{2}\}`$, $`\{\frac{3}{2}\frac{1}{2}\}`$ and $`\{\frac{1}{2}\frac{1}{2}\}`$ (and also the corresponding reversed configurations); ii) this procedure alone is not ergodic in the sense that it does not take into account transitions between spin states with different spin magnitudes (it keeps fixed the number of $`\pm \frac{3}{2}`$ and $`\pm \frac{1}{2}`$ spins in each configuration and, for instance, the quadrupole $`Q=_iS_i^2`$ is always a constant). While the former is just a generalization of the bond probabilities activation for systems with more degrees of freedom, the latter is really a problem since we can not generate all possible configurations for the model. A mixed cluster algorithm has already been proposed to overcome such a non-ergodicity in the case of the spin-1 Blume-Emery-Griffths model. However, a natural hybridization procedure, based on the discussion of the last section, and also from embedding algorithms  proposed to study of spin-1 models  can be worked out here by simply alternating one Metropolis sweep with $`p`$ Wolff steps where $`p`$, in principle, can depend on the system size. The inclusion of alternate single spin-flip sweeps will make this hybrid algorithm ergodic and much simpler than a possible generalization of the mixed cluster procedure to the present spin-$`3/2`$ model. In order to test the efficiency of this hybrid algorithm we have done extensive simulations for the $`L=4`$ lattice where we can compare the results of the simulations with the exact ones. We ran a total of $`1.2\times 10^7`$ hybrid MCS per spins each one including $`p`$ Wolffs intercalated by one Metropolis sweep. The results are shown in Fig. 6 for the ’perfect’ random generator. One can readily see that all the results are in general compatible to the exact ones within the error bars. However, by including some Wolff steps the mean values initially oscillate for small $`p`$, have a better agreement for $`p5`$ and finally deviate for large $`p`$. We also note that the errors are almost the same for $`p=0`$ and $`p=9`$ and are slightly smaller around $`p5`$. The slight deviation of the error as a function of $`p`$ reflects the fact that we have a reasonable number of MCS per spin to obtain a good statistics for this small lattice, even with just the Metropolis algorithm (this is not the case for larger values of $`L`$, as we shall see below). Moreover, the worse results for large values of $`p`$ can be ascribed to the non-ergodicity of the present simple Wolff algorithm (nothing is gained by increasing the number of Wolff steps since we get stuck in the configurations having the same number of $`\pm \frac{3}{2}`$ and $`\pm \frac{1}{2}`$ spins). The overall picture then suggests the use of this hybrid algorithm with $`p5`$ (although $`p`$ can also vary with $`L`$). In order to test this assumption we applied this procedure to the $`L=128`$ lattice (and close to the critical temperature) and obtained the magnetization cumulant $`U`$ with $`5`$ Wolff steps. The corresponding results are also shown in Fig. 4. There is, in this case, no sensitive difference in the data by taking $`3\times 10^6`$ or $`6\times 10^6`$ MCS with the hybrid algorithm. Surprisingly, the statistical errors are now almost two orders of magnitude smaller than those with just the Metropolis algorithm (the errors in Fig. 4 for $`L=128`$ are in fact much smaller than the symbol sizes). To get this same precision with only Metropolis one would have to compute an order of $`10^8`$ configurations for this lattice size. The relative error of the energy, magnetization and its cumulant, and the quadrupole for other values of $`L`$ are shown in Fig. 7. While now they are almost constant for $`M`$ and $`U`$ (exhibiting lack of self-averaging), they decrease for $`E`$ and $`Q`$ (behaving now as self-averaging quantities). It is also important to notice that the hybridization process is again almost insensitive to the quality of the random number generator. The above results strongly indicate $`p=5`$ as a good trial also for other lattice sizes.

In Fig. 8 we present the reduced pseudo-critical temperature $`t_c`$ as a function of $`L^1`$ obtained from the hybrid algorithm described above with the ’perfect’ random generator through the crossings of the fourth-order cumulant of the magnetization $`U_L`$. It was possible, in this case, to get a good resolution from crossings of $`24L128`$ with three different smaller lattices, namely $`U_8`$, $`U_{12}`$ and $`U_{16}`$. The quality of the results are now apparent and yields the extrapolated value $`t_c=3.28799(7)`$. Just for completeness, in Fig. 8 we also give the corresponding values by taking the congruential generator with the present hybrid algorithm. As it is faster, we were able to use, with the same computing time, lattices as large as $`L=192`$ to get $`t_c=3.28789(7)`$. We have then, so far, the best estimate for the critical temperature of the two-dimensional spin-3/2 Ising model: $`t_c=3.28794(7)`$.

It is worthwhile now to address some comments regarding the universality of these models. Regardless the number of states each spin can assume, all $`d`$-dimensional systems are in the same (Ising) universality class. This fact is apparent in Figs. 5 and 8, (mainly the latter one) where the temperatures are all along a straight line as a function of $`L^{1/\nu }`$ with $`\nu =1`$ in two dimensions. However, two more universal quantities can be readily observed from the present simulations. First, the magnetization fourth-order cumulant $`U^{}`$ at the transition temperature can also be estimated from our data to give $`U^{}0.612(1)`$, a value expected for $`d=2`$ Ising systems undergoing a second-order phase transition. This result comes from Fig. 9 where each point was obtained by fixing $`t`$ at our estimate $`t_c`$ and looking the cumulant there for different lattices. Second, a quantity which is studied less often, is the probability distribution of the magnetization $`M_L`$, $`P^{}(b_oL^{\beta /\nu }M_L)`$ which, for large enough systems at the critical temperature is a universal function. In this equation $`b_o`$ is a non-universal constant chosen to give a unity variance for the distribution $`P^{}`$. Fig. 10 shows the fixed-point order parameter distribution for the two-dimensional Ising universality class obtained from models with spin-1/2,1,3/2 at the critical temperature and lattice size $`L=32`$. For each model $`10^7`$ steps were performed with Metropolis algorithm and using the R1279 random number generator. The quality of this match clearly reveals the hallmark of the $`P^{}`$ distribution for the Ising universality class.

## 4 Discussion

The results shown in the previous sections supply strong support for the use of hybrid algorithms as a means of effectively speeding up simulations and also improving the quality of the results. Another advantage is that efficient, parallel implementation of the hybrid algorithm on distributed memory machines is straightforward. A Wolff process running on one processor can “feed” states to other processors which then perform multiple Metropolis updates. The number of Metropolis updates can be varied to maximize load balancing. Data are gathered together from all states which have been generated and then used to construct histograms. This procedure can be enhanced still further by the inclusion of microcanonical updates which require no random numbers! One hybrid update would then consist of, e.g. 1 Wolff update plus 5 Metropolis updates plus 10 Microcanonical updates.

Although we have described hybrid algorithms for one of the simplest models in statistical mechanics (Ising), we believe that the lessons drawn from these studies will be more broadly applicable. For example, continuous spin systems may be (randomly) projected onto Ising models which can be easily simulated using these hybrid algorithms. Histogram analysis of the data can also be used in a similar fashion to produce extremely high resolution results. Of course, the relative performance of each component of the hybrid algorithm will depend upon the specific model, so that “tuning” will be required for each study. Furthermore, these single ensemble hybrid methods can be combined with other algorithms to further improve performance. For the Ising model the microcanonical method is extremely fast and can be easily included. For a classical Heisenberg model, the over-relaxation method provides an effective microcanonical simulation component to a hybrid algorithm.

In summary, we have demonstrated that hybrid Monte Carlo spin-flip algorithms, which include “slower” Metropolis steps, can be made to be effectively faster than cluster flipping algorithms. Furthermore, and perhaps more significantly, they yield substantially more accurate results than does the simple Wolff algorithm (for the spin-1/2 model) or single Metropolis algorithm because the alternation of updating methods breaks up random number correlations.

acknowledgments

This research was supported in part by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq-Brazilian Agency) and NSF grant DMR-0094422.


text:
# Complex Systems: a Physicist’s Viewpoint

## 1 Introduction

In recent years physicists have been deeply interested in studying the behavior of complex systems. The result of this effort has been a conceptual revolution, a paradigmatic shift that has far reaching consequences for the very definition of physics.

In the nutshell, physics is an experimental science in which theoretical predictions are compared to experiments. It differs from other sciences of nature in the crucial role of mathematics: physicists describe the phenomena by using a mathematical model and their predictions are obtained by mathematical reasoning.

In this definition of physics the word prediction plays a crucial role. Naively one could think that the meaning of the word prediction is quite evident and no further explanation is needed. In reality its meaning has changed already in the past and it is still changing now.

Let me recall using modern language what a prediction was at the time physics was born. We assume that the state of the system at the time $`t`$ is represented by a a vector $`X(t)`$ with $`N`$ components.

In the case of classical mechanics the state of a point-like particle is described by the three coordinates and by the three components of the velocity. The corresponding vector $`X(t)`$ has six components. An experiment consists in determining or imposing the position and the velocity of the particle at a given time and measuring it a later time. It is crucial that (at least in principle) the experiment can be done many times by imposing the same initial conditions. If the difference between the initial and the final time is the same, the system is always found in the same final state. When repeated, the experiment always gives the same result. The theoretical task is to predict this unique, reproducible result.

This can be done by using the equations of motion which gives the time evolution of the vector $`X(t)`$. In the case of classical mechanics, they have the form of Newton’s second law

$$\frac{dp}{dt}=F(x),\frac{dx}{dt}=\frac{p}{m}$$

where $`F(x)`$ is the force. Once that the form of the force has been established it is a well posed mathematical problem to compute the final conditions as function of the initial one.

The same perspective applies not only to simple objects but also to much more complicated situations. The description of the system may need more than one position and one velocity. The equations of motion, which are in general of the form

$$\frac{dX}{dt}=G(X(t)),$$

in principle allow the prediction of the final state for each initial state. The theory (which consists in assuming a particular form of the function $`G(X)`$) can be tested by measuring $`X`$ at the initial time, predicting $`X`$ at the final time and verifying the prediction.

In reality, the mathematical difficulties may be substantial even if the system is apparently simple. To predict the evolution of a three body system (like Sun, Earth and Moon) is a rather difficult task and it can be done only with a lot a work. At the beginning of the nineteenth century Laplace was writing that an infinitely intelligent mathematician would be able to predict with certitude the future by observing the state of the universe now and using his knowledge of the laws of motion.

In this classic, deterministic framework the word prediction has a clear meaning. Unfortunately we are not infinitely intelligent mathematicians and besides an infinitely able experimentalist is needed to measure the state of the system with extremely high accuracy. The Newton-Laplace point of view could be applied only to a restricted class of phenomena. In order to consider other phenomena which could not be studied from this point of view, it was necessary to change the general philosophy, by introducing probabilistic concepts and probabilistic predictions.

There have been three revolutions in physics which all went in this direction  and as a consequence they have changed the meaning of the word prediction. They are:

* (1) The introduction of statistical mechanics and of the first probabilistic reasoning by Maxwell, Boltzmann and Gibbs in the second half of the last century.
* (2) The discovery of quantum mechanics at the beginning of this century.
* (3) The study of complex systems and the related techniques that have been developed in these last years.

As an effect of these revolutions, the word prediction acquired a weaker meaning. Predictions in the context of the new paradigm are not acceptable with the old one (and sometimes the supporters of the old point of view try to deny to them a scientific validity). The positive consequence of the process is that the scope of physics becomes much larger and the constructions of physics find many more applications.

My aim is to present point (3). However, it is better to start by discussing point (1) in order to understand in which way the study of complex systems forces us to use the word probability in a wider context. For our purpose we can neglect point (2), because extends the concept of probability in a different direction 111In quantum mechanics we have to deal with the fact that also in principle the results of a single experiment are not reproducible..

## 2 Introducing Probability in Physics

At the time of Boltzmann  and Gibbs, the main motivation for leaving the Newton-Laplace point of view was not that it was wrong, but that it was useless in many cases. Let us consider the following experiments which apparently can be repeated many times with the same result:

* We bring water to 100 centigrade: it boils.
* We bring water to 0 centigrade: it freezes.

When we study this phenomenon from a microscopic point of view, we face a serious problem. The experiment should be done by measuring the positions and the velocities of all the billions of billions of atoms and we should later use this information to compute the trajectories of all the atoms. This difficulty was bypassed by discovering that for a system composed by many many particles, practically all the initial conditions (at fixed total energy) lead to the same macroscopic behavior of the system. The task of measuring everything is not only impossible, it is also useless. If we neglect the possibility of some very special initial configuration - which happens with extremely low probability \- the system always behaves in the same way.

The final predictions are the following:

* Water practically always boils at 100 degrees. The probability is so high that it is extremely likely that this kind of predictions has never failed in the whole history of the universe.
* If we measure the velocity of a single molecule of water we cannot predict anything precise. However we can compute the probability that the molecule has a velocity $`v`$ ($`P(v)\mathrm{exp}(Av^2)`$ where $`A=m/(2kT)`$). Therefore if we measure the velocity of a molecule $`M`$ times the observed distribution of velocity should approach the theoretical curve when $`M`$ goes to infinity.

In general, we could say that for these large systems there are two kinds of quantities:

* Some quantities can be predicted with certitude and have always the same value inside a small interval (usually of order $`N^{1/2}`$ for a system with $`N`$ degrees of freedom).
* Other quantities do not take always the same value. In this case a probabilistic prediction is possible.

## 3 Deterministic Chaos

In the previous discussion, the main motivation for introducing probability was the large number of particles present. Only in the second half of this century was it realized that there are many systems with a small number of particles, or more generally with a small number of degrees of freedom, for which it is also necessary to use probabilistic arguments .

In these systems (which sometimes are called chaotic) the trajectory cannot be predicted for large times, because it is too sensitive to the initial conditions. The difference in the position of the trajectories of two systems which have a very similar initial conditions increases exponentially with time. In other words, even a very small incertitude on the initial condition leads to a total loss of knowledge after a characteristic time $`\tau `$. More precisely if

$$X_1(0)X_0(0)ϵ$$

the difference $`X_1(t)X_0(t)`$ grows as function of the time $`t`$ as $`ϵ\mathrm{exp}(t/\tau )`$ until it reaches a value of order 1.

Deterministic predictions are possible here on a time scale smaller than $`\tau `$, but they become impossible at later times because we cannot measure the initial conditions with infinite precision.

A very simple example is provided by two (or more) balls on a billiard table without friction. For generic choice of the initial condition the two balls will collide from time to time, and after each collision it becomes more difficult to predict the position of the balls.

Fortunately, in many cases the probability distribution at large times for finding the system in a given configuration is independent on the initial conditions and it can be predicted. So here also we can compute only the probability distribution of some variables, not the exact evolution.

The new framework given by statistical mechanics should now be clear. We have a system, we know the equation of motion. In principle, an exact knowledge of the initial conditions allow us to compute the exact evolution of the systems. In reality this task may be impossible because we do not know the initial conditions with sufficiently high accuracy. It is also possible that the computation is terribly complicated. We can however make progress and get much more insight if we restrict ourselves to the task of predicting the probability distribution of the system at large times (in short its behavior). Given the equation of motion $`G`$ (see eq. 1) we define the probability $`P_G(X)`$  as

$$𝑑Xf(X)P_G(X)=\underset{T\mathrm{}}{lim}\frac{_0^Tf(X(t))𝑑t}{T},$$

where $`f(x)`$ is an arbitrary test function. For ergodic systems the probability $`P_G(X)`$ does not depend on the choice of the initial point for almost all the choice.

The task of statistical mechanics consists in the computation for given $`G`$ of the probability $`P_G(X)`$ and of its properties.

## 4 Complex Systems

It was recently realized that this approach described in the previous section (i.e. to predict the behavior of the system from the knowledge of the equations of motion) fails in the case of complex systems for reasons that are very similar (albeit in a different contest) to those which led to the abandonment of the Newton-Laplace view-point.

There are many possible definitions of a complex system. I will use the following one. A system is complex if its behavior crucially depends on the details of the system. This dependence is often very difficult to understand 222 A single complex system may also display different types of behavior and a small perturbation is enough for switch from one behavior to another. For example an animal may: sleep, dream, run, hunt, eat, play…. In other words, the behavior of the system may be extremely sensitive to the form of the equations of motion and a small variation in the equations of motion leads to a large variation in the behavior of the system. More precisely for two systems with $`N`$ degrees of freedom with equations of motion $`G_0`$ and $`G_1`$, where

$$G_1(X)G_0(X)=O(ϵ),$$

we could have that for small $`ϵ`$ but large $`ϵN`$,

$$P_{G_1}(X)P_{G_0}(X)=O(1).$$

In this happens it is practically impossible to compute the asymtotic probabily as function of $`G`$.

My aim is to show that although a system is complex, it is still possible to get predictions for its behavior, but these predictions will be now of a probabilistic nature. I will illustrate this point in a particular example. The prototype of a complex system in physics is spin glass . I will choose here something more familiar to most people, a large protein.

A protein is composed by a long chain of a few hundred amino-acids and its chemical composition is specified by the sequence of aminoacids (primary structure). In physiological situations, normally a protein folds in an unique way (tertiary structure). During the folding, the protein goes to the configuration of minimal energy (more precisely, to the configuration of minimal free energy). However, there are many quite different configurations of the protein which have an energy near the minimal one. There are proteins, which are called allosteric, which have two configurations which practically the same energy). Sometimes, if we perturb the protein by a small amount, e.g. changing the pH or changing only one aminoacid in the chain, the folding changes dramatically. We can also select one of the two different configuration by the binding of a ligand.

The variations in the shape of mioglobulin in different conditions are crucial for its respiratory functions. More generally, the possibility of switching from one configuration to another is one of the fundamental mechanisms which allows proteins to work as enzymes. For example the change of configuration of allosteric proteins it is used for doing useful work, e.g. for contracting muscular fibers.

The existence of allosteric proteins implies that a small variation in the form of the potential among the components of the protein changes the folding of the protein dramatically. It is thus evident that a small error in the form of the interaction among the different components of the protein or of the interaction with the solvent would lead to quite wrong results.

Natural proteins are on the borderline of what can actually be understood; we have not yet computed their folding properties, but this does not seem completely impossible, especially for the smallest ones. However, to compute the folding properties of much large proteins seems to be a hopeless task.

If an approximate knowledge of the equations of motion does not allow us to compute the behavior of a complex system, we can give up this task and restrict ourselves to compute the probability of a behavior when the equation of motions are chosen randomly inside a given class.

In other words we can study the same problem in two different ways:

* In the old approach we consider a given protein (e.g. mio globulin) and we compute the form of the two foldings with minimal energy and their energy difference as a function of the chemical properties of the solvent.
* In the new approach we consider a class of proteins (to which mioglobullin belogns) and we compute the probability distribution of the energy difference among the two foldings with minimal energy .

Of course, this second computation will not tell us the properties of mioglobulin. It will give us predictions that cannot be tested on a single protein (in the same way the probability distribution of the velocities cannot be tested by measuring only one velocity); we need to repeat the experiments on many different proteins.

Taking the new point of view, we often gain insight. The behavior of a given complex system may be obtained only after a very long computation on a computer and often we cannot understand the deep reasons of the final result. The study of the probability distribution of the behavior can sometimes be done analytically and we can follow all the steps that lead to the conclusions.

The advantage of the second point of view is that its results may be easily generalized to other systems. For example long RNA molecules also fold in a characteristic way, in the same way that proteins do. A general explanation which shows that some characteristic of folding are shared by all polymers composed by different monomers may be applied to RNA as well as to protein

There are simple models (e.g. spin glasses ) in which this program works. The computation of the behavior of a given system is extremely difficult and sensitive to the minimal detail and can be done only with lengthy computer simulations. The analytic computation of the probability distribution of the behavior for a generic system can be done analytically and these probabilistic predictions have been successfully verified.

The extension of this program to other systems leads to two different types of difficulties.

* The generalization of the methods which allow the computation of the probability distribution of the behavior to models which are not as simple as spin glasses, is technically difficult and progress in this direction are rather slow.
* This point of view is quite different from the traditional one (mathematically it corresponds to the use of imprecise probabilities). The predictions are not made for the properties of a given system, but for the probability distribution of those properties which change with the system. We are not enough familiar with this new point of view to appreciate all its potentialities.

As usual a change in the paradigm leads to a change in the questions posed. We do not have to ask how a particular system behaves. We have to ask which are the general features of the behavior of a system belonging to a given class.

In spite of these difficulties, this new approach seem to be absolutely necessary to get some understanding in really complex systems. We have seen that the arguments are very similar to those described before. In the two cases they can be summarized as follows

* (1) The extreme sensitivity of the trajectory on the initial conditions forces us to study the probability distribution of the systems at large times. This must be done in spite of the fact that an exact knowledge of the initial condition determines the trajectory at all times. The behavior of the system can be computed from the equations of motion.
* (2) The extreme sensitivity of the behavior on the detailed form of the equation of motion forces us to study the probability distribution of the behavior of the system. This must be done in spite of the fact that an exact knowledge of the initial conditions determines the behavior. The probability distribution of the behavior of the system can be computed from the probability distribution of the equations of motion. In other terms, given a probability distribution $`\rho (G)`$ we can compute the probabillity $`𝒫[P_G(X)]`$, which of course depends on $`\rho `$.

The advantage of this new point of view is to expand the range of application of physical reasoning. We shall see now how it puts the relations between physics and biology in a new perspective.

## 5 Possible Biological Applications

In biology problem we face a very difficult of how to use the immense amount of information we gather at the molecular level in order to understand the behavior of the whole organism.

Let us consider an example. One of the simplest organisms is E. Coli . Its genome (which is known) codes for about 3000 different proteins which interacts among themselves. Some proteins promote the production of other proteins, while other proteins have a suppressor effect.

In principle we can gather the information on the properties of the proteins, on the interaction among the proteins (for simplicity let us neglect all the other chemical elements of the cell). At the end we can use this information to construct a model of the living cell. The model may lead to a system of few thousand (or more) coupled differential equations which can be studied by a lengthy computations on a computer.

If the information is accurate enough, the model will describe a living cell. However we know that in real life many mutations are lethal, and therefore it is quite likely that imprecision in the form of the interactions among proteins will lead to a non-living cell or to a living cells with a quite different behavior. It is also clear that such a gigantic computation (although welcome) will not capture the essence of life. Indeed it should be repeated nearly from scratch if instead of E. Coli we consider a different organism.

The previous discussion suggests that we could give up the aim of deriving the properties of a given organism from its chemical components, and ask different questions, e.g. what are the general properties of a living organism and how they do change from organism to organism.

An example of this last of question has been studied by Kaufmann . There is a dependence of the number different cell types on the number of genes in all kind of living organisms.

$$\mathrm{\#}celltypes(\mathrm{\#}genes)^{1/2}$$

As Kaufmann points out, this fact must have a general explanation and he proposes one based on mathematical considerations of evolution.

An other empirical law, which calls for a general explenation is

$$\mathrm{\#}speciesinaisland(\mathrm{\#}genes)^{1/4}$$

In these two cases we can hope to compute the exponents. However in many other cases this point of view is so different form the traditional one, that it is very hard to find the right questions and also to answer to them.

A field in which this approach seems to be mandatory is the study of the origin of life. Here it is crucial to understand which are the enzymatic capabilities of the prebiotic material (maybe long randomly assembled RNA chains) before evolution starts. We would like to know which of the properties of these long chains had to be selected and which ones where already present before selection.

This approach is quite different from a redutionistic one, as far as it puts the stress on the behavior of the whole system. For the moment the largest impact of these ideas in biology has been in the field of neural networks, but in the other fields the progress is quite slow. We need to sharpen our theoretical tools in order to be able to predict the typical behavior for different classes of systems. We are seeing now only the beginning of these investigations in mathematical physics.

When the physical instruments become more robust and our theoretical command increases, the interaction with biology will become more easy. I am convinced that in the next century a much more deep understanding of life will come from this approach.