text:
# On the Critical Mass: the case of white dwarfs

## 1 Introduction

The problem of critical mass against gravitational collapse is one of the most important in relativistic astrophysics and has been the subject of a very large number of theoretical investigations with many observational consequences, the most important one being the paradigm for the identification of black holes in binary X-ray sources. The concept of critical mass appears in a large variety of astrophysical settings: white dwarfs, neutron stars, galactic nuclei. The occurrence of gravitational collapse in galactic nuclei is still matter of great debate both from the observational and theoretical point of view: the physical conditions leading to the onset of gravitational collapse, even the matter constituents of such a process, are still far from settled although the existence of extremely massive relativistic systems in the range of $`10^610^8M_{}`$ in galactic nuclei is certain. The occurrence of gravitational collapse in neutron stars has obtained grandiose observational verification in binary X-ray sources although, from a theoretical point, it still presents problems of unsurpassed difficulty due to theoretical and observational unknowns in the behavior of matter compressed at supra-nuclear densities. In this case only a firm upper limit on the numerical value of the critical mass has been obtained from first principles: the validity of general relativity, the causality non-violation principle and the behavior of matter at a fiducial nuclear density.

In this respect the case of white dwarfs appears most interesting: it is amenable to a complete theoretical treatment which is also testable by refined observational measurements thanks to the superb observations of binary pulsars. The critical mass for white dwarfs occurs at densities relatively low in comparison with nuclear densities and many of the intricacies of the behavior of matter in bulk at supra-nuclear densities can be completely neglected. On the other hand the theoretical analysis of the white dwarf critical mass can become the arena for a detailed treatment of gravitational and electromagnetic effects at relativistic regimes in the fully tested region of validity of these field theories: general relativity, quantum theory and electromagnetism. I will devote this contribution to this topic which has reached complete maturity and in principle can give important observational clues for the understanding of astrophysical high-precision measurements.

## 2 Historical Background

That white dwarfs cannot have an arbitrary large mass and that a critical mass against gravitational collapse must exist in their configuration of equilibrium reachable by increasing their central densities was clearly understood in the thirties independently by Subranaman Chandrasekhar and by Lev Davidovich Landau. Their treatments differ enormously in style.

After publishing his result in an article, Chandrasekhar gave his derivation in a 343 page classic book which summarized: a) some of the results originally presented in German in the classic textbook Gaskugeln by Emden, b) basic aspects of quantum mechanics, developed in the classic book of R.H. Fowler and in the work of P.A.M. Dirac, and 3) the detailed computations of the mass-radius relation for white dwarfs and their agreement with observations with the clear statement of the existence of a critical mass equal to $`1.45M_{}`$.

The work of Lev Landau is totally different in its essential approach: Landau also uses all the expertise which could be gained from the classic work of Emden and quantum theory, at that time just formulated, and the corresponding relativistic corrections. Similarly he reaches the conclusion of the existence of a critical mass for a white dwarf, but now, everything, in a dramatic paper of no more then three pages. There by the stroke of crucial shortcuts in the physics of self-gravitating degenerate matter in relativistic regimes, Landau reaches that unexpected conclusion. Landau’s paper is not only beautiful for the masterful interplay of physics and astronomy, but can also be considered the first paper introducing a new style of work in the scientific analysis of astronomical phenomena.

Another non-negligible difference must be mentioned distinguishing the work of Landau and Chandrasekhar: Landau did not believe the possible existence in nature of the new concept of a critical mass he had just obtained, while Chandrasekhar never revealed such a critical attitude towards it.

At the time of presentation of his work in Cambridge, Chandra met with very strong resistance and criticism from Arthur Eddington, so much so that his paper published in the proceedings of the Royal Astronomical Society was preceded by an article of strong criticism by Eddington himself which clearly stated, with reference to Chandra’s work, “… but its physical foundation does not inspire confidence, since it is a combination of relativistic mechanics with non-relativistic quantum theory.”

I was always interested in understanding some of the reasons which motivated Eddington in his criticism. After all, at the time, Eddington was one of the greatest experts in the theory of general relativity and was the author of fundamental textbooks both in astronomy and general relativity: for me more than his expressed criticism toward Chandrasekhar’s work, it was interesting to understand whether there was any fundamental theoretical issue adopted by Chandrasekhar which could be disturbing to the eyes of a relativist.

The answer to Eddington’s criticism came finally in 1974 from Chandrasekhar himself in his outstanding article “Why are the stars as they are” presented at the Varenna Summer school on the Physics and Astrophysics of Neutron Stars and Black Holes.

Nevertheless reconsidering recently the issue of the white dwarf critical mass, I have found a basic unproven assumption in Chandrasekhar’s work. The distribution of the degenerate electrons is assumed to be uniform in the stars and unaffected by Coulomb interaction with the nuclei all the way to infinite densities. Indeed in Chandrasekhar’s work the value of the critical mass is reached for an infinite central density of the star. On this issue Landau had made two important remarks:

1. In 1938, following the earlier suggestions by Hund, Landau pointed out “it is well known that matter consists of nuclei and electrons. Nevertheless it can be shown that in bodies of very large mass, this usual “electronic” state of matter can become unstable. The reason for this lies in the fact that the “electronic” state of matter does not lead to extremely great densities, because at such densities electrons form a Fermi gas having an immense pressure. On the other hand, it is easy to see that matter can go into another state which is much more compressible—the state where all the nuclei and electrons have combined to form neutrons.”. Later Yakov Borisovich Zel’dovich pointed out, quantitatively, that this “neutron condensation” for inverse beta decay starts to occur for hydrogen already at $`10^7\frac{g}{cm^3}`$. Then Harrison et al pointed out that all configurations of white dwarfs after the onset of the inverse beta decay are indeed unstable and the critical mass is reached not at an infinite density but at precisely the finite density marked by the onset of inverse beta decay. This result was confirmed by the analysis of pulsational modes of White Dwarfs by Chandrasekhar  (see also Harrison et al).
2. The second contribution of Landau is contained in the treatise with Yevgeny Lifshitz, where Landau gives an estimate of the effects of the electromagnetic interaction between nuclei and electrons and explicitly shows how these interactions can be neglected only in the asymptotic regime for very large densities. The fiducial densities where indeed the Coulomb interactions can be neglected occur at $`\rho >10^7g/cm^3`$ for Z=56, namely at densities too high if compared to the density regime at which the stable white dwarfs and their critical mass occur. The interactions between nuclei and degenerate electrons therefore are not negligible for white dwarfs: they were certainly difficult to evaluate at the time of the Chandrasekahr’s original work.

Actually an interesting alternative approach, with respect to the evaluation of the interaction term between the electrons and the nuclei, was advanced two years earlier than the work of Chandrasekhar and Landau by the Soviet physicist Yacov Ilich Frenkel. This paper, quoted by Chandrasekhar in his classic book, has been in the past and is still today surprisingly ignored. Frenkel proposed to use a relativistic Thomas-Fermi model within a Wigner-Seitz approximation in order to describe stellar matter. The astrophysical motivations were not clear to Frenkel, though the theoretical formulation he proposed would have deserved a much more thorough examination.

An important turning point came in this field due to the fundamental work of Feynman, Metropolis and Teller. While at los Alamos, they considered the equation of state of compressed matter.

It was Ed Salpeter who first applied the Frenkel equations and the Feynman-Metropolis-Teller approach to the treatment of white dwarfs: Salpeter’s treatment was somewhat affected by the intricacies of the numerical integrations and the drastic approximations he adopted in the numerical solutions.

We ourselves rediscovered the Frenkel Equations and the relativistic Fermi-Thomas treatment applied to an appropriate Wigner-Seitz problem following the study of self-gravitating systems  and the riformulation of the treatment of these systems by using a relativistic generalization of gravitational Fermi-Thomas like equations   . Today with relatively simple numerical computations, we can reconsider the problem and present a comprehensive treatment by first considering the issue of the electromagnetic interactions between nuclei and electrons and turning later to the issue of inverse beta decays. The new treatment leads to the disappearance of all scaling laws with the chemical composition of the stars, which have been commonly assumed in all treatments up to now. A new value of the critical mass of white dwarfs is obtained, smaller than the one originally proposed by Chandrasekhar and by Salpeter, and one which strongly depends on the chemical composition of the star. These results appear in principle to be observationally testable. Particularly attractive is the possibility of explaining the observed masses of binary pulsars: these systems are clearly neutron stars which could be formed by the onset of instability of stellar cores very close to the critical mass of white dwarfs as presented by Chandrasekhar. The evaluation of their masses will be clearly affected by our computations.

## 3 The Relativistically Compressed Thomas-Fermi Atom

In the well known Feynman-Metropolis-Teller treatment of the classical Fermi-Thomas model an accurate numerical analysis is made of atomic configurations submitted to external pressure. Following Ferreirinho et al.  and Ruffini and Stella et al. we reconsider such a treatment in the framework of a relativistic generalization of the Fermi-Thomas model, which is here briefly derived and solved for different atomic species in different states of compression. Results so obtained are applied to determine the equation of state of cold degenerate stars (white dwarfs) to show how the relativistic and electromagnetic effects acting on microscopic scales affect the overall structure of the star. Finally, the equation of state obtained is used to determine the equilibrium configurations of white dwarfs. A comparison is made with the classic work of Chandrasekhar and Salpeter to show the differences introduced by the use of the new equation of state, and to show how those treatments correspond respectively to a zero-order and first-order approximation to the one we present here.

Let us consider the spherically symmetric problem of a nucleus with Z protons and A nucleons interacting with a fully degenerate gas of Z electrons. The equation of electrostatics for this problem reads, with usual meanings for the symbols appearing in it,

$$\mathrm{\Delta }V(r)=4\pi en_e(r)4\pi en_p(r),$$
(1)
where we have introduced the number density of electrons $`n_e(r)`$ and of protons $`n_p(r)`$.

In Eq. (1) the quantities $`V(r)`$ and $`n(r)`$ are of course not independent, being related by the equilibrium condition for a relativistic gas in a Coulomb external potential (see )

$$c\sqrt{p_F^2+m^2c^2}eV(r)=\text{const}E_F,$$
(2)
where $`E_F`$ stands for *Fermi energy* of the electrons.

To put the Eq. (1) in a simple and dimensionless form we introduce the new function $`\mathrm{\Phi }(r)`$, related to the Coulomb potential by

$$\mathrm{\Phi }(r)=V(r)+E_F/e$$
(3)
and the corresponding dimensionless function $`\chi (r)`$, implicitly defined by

$$\mathrm{\Phi }(r)=\frac{Ze\chi }{r}$$
(4)
Similarly we have the dimensionless expression for the radius $`r=bx`$, where

$$b=(3\pi )^{2/3}\frac{\mathrm{}^2}{me^2}\frac{1}{2^{7/3}}\frac{1}{Z^{1/3}}.$$
(5)
We then have

$$p_F=2mc\left(\frac{Z}{Z_{cr}}\right)^{2/3}\left(\frac{\chi }{x}\right)^{1/2}\left[1+\left(\frac{Z}{Z_{cr}}\right)^{4/3}\frac{\chi }{x}\right]^{1/2},$$
(6)
where

$$Z_{cr}=\left(\frac{3\pi }{4}\right)^{1/2}\left(\frac{\mathrm{}c}{e^2}\right)^{3/2}2462.4.$$
(7)
From

$$n_e=\frac{p_F^3}{3\pi ^2\mathrm{}^3}$$
(8)
we then obtain

$$n_e=\frac{Z}{4\pi b^3}\left(\frac{\chi }{x}\right)^{3/2}\left[1+\left(\frac{Z}{Z_{cr}}\right)^{4/3}\frac{\chi }{x}\right]^{3/2}.$$
(9)
Let us express the second term of the source of Eq. (1) in terms of dimensionless quantities: we assume a homogeneous nucleus, with a radius given by the approximate formula

$$r_{nuc}=1.2A^{1/3}\mathrm{Fermi}$$
(10)
The number density of protons is therefore

$$n_p=\frac{3Ze}{4\pi r_{nuc}^3}\mathrm{\Theta }\left(x_{nuc}x\right)$$
(11)
Finally we can write Eq. (1) in the form

$$\frac{d^2\chi }{dx^2}=\frac{\chi ^{3/2}}{x^{1/2}}\left[1+\left(\frac{Z}{Z_{cr}}\right)^{4/3}\frac{\chi }{x}\right]^{3/2}\frac{3x}{x_{nuc}^{}{}_{}{}^{3}}\mathrm{\Theta }\left(x_{nuc}x\right)$$
(12)
where $`x_{nuc}`$ is the dimensionless size of the nucleus ($`r_{nuc}=bx_{nuc}`$). Eq. (12) is what we call the “*generalised dimensionless Fermi-Thomas equation*.”

The first boundary condition follows from $`\chi r\mathrm{\Phi }`$ and therefore $`\chi \stackrel{r0}{}0`$, and so

$$\chi (0)=0.$$
(13)
The second condition comes from the normalization condition

$$N=_0^{r_0}4\pi n_er^2𝑑r=Z_0^{x_0}\frac{\chi ^{3/2}}{x^{1/2}}\left[1+\left(\frac{Z}{Z_{cr}}\right)^{4/3}\frac{\chi }{x}\right]^{3/2}x𝑑x$$
(14)
with atom size $`r_0=bx_0`$. Using Eq. (12) we obtain

$$N=Z_0^{x_{nuc}}x\chi ^{\prime \prime }𝑑x+\frac{3Z}{x_{nuc}^3}_0^{x_{nuc}}x^2𝑑x+Z_{x_{nuc}}^{x_0}x\chi ^{\prime \prime }𝑑x$$
(15)
which gives the simple relation

$$N=Z\left[x_0\chi ^{}(x_0)\chi (x_0)+1\right].$$
(16)
For a neutral atom we have $`N=Z`$ and the Eq. (16) reads simply

$$x_0\chi ^{}(x_0)=\chi (x_0)$$
(17)
Note that the physical quantities of interest, such as the Coulomb potential and the density of electrons do not show any singularity in the center, nor on the surface of the nucleus, since they only depend on the function $`\chi `$ and his first derivative. The only discontinuity appears in the second derivative of $`\chi `$ due to our rough assumption of a homogeneous nucleus.

## 4 The equation of state

Under these hypothesis, we determine the equation of state of compressed matter simply by computing for each value of the atomic compression parameter $`x_0`$ the corresponding values of pressure and mass density: the first is generated by the degenerate gas of electrons and the second is given by the mass density of nuclei. We assume that the nuclei are arranged in a Wigner-Seitz lattice, each cell being filled by a relativistic gas of degenerate electrons. The shape of the lattice cell is unimportant in the lowest approximation, so we will assume a spherical cell, with a radius equal to the atomic radius. The average density of cells is thus

$$\overline{\rho }(x_0)=\frac{Am_0}{V_{cell}}=\frac{3Am_0}{4\pi (bx_0)^3},$$
(18)
where $`A`$ and $`m_0`$ are respectively the number and mass of the nucleons.

In this scheme the pressure is also easily determined, being simply the pressure generated by a Fermi gas of number density equal to the number density of electrons at the surface of the atom $`n_e(x_0)`$, which can be expressed in terms of dimensionless quantities as follows

$$n_e(x_0)=\frac{Z}{4\pi b^3}\left(\frac{\chi (x_0)}{x_0}\right)^{3/2}\left[1+\left(\frac{Z}{Z_{cr}}\right)^{4/3}\frac{\chi (x_0)}{x_0}\right]^{3/2}.$$
(19)
Our treatment differs from the zero-order approximation of Chandrasekhar (uniformly distributed electrons) and also from the first-order approximation of Salpeter in which the Thomas-Fermi equation for each value of Z and $`x_0`$ was not solved and instead an expansion of the number density of electrons was introduced

$$n(x)=n_0+ϵ(x).$$
(20)
This expansion applies only when the electrostatic interaction within nuclei and electrons is small with respect to the Fermi energy of electrons, i.e. in the limit of very high densities.

## 5 Equilibrium Configurations for White Dwarfs

We can now compute the equilibrium configurations for white dwarfs. The results are shown in figures 1 and 2. We obtain, for the same central density, less massive configurations than those obtained by both Chandrasekhar and Salpeter. Namely for a fixed mass density, i.e. for a fixed value of A and of the compression parameter, our pressure is systematically lower than the one generated by a uniform shell of electrons.

Numerical integration of the models is performed here for stars of oxygen and iron in the framework of a Newtonian theory of gravitation, as well as in General Relativity. We show that the first one is actually a fairly good approximation up to central densities of the order of $`10^9g/cm^3`$.

It is shown that deviations from classical models are very important at the lower densities, up to 30% from the Chandrasekhar model and 10% from the Salpeter one.

We stop the integration at the critical density corresponding to the onset of inverse beta decays and we thus find the maximum mass for equilibrium configurations. This critical configuration marks also the onset of gravitational instability (see Harrison et al. ).


format: text

text:
# Sudden Hadronization in Relativistic Nuclear Collisions

\[

## Abstract

We formulate and study a mechanical instability criterion for sudden hadronization of dense matter fireballs formed in 158$`A`$ GeV Pb–Pb collisions. Considering properties of quark-gluon matter and hadron gas we obtain the phase boundary between these two phases and demonstrate that the required deep QGP supercooling prior to sudden hadronization has occurred.

PACS: 12.38.Mh, 25.75.-q, 12.40.Ee,

\]

Hot and dense hadron matter fireball is formed in central collisions of relativistic 158$`A`$GeV heavy Pb-nuclei with Pb-target, comprising a new state of matter . Driven by internal pressure, a fireball expands and ultimately a breakup (hadronization) into final state particles occurs. Early on in the hadron production data analysis it was discovered that strange hadrons emerge from a source in which strange $`s`$ and antistrange $`\overline{s}`$ quarks have same size phase space . This is the case if final state hadrons are directly produced by the deconfined quark matter phase.

The required sudden fireball breakup could arise if a fireball made of a the new form of matter significantly supercools, and in this state encounters a strong mechanical instability . Despite extensive ensuing study, the mechanisms determining when and how sudden hadronic particle abundances are formed (chemical freeze-out) have not been fully understood . However, there is growing evidence that the fireball breakup occurs over a relatively short period of time .

We propose and study here a natural mechanical instability criterion ensuing the fireball expansion into the metastable supercooled state. We consider the exploding fireball dynamics in its center of momentum frame of reference. The surface normal vector of exploding fireball is $`\stackrel{}{n}`$, and the local velocity of matter flow $`\stackrel{}{v}_\text{c}`$. The rate of momentum flow vector $`\stackrel{}{𝒫}`$ at the surface is obtained from the energy-stress tensor $`T_{kl}`$ :

$$\stackrel{}{𝒫}=P^{(i)}\stackrel{}{n}+(P^{(i)}+\epsilon ^{(i)})\frac{\stackrel{}{v}_\text{c}\stackrel{}{v}_\text{c}\stackrel{}{n}}{1\stackrel{}{v}_\text{c}^{\mathrm{\hspace{0.17em}2}}}.$$
(1)
The upper index $`(i)`$ refers for the intrinsic energy density $`\epsilon `$ and pressure $`P`$ of matter in the frame of reference, locally at rest, i.e. observed by a co-moving observer. We omit the superscript $`(i)`$ in the following. For the fireball expansion to continue, $`𝒫|\stackrel{}{𝒫}|>0`$ is required. For $`𝒫0`$ at $`v_\text{c}0`$, we have a conflict between the desire of the motion to stop or even reverse, and the continued inertial expansion.

When the flow velocity remains large but $`𝒫0`$, the intrinsic pressure $`P`$ must be negative. As illustration consider the fireball to be made of a quark-gluon liquid confined by an external vacuum pressure $``$. The total pressure and energy comprise particle (subscript p) and the vacuum properties:

$$P=P_\text{p},\epsilon =\epsilon _\text{p}+.$$
(2)
Eq. (1) with $`\stackrel{}{𝒫}=0`$ thus reads:

$$\stackrel{}{n}=P_\text{p}\stackrel{}{n}+(P_\text{p}+\epsilon _\text{p})\frac{\stackrel{}{v}_\text{c}\stackrel{}{v}_\text{c}\stackrel{}{n}}{1v_\text{c}^2},$$
(3)
and it describes the (equilibrium) condition where the pressure of the expanding quark-gluon fluid is just balanced by the external vacuum pressure.

Expansion beyond $`𝒫0`$ is in general not possible. A surface region of the fireball that reached it but continues to flow outwards must be torn apart. This is a collective instability and thus the ensuing disintegration of the fireball matter will be very rapid, provided that much of the surface reaches this condition. We adopt the condition $`\stackrel{}{𝒫}=0`$ at any surface region to be the instability condition of an expanding hadron matter fireball.

Negative internal pressure $`P<0`$ is a requirement. At this stage the fireball must thus be significantly supercooled. The adiabatic transfer of internal heat into accelerating flow of matter provides the mechanism which leads on the scale of $`\tau =\mathrm{2\hspace{0.17em}10}^{23}`$ s to the development of this ‘deep’ supercooling

It is possible to determine experimentally if the condition $`P<0`$ has been reached. Namely, the Gibbs-Duham relation for a unit volume:

$$P=T\sigma +\mu _\text{b}\nu _\text{b}\epsilon ,$$
(4)
relates the pressure, to entropy density $`\sigma =S/V`$ , energy density $`\epsilon =E/V`$ , and baryon density $`\nu _\text{b}=b/V`$ , $`V`$ is the volume, $`T`$ is the temperature, and $`\mu _\text{b}`$ the baryochemical potential. Dividing by $`\epsilon `$ we obtain:

$$\frac{PV}{E}=\frac{T_\text{h}}{E/S}+\frac{\mu _\text{b}}{E/b}1.$$
(5)
The microscopic processes governing the fireball breakup determine how the quantities entering the right hand side of Eq. (5) are changed as hadrons emerge. Understanding this we can determine, if the intrinsic fireball pressure prior to breakup, has been negative.

The energy $`E`$ and baryon content $`b`$ of the fireball are conserved. Entropy $`S`$ is conserved when the gluon content of a QGP fireball is transformed into quark pairs in the entropy conserving process $`G+Gq+\overline{q}`$. Similarly, when quarks and antiquarks recombine into hadrons, entropy is conserved in the range of parameters of interest here. Thus also $`E/b`$ and $`S/b`$ is conserved across hadronization condition. The sudden hadronization process also maintains the temperature $`T`$ and baryo-chemical potential $`\mu _b`$ across the phase boundary. What changes are the chemical occupancy parameters. As gluons convert into quark pairs and hadrons $`\gamma _g0`$ but the number occupancy of light valance quark pairs increases $`\gamma _q>\gamma _{q_0}1`$ increases significantly, along with the number occupancy of strange quark pairs $`\gamma _s>\gamma _{s_0}1`$.

The sudden hadronization picture differs from, e.g., the droplet-driven reequilibration transformation , in which chemical equilibrium of valance quark pair abundances is maintained. Instead a change (reheating) of statistical parameters $`T,\mu _b`$ occurs, along with a possible formation of a mixed phase required for volume expansion. In such reequilibration hadronization picture, entropy increase can also occur . To draw the line between these two hadronization pictures (equilibrium/sudden), we need to determine the quark pair occupancy parameters $`\gamma _i`$. Our study of final state hadron abundances strongly favors $`\gamma _q\mathrm{exp}(m_\pi /2T)>1,\gamma _s\gamma _i`$ . Evaluating Eq. (5) using the results of our data analysis, we indeed obtain $`P_f<0`$. The magnitude of $`|P_f|`$ can vary between a few percent (in terms of energy density $`E_f/V`$), up to 20% for the latest published result . The precise value, which arises from several cancellations of larger numbers is sensitive to the strategy of how the currently available experimental data is described, e.g., if strangeness conservation is implemented, and if so, if differentially at each rapidity, or as an overall conservation law; how many high mass resonances can be excited in hadronization process, etc …

Importantly, we have not been able to obtain a scheme of hadron production analysis which describes the data with ‘$`\chi ^2/\text{dof}<1`$ and would not imply $`P<0`$ for the hadronizing fireball matter. On the other hand, if we do force the hadronizing particles to be in chemical equilibrium, we find $`\chi ^2/\text{dof}>2.5,\text{dof}=10`$ in our analysis which agrees for this limit with , and in this case we find $`P>0`$.

Understanding in detail the breakup condition $`𝒫0`$ requires that we model the shape and direction of flow in the late stage of fireball evolution, obviously not an easy task. However, considering $`\stackrel{}{n}\stackrel{}{𝒫}0`$, we find the constraint:

$$\frac{PV/E}{1+PV/E}=\kappa \frac{v_\text{c}^2}{1v_\text{c}^2},\kappa =(\stackrel{}{v}_\text{c}\stackrel{}{n})^2/v_\text{c}^2.$$
(6)
For an exactly symmetrical, spherical expansion the two vectors $`\stackrel{}{v}_\text{c}`$ and $`\stackrel{}{n}`$ are everywhere parallel, thus $`\kappa 1`$. However, in 158$`A`$ GeV Pb–Pb reactions the longitudinal flow is considerably greater than the transverse flow , and we note $`\kappa 0`$ for a longitudinally evolving cylindrical fireball. For the Pb–Pb collisions considered here, our analysis suggest $`0.1<\kappa <0.6`$.

We now substitute, in Eq. (6), the fireball matter properties employing the Gibbs-Duham relation, Eq. (4), and arrive at:

$$\frac{E}{S}=\left(T_\text{h}+\frac{\mu _\text{b}}{S/b}\right)\left\{1+\kappa \frac{v_\text{c}^2}{1v_\text{c}^2}\right\}.$$
(7)
Eq. (7) establishes a general constraint characterizing the fireball breakup condition.

The solid line, in figure 1, shows the behavior of $`v_\text{c}(T_\text{h})`$ constraint arising from Eq. (7) for the example $`E/S=0.184\pm 0.05`$ GeV (error range shown by dotted lines), $`\kappa =0.6`$. Outside of the region bounded by the solid line (i.e., for greater $`T_\text{h}`$ and $`v_\text{c}`$), the flow expansion can occur as the internal particle pressure is greater than the confining pressure. Also shown in figure 1 is hadron production analysis result and its statistical error, the systematic error is of same magnitude. The agreement of theory and experiment results from the choice of non-spherical flow with specific freeze-out shape described by the average value $`\kappa =0.6`$, see Eq. (6). Figure 1 illustrates the great sensitivity to the analysis on the freeze-out constraint. The dashed horizontal line, in figure 1, is the velocity of sound of the interacting quark-gluon liquid, which barely differs from $`1/\sqrt{3}`$ .

We have so far not used in the discussion any key specific property of the equations of state of the matter filling the fireball. However, our results imply that the matter inside the fireball is deeply supercooled. Can this be the deeply supercooled liquid of quarks and gluons? In fact a study of QGP equations of state employing properties of QCD interactions and thermal QCD , fine tuned to agree with the properties of lattice QCD results  suggest that. We extend this study to consider the phase boundary. The thin solid line in the $`T,\mu _\text{b}`$ plane in figure 2 shows where the pressure of the quark-gluon liquid equals the equilibrated hadron gas pressure. The hadron gas behavior is obtained evaluating and summing the contributions of all known hadronic resonances considered to be point particles. When we allow for finite volume of hadrons , we find that the hadron pressure is slightly reduced, leading to some (5 MeV) reduction in the equilibrium transition temperature, as is shown by the dashed line in figure 2 . For vanishing baryo-chemical potential, we note in figure 2 that the equilibrium phase transition temperature is $`T_{\text{pt}}172`$ MeV, and when finite hadron size is allowed, $`T_{\text{fp}}166`$ MeV, The scale in temperature we discuss is result of comparison with lattice gauge results. Within the lattice calculations , it arises from the comparison with the string tension.

The dotted lines, in figure 2, correspond to the condition Eq, (3) using the shape parameter $`\kappa =0.6`$, Eq. (6), for (from right to left) $`v_\text{c}^2=0,1/10,1/6,1/5,1/4`$ and $`1/3`$. The last dotted line corresponds thus to an expansion flow with the velocity of sound of relativistic noninteracting massless gas. The thick solid line corresponds to an expansion with $`v_\text{c}=0.54`$. The hadron analysis result is also shown . Comparing in figure 2 thin solid/dashed with the thick line, we recognize the deep supercooling as required for the explosive fireball disintegration. The super-cooled zero pressure $`P=0`$ QGP temperature is at $`T_{\text{sc}}=157`$ MeV, (see the intercept of the first dashed line to the right in figure 2) and an expanding fireball can deeply super-cool to $`T_{\text{dsc}}147`$ MeV (see the intercept of thick solid line) before the mechanical instability occurs.

Deep supercooling requires a first order phase transition, and this in turn implies presence of latent heat $``$. Physical consistency then requires presence of external (negative) vacuum pressure $``$. More precisely, the vacuum contribution to the physical properties of deconfined matter can be derived from $`\mathrm{ln}𝒵_{\text{vac}}V\beta `$ :

$`P_{\text{vac}}`$ $`=`$ $`{\displaystyle \frac{T}{V}}\mathrm{ln}𝒵_{\text{vac}}=,`$ (8)
$`\epsilon _{\text{vac}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{ln}𝒵_{\text{vac}}}{V\beta }}=\left\{1+{\displaystyle \frac{\mathrm{ln}/_0}{\mathrm{ln}\beta /\beta _0}}\right\}.`$ (9)
The temperature, $`T=1/\beta `$, dependence of the vacuum pressure has been considered within the model of color-magnetic vacuum structure . Near to the phase transformation condition, the variation of $``$ with $`\beta `$ is minimal (see figure 2 in ), and thus the logarithmically small last term in Eq. (9) can be, in principle, ignored.

We now combine the theoretical properties of the QGP equations of state with the dynamical fireball properties in order to constrain $``$. Reviewing Eq. (5), we obtain:

$$\frac{PV}{E}\epsilon _{\text{QGP}}+P_\text{p}=,$$
(10)
To evaluate $``$, we note that lattice results for $`\epsilon _{\text{QGP}}`$ are well represented by $`\epsilon _{\text{QGP}}=aT^4`$, with $`a11`$, value extrapolated for the number of light quark flavors being $`n_f=2.5`$ at the hadronization point . We obtain, for the fireball formed in Pb–Pb reactions,

$$0.211T_\text{h}^40.17\text{ GeV/fm}^3.$$

Is our picture of fireball evolution compelling? We found that particle production occurred at condition of negative pressure expected in a deeply supercooled state and have shown internal consistency with (strange) hadron production analysis involving chemical non-equilibrium. Moreover, these chemical freeze-out conditions agree with thermal analysis , allowing the conjecture that the explosive quark-gluon fireball breakup forms final state hadrons, which do not undergo further reequilibration. However, we noted that the chemical equilibrium reaction picture differs from ours only in terms of its statistical significance ($`\chi ^2/\text{dof}>2.5`$, ). It produces a chemical freeze-out temperature of $`T=168`$ MeV, just the value we found for an equilibrium phase transition implicit in the assumption of chemical equilibrium. The higher chemical freeze-out temperature produces greater population of excited hadronic states. Their decays deform hadron spectra, and this allows for a second evolution stage with a thermal freeze-out temperature at or below $`120`$ MeV.

Strongly in favor of the here described sudden QGP hadronization resulting in chemical nonequilibrium reaction picture is the presence of a hadron multiplicity excess, related to an entropy excess . This is seen both as multiplicity per baryon, and as increase of multiplicity comparing $`pA`$ to $`AA`$ reactions. We could not describe this effect by admitting other physical models such as, in medium, change of hadron masses. We have found that invariably the statistical significance of the analysis decreases as we modify individual hadron properties in an ad-hoc fashion, while maintaining chemical equilibrium of hadron abundances. Our experience shows that the only theoretical description of hadron production data that works ($`\chi ^2/\text{dof}<1`$) requires excess of valence quark pair abundance, irrespective of the detailed strategy of data analysis. This agrees well with the dynamical study of nuclear collisions within the UrQMD model concludes that at CERN energies the chemical non-equilibrium is required to characterize the numerical results in terms of a statistical model .

Also in favor of our result is the conclusion of Csörgő and Csernai : who required as verification for the presence of a deeply supercooled state of matter and sudden hadronization: i) short duration and relatively short mean proper-time of particle emission, now seen in particle correlations , ii) clean strangeness signal of QGP ; iii) universality of produced particle spectra which are the remarkable features of strange particle production  v) no mass shift of the phi-meson; despite extensive search such a shift has not been found by the NA49 collaboration .

In summary, we have introduced a constraint, Eq. (7), which relates the physical and statistical properties of the hadronic fireball at the point of sudden breakup. We obtained this result from mechanical stability consideration, employing only properties of the energy-stress tensor of matter, and the Gibbs-Duham relation, Eq. (4). We showed that this constraint is consistent with analysis results obtained considering the experimental particle production data for Pb–Pb collisions at 158$`A`$ GeV. We further studied the behavior of the phase transition between hadron gas and quark-gluon liquid, and have determined the magnitude of the deep supercooling occurring in the fireball expansion. Employing a lattice-QCD based estimate on number of degrees of freedom in the energy density of the QCD thermal matter, we obtained a constraint on the magnitude of latent heat/vacuum pressure $`0.17\text{ GeV/fm}^3`$. We conclude that both in theoretical study of the data as well as for reasons of principle the deciding factor about the sudden nature of the phase transformation is the absence of chemical equilibrium.

We thank T. Csörgő for constructive comments. Work Supported in part by a grant from the U.S. Department of Energy, DE-FG03-95ER40937 . LPTHE, Univ. Paris 6 et 7 is: Unité mixte de Recherche du CNRS, UMR7589.


text:
# Multiagent Control of Self-reconfigurable RobotsA condensed version of this paper will appear in Proc. of the Intl. Conf. on Multiagent Systems (ICMAS2000), July 2000.

## 1 Introduction

Modular self-reconfigurable (or metamorphic) robots  consist of many simple identical modules, that can attach and detach from one another to change their overall topology. These systems can dynamically adapt their shape to suit the needs of the task at hand, e.g., for manipulation, locomotion, and adaptive structures responding to environmental stresses.

From a planning and control viewpoint, metamorphic robots pose several interesting research challenges. Self-reconfiguration, or how to change shape automatically, is a new and so far little studied problem for robots. Decentralized control is a useful approach to this problem especially when the robot has a the large number of modules, each of which is a self-contained unit with its own processing, sensing and actuation. With the tight physical interactions due to contact between neighboring modules and constraints arising from actuator geometry and power limitations, modular metamorphic robots pose an interesting challenge for multiagent control. They thus require a more physically grounded approach than many studies of artificial life  which tend to concentrate on how complex natural organisms achieve sophisticated crowd behaviors or deal with abstract agents that currently can not be physically constructed . Alternatively, biological and chemical techniques for distributed construction of shapes  provide examples of local behaviors producing complex shapes, but are not yet capable of producing general programmable robotic systems.

One approach to reconfiguration uses a precise specification of the desired locations of all the modules, and then solves the combinatorial search required to identify motions for the modules according to some criterion, such as minimizing the number of moves or power consumption. Such searches are generally intractable for robots consisting of many modules, but can be addressed approximately using heuristics . Unfortunately, in many practical applications, defining an exact target shape may not be suitable or even possible. This may arise when some modules fail or the nature of the environment or task is uncertain, for example, when grasping an object of unknown size or shape.

Instead, we use multiagent control to achieve suitable reconfiguration as a side-effect of creating a structure with the properties (structural, morphological, etc) required for the task. When the properties can be expressed largely in terms of the local environment for each module, agent-based control often achieves a suitable shape without any need to precisely specify the exact position of each module. Moreover, an agent-based architecture is well-suited to decomposing control problems based on the different physical phenomena dominant at different scales, especially for metamorphic robots with a large number of tiny modules. For example, micromachined robots  are dominated by friction and other surface forces rather than gravity. Even smaller structures  are subject to randomly fluctuating forces, i.e., Brownian motion. In contrast to larger robots, these tiny machines readily move objects many times their own weight and have comparatively high speeds and strengths . Biology has numerous examples of different structures appropriate for different scales . In such cases, different types of agents could be responsible for behaviors of individual modules, small groups of modules and so on, forming a hierarchical or multihierarchical correspondence between physical structures of the robot and its environment, the levels of specification for the desired task and the controlling software .

A related line of work applies multiagent control to teams of robots cooperating to achieve a common task . These methods usually apply to independently mobile robots that are not physically connected and have little or no physical contact. In most current modular robot systems the modules remain attached to one another forming a single connected whole , and giving rise to a number of tight physical motion constraints that do not apply to teams of independent robots. On the other hand, the physical contact between modules allows modules to locate their neighbors without complex sensory processing as would be required, for example, to visually identify a physically disconnected member of a team. Hence the techniques for coordinating teams are not directly applicable to modular robot controls. Furthermore, more traditional self-reconfiguration algorithms require an a-priori exact description of a target shape for the given task, which may not be suitable when the robot operates in uncertain environments.

Specifically, we explore the use of simple, purely local rules to produce control algorithms for accomplishing tasks such as dynamic adaptation under changing external conditions (e.g., added weight) and grasping objects. We will use the term behaviors for the specific control algorithms. The assumptions we make are as follows:

* Modules have limited computational capabilities, consisting of a limited memory and a simple finite-state machine (FSM). State transitions are driven by the states, the relative positions of a module and its neighbors, and some external sensor information. The state transitions can update the memory of a module and its neighbors.
* Communication is limited to immediate neighbors, and a limited number of bits are exchanged at each step. We specifically did not allow modules to broadcast messages globally, because the power required would likely not scale well if the size of modules shrinks. Morover, fully utilizing broadcasts would require more complex knowledge and processing on the part of the modules.

The remainder of this paper presents the design and evaluation of control algorithms for metamorphic systems that coordinate their actions locally to achieve emergent, global behaviors. The next section describes the robot platform used in our work and general control primitives. The following two sections present the results of our algorithms applied to specific tasks and provide an analysis of the behavior.

## 2 The Robot System

This section describes the modular robot used in our work, its simulator and the control primitives combined to create the behaviors described below.

### 2.1 Experimental Robot Platform

Several modular metamorphic robot designs have been proposed . For our study, we focus on one such design: Proteo , which is a modular self-reconfigurable, or metamorphic, robot being developed at Xerox PARC. In the current design, modules are rhombic dodecahedra (polygons with twelve identical faces, each of which is a rhombus). Modules attach to one another along their faces forming general three dimensional solids. To achieve a change of shape modules on the outer surface roll over the substrate of other modules to new positions. This type of reconfiguration has been called “substrate reconfiguration”, in contrast with other existing reconfiguration classes .

The rhombic dodecahedron can be thought of as the 3D analog of the hexagon and has several useful properties. Notably, these polyhedra result in maximium internal volume for a given surface area , meaning more room for packing electronic and mechanical components. Also, all module-on-module rotations are always 120 degrees, unlike the cube which requires 180 degree rotations in certain cases. A main disadvantage is that twelve faces per module need to be actuated, increasing the complexity and expense of the hardware. For a target module size in the centimeter scale, current actuator technologies result in large weight to power ratios and high cost per module. A manually actuated version of Proteo has been built and is shown in Figure 1. Modules connect along their faces, as illustrated in Figure 2.

In addition to actuators, the modules can communicate directly with their neighbors. Depending on the application, they can also include various sensor elements, e.g., to detect forces imposed by the environment or the weight of other modules.

For the experiments reported below, we used a simulation of the Proteo robot platform. This simulation includes the physical motion constraints on the modules, allowing them to rotate onto neighboring modules only when that motion is not obstructed by other modules. This simulation has the modules operating asynchronously. Specifically, every module is given a chance to execute its behavior exactly once each time cycle. The behavior executions for the modules are ordered in a different, random way for each time cycle. When the behavior code is executed each module describes what movement (if any) it wants to perform at the end of the cycle. After all the behavior code is executed, the simulator attempts to perform the requested movement for each module in turn, subject to geometrical constraints. This means that a module might be denied the movement it has requested.

The control programs in the modules are not directly notified of the outcome of their movement requests (by means of an exception, for example). While this capability could certainly be added to the system, we found it is not necessary for applying our coordination techniques to the tasks described below. This lack of notification means modules can become “stuck” and attempt the same movement in vain for several time cycles. However, the overall behavior is still achieved because at the same time other modules are able to complete their moves, and can eventually remove the obstacle preventing motion. Alternatively, the motion of these other modules can sufficiently change the structure that the stuck module decides on another motion.

Another programming issue was how a module should sample randomly from the available moves. When a completely random move is desired, the module simply picks with equal probability from the available moves to neighboring positions. When a directed random move is needed, a different method is used, which picks only from the available moves that go in a direction that has a positive dot product with a specified general direction in space. The latter approach is used when there is some bias direction inherent to the specific application (e.g. the direction of the “ball” in the grasping behavior described in Section 3.4). While such biased choices do not necessarily result in the most direct motion toward the goal, the randomness reduces the likelihood of becoming stuck and increases robustness in case of failed modules. From a global planning perspective, this randomness can be viewed as a simple search procedure that avoids the potentially intractable combinatorial search of finding an optimal set of motions for all the modules.

The simulation does not model all the details of the hardware. For example, it does not include limits on power consumption or actuator strength that could be significant in actual hardware implementations. The simulator also assumes the modules operate correctly and have accurate local sensory information, e.g., concerning contact with neighbors or other objects or mechanical stress. The results of the simulation are illustrated using the schematic representation of the modules included in Figures 1 and 2.

### 2.2 Control Primitives

While the Proteo system illustrates the issues facing modular robotics and provides a useful testbed, our methods are not tied to any specific design, but instead are applicable to metamorphic robots in general. Specifically, the agents make use of a few primitive constructs that are combined to form a control program for each module. These primitives allow simple communication and coordination among the modules, and could be implemented by a variety of modular robots.

Our approach uses the following primitives:

Growth is the process that creates structures. It is important to provide a mechanism for focusing the movement of modules to specific spots. Otherwise, if the movement is too random, it could take a long time before any reasonably good structure is grown.

Seeds are the main agents that cause growth. A seed is a module that attracts other modules in order to further grow the structure. As more modules are attracted to a seed, they can in turn become seeds themselves, and thus propagate the growth process.

Scents are the means of global communication among modules. Scents are propagated through the system in a distributed breadth-first fashion as follows. Each module keeps track of the scent strength at its location with a value in its memory. Those modules that emit the scent set their own value to zero. Other modules examine the scent values of their neighbors at each time step, and then set their own scent value to one more than the minimum value among their neighbors. Thus, smaller values indicate a stronger scent, and a scent gradient is created throughout the structure. Scent values are an approximation of the minimum distance to a scent-emitting module in the system at any given time. The approximation arises because scents are broadcast via neighbors and thus can take many steps to propagate throughout the entire system, and during this time modules can move. Nevertheless, as discussed further in Section 4, the scents generally propagate much faster than substantial rearrangements of the modules so this approximation worked quite well. We typically use either one or two different scents in our experiments. Seeds emit a scent to attract modules, and modules that are searching for seeds move along the surface gradient of the scent to find them.

The Mode of a module is its present FSM state. For instance, seeds are usually denoted by the SEED mode, modules that search for seeds are in the SEARCH mode, modules that are part of the finished structure are in the FINAL mode, and modules that cause branching are in NODE mode. The mode in turn will determine the rules of behavior of a module.

A final component of our control method is the use of random motions when different choices appear equally useful. This allows modules to continue performing even when the scent gradients do not completely specify the best direction. Because such random motions are less efficient than following gradients, our algorithms arrange for relatively few states to search before useful gradient information is found.

A simple example is shown in Figure 3. A single SEED module, with scent $`S=0`$, is connected to a chain of modules in SEARCH mode with successively larger scent values. In this configuration, only the SEARCH module at the end of the chain is free to move. Its neighbor’s smaller scent value (i.e., $`S=2`$) indicates that the module at the end of the chain can find a position with even smaller scent value by rolling over the faces of its neighbor. One possible sequence of steps following from this initial configuration is shown in Figure 4. For simplicity, this figure shows the remaining modules remaining stationary.

In this example, the moving module goes directly toward the seed. However, the gradient does not specify precisely which other face is best. Instead, the free module moves randomly over the available faces of its neighbor until it encounters the module with $`S=1`$. At this point, the moving module updates its scent value to $`S=2`$ and then continues its search for the seed by moving over the faces of the $`S=1`$ module. Because there are only a few faces to explore before new gradient information is found, this use of random search is likely to complete after only a few moves.

## 3 Results

This section presents experiments on Proteo with some environmental sensing and using our agent-based control algorithms within each module. The experiments were restricted to no more than a few hundred modules to keep simulation times reasonable (i.e., at most several hours).

### 3.1 A Chain

A chain is often the preferred locomotion configuration of modular robots for moving over steps, snaking into holes or squeezing through narrow passages. We describe a simple local algorithm that, starting from an arbitrary connected configuration, creates a single, one-module-thick chain. Initially, all modules are in the SLEEP mode. The other possible modes are: SEARCH, FINAL, and SEED. A module is arbitrarily picked to serve as the initial seed to start the chain. This initial seed can be determined centrally, or each module in SLEEP mode can be given a small probability to become a seed at each step. This module will pick a direction of growth and emit a scent to start attracting other modules toward it. When a module moves next to the seed module and is connected to it in the chosen direction, the chain grows by one module and that new module then becomes the new seed that will further grow the chain.

The control rules for each module are:

* If in SLEEP mode, if a scent is detected, go to SEARCH mode.
* If in SEARCH mode, propagate scent and move along scent gradient.
* If in SEED mode, emit scent, and if a module has appeared in the direction of growth, set that module to SEED mode, and go to FINAL mode.
* If in FINAL mode, propagate scent.

These rules amount to a simple finite-state machine for each module, shown in Figure 5. An example configuration of a group of modules following these rules is given in Figure 6. As the algorithm continues, the modules form a single chain.

### 3.2 Recursive Branching for Locomotion and Manipulation

A variety of applications for modular robots can benefit from a branching structure with “limbs” at several levels. In such a hierarchical structure, each level has varying degrees of precision, range of motion and strength. For example, one level of branching gives a structure that can be used as an artificial hand or as a “spider” for locomotion. Adding extra levels of branching results in “fingers” or “toes” with lower strength and range of motion but higher level of precision. Human limbs exhibit only two levels of branching, but a robot could grow as many levels as needed to achieve increasingly complex manipulation and locomotion tasks .

The algorithm is as follows. All modules are initially in the SLEEP mode. Growth is initiated when a module randomly decides to switch to the NODE mode. In this mode it attracts other modules (emits a scent) and randomly chooses which ones will become seeds. When it has spawned a certain number of seeds, the node becomes inactive. This experiment uses two scents: a “regular” scent is used to grow the structure (as in the Chain example), and a second “node” scent is used by modules to determine how far the nearest node is. If this distance is too big, a module can decide to become a node. The rules are:

* If in SLEEP mode, if regular scent is detected, go to SEARCH mode.
* If in SEARCH mode, propagate both scents and move along the regular scent gradient.
* If in SEED mode, propagate the node scent, and emit a regular scent; if there is a neighboring module in the direction of growth of the branch, set that module to be a seed, and go to FINAL mode.
* If in FINAL mode, propagate both scents; if the node scent is weak (i.e., has a value greater than some threshold $`T`$), go to NODE mode.
* If in NODE mode, emit both scents, spawn seeds in random directions, until a certain count is reached $`B`$, then go to INODE mode.
* If in INODE mode, emit a node scent and propagate the regular scent.

An example structure resulting from these rules is shown in Figure 7. The threshold $`T`$ determines the distance between nodes, i.e., how far modules continue building one branch before starting new branches. The count $`B`$ determines the branching ratio throughout the structure.

### 3.3 Dynamically Adapting to External Forces

One important task for reconfigurable robots is to adjust themselves in response to environmental forces. For instance, a collection of such robots supporting a weight on a set of “legs” should be able to change the location and density of the legs when the weight shifts. As a simple example of this task, we examined the formation of legs to support a flat structure with an additional imposed force whose location could change over time. For simplicity, we neglect the weights of the modules themselves, i.e., we assume their weights are small compared to the weight of the additional objects they are supporting on top of the flat structure or “table”.

For this behavior, we form some of the modules into a “board” (the top of the table) which supports the additional weight. These modules, which never move, are of two types. The first, always in FIXED mode, only transmit scent. The second type are “root” modules which are uniformly spaced over the board and can communicate with each other, e.g., through signals sent locally through the fixed modules. We assume these signals propagate much faster than the physical movements of the modules. The roots are grouped into regions on the board (in the experiments reported here, we use two regions: the two halves of the board).

We assume modules are equipped with a force sensor that allows them to measure the weight they are supporting. Furthermore, we assume the modules can either determine the direction gravity acts or this direction toward the ground is prespecified when creating the roots in the board. Roots within a region monitor the total weight supported in that region by communicating their sensor readings.

The root modules on the board can be in one of three modes: ROOT, IROOT or AROOT. Initially the root modules are in IROOT mode, which monitor the weight they are supporting and may start growing a leg. AROOT mode emits a scent to attract other modules to start growing a leg, in much the same way as growing a single chain. Once the leg starts growing, the root shifts to ROOT mode, where it remains until the weight it supports drops below a specified threshold in which case it probabilistically causes its leg to disband. Thus the root modules grow or disband legs probabilistically according to how much weight is experienced in their part of the table. We refer to root modules in ROOT or AROOT mode as “active” roots.

The fixed and root modules are set at the beginning and don’t change throughout the experiment. All other modules are initially in SLEEP mode. The structure first grows legs, then disbands some of them and creates new ones according to shifts in weight (controlled interactively by the user of the simulator).

The rules for this behavior for the root modules are:

* If in IROOT mode, if average weight per active root in its region is above a certain threshold $`F_{\mathrm{max}}`$, with a small probability $`p_{\mathrm{max}}`$: go to AROOT mode.
* If in AROOT mode, emit scent.
* If in ROOT mode, if average weight per active root in its region is below a certain threshold $`F_{\mathrm{min}}`$, with a small probability $`p_{\mathrm{min}}`$: set all FINAL and SEED neighbors to DISBAND mode, go to IROOT mode.

Although similar adaptive behaviors can occur with a range of values for the thresholds, one important issue is to prevent the disbanding of a single leg supporting a minimal weight in a region. In our experiments, this minimal value is one unit of weight, so we should have $`F_{\mathrm{min}}1`$. The rules for the remaining modules are:

* If in SLEEP mode, if a scent is detected, go to SEARCH mode.
* If in SEARCH mode, propagate scent and move along its gradient; if there is a neighboring SEED or AROOT module and location is towards the ground with respect to this neighbor

+ if this neighbor is a SEED, set it to FINAL mode, otherwise set it to ROOT mode
    + go to SEED mode
* If in SEED mode, if touching the ground, go to FINAL mode.
* If in FIXED or FINAL mode, transmit scent.
* If in DISBAND mode, set all FINAL and SEED neighbors to DISBAND mode, go to SEARCH mode.

Figure 8 illustrates how, once created, a supporting structure can adapt in response to changes in external forces. The behavior dynamically adapts the location and number of its legs to accommodate changes in the supported weight. In this example, we suppose weight $`W_1`$ is applied to the first half of the table and $`W_2`$ applied to the second half, with an arbitrary choice of units giving $`W_1+W_2=10`$. In each half of the table, active roots communicate their sensory information to determine how much of this weight is supported by those in their half of the table. The result is the average weight supported by each active root in each half of the table, i.e., $`w_1=W_1/N_1`$ for the first half (where $`N_1`$ is the number of active roots in the first half), and similarly for $`w_2`$ as the average weight supported per active root on the second half.

The probabilistic growing and disbanding of legs avoids possible oscillations in root behavior: if deterministic, when the weight shifts, many roots could decide to grow legs, then disband them on the next step, thus oscillating between the two without doing anything useful. Such oscillations are common for systems with synchronous updates : randomization is a simple technique to prevent spurious synchronization of agent activity.

After a root $`r`$ decides to grow a leg, some time is required for the modules to move to that root and produce the leg. However, because a root’s decision to grow a leg is based on the weight in its region averaged over the active roots, once root $`r`$ starts the process of growing a leg by emitting scent, the average weight per active root drops immediately and proportionately with the number of supporting modules. This change, rapidly communicated among the roots, gives feedback to other roots in the region and prevents other legs from growing nearby unnecessarily. In practice, this implies that the weight must be shifted slowly: the system will not be able to respond quickly to sudden changes.

An indication of how fast the modules react to external conditions is the number of steps after a weight shift until a new stable configuration is achieved. In the examples shown here, this time is on the order of 200 cycles. This includes the time needed to disband some legs and construct some new ones. Although the modules know about the shifted weight immediately, some time is needed until some roots decide to disband their legs, and some inactive roots decide to grow legs. These decisions are probabilistic: each root decides to change its state (go from active to inactive or vice versa) with a probability we chose to be on the order of $`1/R`$ where $`R`$ is the total number of roots in the structure. Thus it is unlikely that many roots will decide to switch state at the same time step, thereby avoiding unstable or oscillating system behavior. We depend on the rough assumption that decisions to disband or grow legs have immediate effect on the strain felt by the modules in the system. This assumption is questionable for real-world applications, especially with larger structures. However, it provides a good starting point. In real-world applications the trade-off will be between fast response with a risk of oscillating behavior and a slow, but stable response to external conditions.

### 3.4 Grasping an Object

The previous example showed how simple agent-based control allows the robot to respond to changes in environmental stress. In practice, such robots will be used to perform a globally specified task, leading to the question of how high-level, imprecise specifications can be combined with low-level behaviors.

As an example of such a task, we examined the behaviors required for Proteo to reach out and grasp an object of unknown shape and size and with only roughly specified location by “growing” around it. In this behavior the only modules that do substantial work are the modules in SEARCH mode. Initially all modules are in SLEEP mode, except for eight modules that are set to SEED mode (to create eight branches to grow toward the object). Two types of seeds are used: SEED modules, and TOUCHSEED modules. Initially growth is caused by SEED modules, but once the object is reached, seeds transfer to the TOUCHSEED state. The two types of seeds allow modules in SEARCH mode to know whether to grow in the direction of the object (to reach it) or to grow the structure around the object (to grasp it). The system is initially given the approximate direction the object is in, however, it has no knowledge of the exact direction of the object, how far it is, or what shape and size it has. The modules are assumed to have contact sensors to detect when the object has been reached. The rules are simple:

* If in SLEEP mode, if scent is detected, go to SEARCH mode.
* If in SEARCH mode, propagate scent, and move along the scent gradient; if a seed has been found: if touching the object, then set seed to TOUCH mode, and go to TOUCHSEED mode; otherwise, if not touching the object, and seed is in SEED mode, find a spot next to the seed in the direction of the object, set seed to FINAL mode, and go to SEED mode; otherwise, make a random move.
* If in TOUCHSEED or SEED mode, emit scent.
* If in TOUCH or FINAL mode, propagate scent.

This example incorporates a global goal, specified by the general direction of the object, with local adjustments by the agents to the actual object surface. It thus illustrates an important decomposition of the control task into global specification in terms of general parameters combined with low-level control based on the actual environment of the modules. This decomposition considerably simplifies the overall motion planning since the global specification need not have access to all the details of the environment, knowledge of which modules might have failed, etc.

These two examples could be combined: including stress sensors could also allow the modules to adjust for the object weight if it is to be picked up rather than just grasped.

## 4 Local Minima and Stability

Local control has the danger of getting trapped in locally “comfortable” configurations that are not globally suitable. Our agent approach provides good solutions to this problem both for individual modules and the system as a whole.

A module searches for a seed by moving along a gradient of least-descent. In Proteo, modules can only move by “rolling” over the outside surface of the structure, while scents can propagate through the inside of the structure. Usually this is not a problem: all modules searching for seeds are constantly moving preventing any fixed local extremum. However, a problem may arise when some modules remain fixed. Modules searching for seeds may not be able to follow the internal scent and get through to seeds on the other side because fixed modules effectively block the way. There are two possible solutions: one is to alter the design to allow modules to pass through other modules (by squeezing through), and the other is to use only a surface scent, which the modules can follow easily.

At a higher level, the whole system is performing a search for a configuration that suits a given criteria. If the reconfiguration rules are not set up carefully, however, the system may never converge to a suitable state. For instance, if for each module we only specify preferred immediate neighbor configurations, the most likely result is that surface modules will settle quickly and comfortably, preventing modules on the inside of the structure from moving. Thus, something must guide the overall reconfiguration. We use the concepts of growth and seeds to drive the reconfiguration of the whole system in a focused way.

We found the structures grown are generally stable. During growth, stability depends mostly on the ratio of scent propagation speed to module speed. Since modules generally move slower than they can transmit scent, the response to scents does not introduce oscillating behavior  or irregular growth, and the intermediate configurations are well-balanced at all times.

## 5 Conclusions

We have presented an agent-based approach to self-reconfiguration for modular robots that results in the creation of “emergent” structures with the desired functionality. We use purely local, simple rules and limited sensing. To “grow” stable structures, modules are guided by local attractors (seeds) and global gradients (scents) toward good configurations. Different structures result by varying the number and combination of seeds and scents in the algorithm.

Because the method has a strong random component, the specific shape of the resulting structures is non-deterministic. However, the resulting structure has the desired functionality. This is particularly well-suited to tasks with a strong element of uncertainty in the precise nature of the task or in the environment, in contrast to other reconfiguration algorithms.

We examined two tasks for which the reconfigurable nature of the robots is well suited: adjusting structures to environmental stresses and grasping an object whose detailed shape is not known. The global behavior of the algorithms we developed for these tasks was determined through simulations of the Proteo robot. The control primitives can easily be scaled up to a larger number of modules, and down in the size of individual modules, and are general enough to fit most modular robot designs.

This approach could be used in conjunction with other self-reconfiguration or control methods, as part of an overall hierarchical control scheme to handle increasingly complex tasks. Lately, an important research effort is applying genetic algorithms to the control of autonomous systems  or FPGA programming . The techniques described in this paper are amenable to evolutionary approaches, and could benefit from these results. Specifically, a population of agent programs combining the control primitives in different ways could be tested against variations in the desired task to evolve better behaviors according to various performance metrics . In the context of the primitives presented here, genetic techniques could identify methods to modulate the various parameters in our methods, e.g., the probabilty to form a seed for new growth.

Finally, improvements in hardware design to augment current module capabilities could greatly aid reconfiguration. Scaling down the size of individual modules using micromachines (MEMS)  or molecular assemblies could benefit actuation. For instance, protein motors can carry weights orders of magnitude larger than their own, at speeds that are impressive given their size . This means a module could carry several others around as it moves, something current Proteo modules cannot do. In addition, designs that allow individual modules to change their shape and size could also enhance overall performance. These improvements, however, would require truly interdisciplinary expertise in biology, materials, and electronics and remain, for the most part, a formidable challenge. Nevertheless, as such robotic systems are developed, multiagent control approaches are likely to help them achieve robust behaviors.

### Acknowledgments

We thank Ying Zhang who wrote the original Proteo simulator used in our experiments.


text:
# Observation of Superfluid Flow in a Bose-Einstein Condensed Gas

\[

## Abstract

We have studied the hydrodynamic flow in a Bose-Einstein condensate stirred by a macroscopic object, a blue detuned laser beam, using nondestructive in situ phase contrast imaging. A critical velocity for the onset of a pressure gradient has been observed, and shown to be density dependent. The technique has been compared to a calorimetric method used previously to measure the heating induced by the motion of the laser beam.

\]

Beginning with the London conjecture , Bose-Einstein condensation has been considered crucial for the understanding of superfluidity. Since then, the weakly interacting Bose gas has served as an idealized model for a superfluid . It has a phonon-like energy-momentum dispersion relation that does not allow for the generation of elementary excitations below a critical velocity, thus implying dissipationless flow at lower velocities. The onset of dissipation has been treated in the framework of the nonlinear Schrödinger equation  and shows an intriguing richness. Experiments on liquid helium could not test these theories, because superfluidity and dissipation are even more complex in this system due to the presence of strong interactions and surface effects .

The creation of Bose-Einstein condensates in dilute gases has dramatically changed this situation, allowing for quantitative tests of microscopic theories using the tools and precision of atomic physics experiments . A number of recent experiments have examined phenomenological features of superfluidity in gaseous Bose-Einstein condensates. These include the observation of vortices , a non-classical moment of inertia , and suppression of collisions from microscopic impurities . In previous work we found evidence for a critical velocity in a stirred condensate . In this Letter we study the onset of dissipation with higher sensitivity using repeated in situ non-destructive imaging of the condensate. These images show the distortion of the density distribution around the moving object, thus directly probing the dynamics of the flow field that has been recently treated with different models . The experimental setup was similar to the one used in our previous work . Improvements in the evaporation strategy and decompression techniques allowed us to produce pure condensates with up to 5$`\times 10^7`$ sodium atoms, with densities ranging from 8.4$`\times 10^{13}`$ to 3.5$`\times 10^{14}`$ cm-3, corresponding to chemical potentials from 60 to 250 nK. We determined the Thomas-Fermi radius $`R_z`$ along the axial direction through in situ phase contrast imaging. The sound velocity at the center of the condensate was then evaluated through the relationship $`c_s=2\pi \nu _zR_z/\sqrt{2}`$, with $`\nu _z`$=20.1 Hz being the axial trapping frequency. The macroscopic moving object was a 514 nm laser beam blue-detuned with respect to the sodium transitions, thereby creating an effective repulsive potential for the atoms. The beam is focused on the center of the condensate to a Gaussian $`1/e^2`$ diameter of $`2w=10\mu `$m. It is scanned with an acousto-optic deflector driven by a triangular wave along the axial direction of confinement of the magnetic trap. This results in a motion of the laser beam at constant speed $`v=2af`$, where $`f`$ is the frequency and $`a`$ is the amplitude of the scan. To avoid effects due to the inhomogeneous density distribution in the condensate, the laser beam was scanned at various frequencies with a constant amplitude approximately equal to the Thomas-Fermi radius $`R_z`$ of the condensate along the axial direction.

The frequency range was limited to 200 Hz. Above this frequency the atoms experience an effective time-averaged stationary potential. The condensate was imaged using the phase contrast technique  with probe light red-detuned from the F=1 to F=2 transition. Different detunings (0.4 and 1.7 GHz) allowed us to maintain a large linear dynamic range when imaging condensates at different densities. Fig. 1 shows multiple phase contrast images of the hydrodynamic evolution of the density distribution of a single condensate during the stirring. For low velocities there was little effect on the condensate. However, at higher velocities the density distribution became asymmetrical and the condensate piled up in front of the moving beam, resulting in turbulent flow around the object.

The main result of this Letter is based on quantitative analysis of images like those in Fig. 1. The analysis of the observed flow pattern shows a critical velocity for the onset of a drag force between the macroscopic laser beam and the condensate.

For a weakly interacting Bose-condensed gas at density $`n(𝐫)`$ and chemical potential $`\mu `$, pressure is identical to the mean-field energy density $`P=\mu (𝐫)n(𝐫)/2`$ . A drag force arises due to the pressure difference across the moving object. The chemical potential is given by $`\mu (𝐫,t)=gn(𝐫,t)`$, where $`g=4\pi \mathrm{}^2a/M`$ is the strength of two-body interactions for atoms with mass M and scattering length $`a`$. Therefore the magnitude of the drag force $`F`$ is given by

$$FgSn𝚫n=S\mu \mathrm{\Delta }\mu /g$$
(1)
where $`\mathrm{\Delta }n`$ and $`\mathrm{\Delta }\mu `$ are the differences in density and chemical potential across the stirring object, and $`S`$ the surface area the object presents to the condensate.

If the laser beam is stationary, or moves slowly enough to preserve a possible superfluid state for the condensate, there will be no gradient in chemical potential across the laser focus, and therefore zero force according to Eq.(1). A drag force between the moving beam and the condensate is indicated by an instantaneous density distribution $`n(𝐫,𝐭)`$ that is distorted asymmetrically with respect to the laser beam. Fig. 2 shows this asymmetry as the formation of a bow wave in front of and a stern wave behind the moving laser beam.

To disentangle asymmetries due to the moving beam from the intrinsic inhomogeneity of the density profile we strobed the images at half the stirring period and always viewed the laser beam in the center of the condensate. We defined the asymmetry $`A`$ as the relative difference between the peak column densities in front ($`\stackrel{~}{n}^f`$) and behind ($`\stackrel{~}{n}^b`$) the laser beam $`A=2(\stackrel{~}{n}^f\stackrel{~}{n}^b)/(\stackrel{~}{n}^f+\stackrel{~}{n}^b)`$. We averaged $`A`$ over 9 successive images of a single condensate.

The average asymmetry versus velocity is shown in Fig. 3. Below $`0.6`$ mm/s the asymmetry was zero within the statistical error indicating a suppression of dissipation in the superfluid regime. The asymmetry increased for velocities greater than 0.6 mm/s. The dip near 1.5 mm/s corresponds to a stirring frequency of $``$ 20 Hz for the chosen amplitude. This is close to the axial frequency of the magnetic trap. Near this frequency the laser beam synchronously excited the axial dipole motion of the condensate (see Fig. 3), thus reducing the relative velocity of the stirrer and the condensate and therefore the drag. We repeated the measurements at different stirring amplitudes to disentangle velocity and frequency dependent effects. This confirmed that the observed dip always occurred near 20 Hz, whereas the other features of the observed asymmetries depended only on velocity.

To study the superfluid regime at different atomic densities, we focused on the low velocity region ($`v`$ 2 mm/s) by choosing the stirring amplitudes in such a way that the frequencies were below the trap resonance. The density was varied by adjusting the number of atoms for constant trapping parameters. The power in the laser beam was scaled to keep the ratio of the potential height of the laser to the chemical potential equal to 8 within 10%, thus preserving the effective size of the laser beam.

In Fig. 4 we show measurements of the average asymmetry for two maximum densities $`n_0`$ of $`9\times 10^{13}`$ and $`1.9\times 10^{14}`$ cm-3. Two features are evident. In each data set there is a threshold velocity $`v_c`$ below which the drag force is negligible, and this threshold increases at higher density. Above this critical velocity, the drag force increases monotonically, with a larger slope at low density. The data were fitted to a piecewise linear function that reflects the two regimes of dissipation, one with zero asymmetry below $`v_c`$ and another, above $`v_c`$, where the asymmetry is linear in $`vv_c`$. The fits yield $`v_c=(0.11\pm 0.02)c_s`$ for the higher density and $`v_c=(0.07\pm 0.02)c_s`$ for the lower density. This is in agreement with the expectation that the critical velocity increases with the speed of sound and disagrees with models that predict only a dependence upon the diameter of the object .

One can compare measurements of the asymmetry (proportional to the drag force $`\stackrel{}{F}`$) with the power transferred to the condensate by the moving beam, $`\stackrel{}{F}\stackrel{}{v}`$, using an improved implementation of the calorimetric technique introduced in . For this, we stirred the condensate for times between 100 ms and 8 s, in order to produce approximately the same final temperature. After the stirring beam was shut off, the cloud was allowed to equilibrate for 100 ms. The thermal fraction was determined using ballistic expansion and absorption imaging . We inferred the temperature and total energy using the specific heat evaluated as in . A further improvement was obtained by comparing successive images with and without stirring of the condensate, to subtract out all other sources of residual heating. As a result of the improved fitting and background subtraction procedure, we were able to detect changes of less than 10 nK in the energy of the gas, as we will describe in more detail in . A condensed fraction higher than 90% ensured that the laser beam primarily heated the condensate and not the thermal cloud.

Fig. 5 shows the energy transfer rate to the condensate versus the stirring velocity. Here, as in the phase-contrast measurements, the amplitude of the stirring was kept fixed while scanning the frequency. We see the onset of dissipation near the value obtained using the drag force method, $`v_c`$ 0.5 mm/s = 0.1 $`c_s`$.

The calorimetric measurements can be compared with the drag force inferred from the asymmetric density distribution. Using Eq.(1), the energy transfer rate per atom is written in terms of the asymmetry as

$$\frac{dE}{dt}=\frac{\stackrel{}{F}\stackrel{}{v}}{N}\frac{\mu _0\stackrel{~}{n}_0d}{N}vA(v)$$
(2)
Here we have used $`S=l_zd`$ for the effective surface area of the object, with $`d`$ the diameter of the laser beam, $`l_z`$ the effective column depth, and $`\stackrel{~}{n}_0=n_0l_z`$ for the observed column density.

The comparison between the calorimetric and the drag force measurements is also shown in Fig. 5, with an adjusted scale factor for $`vA(v)`$. The measurements of force (as an asymmetry in the instantaneous condensate wavefunction) and heating (as an integrated increase in thermal fraction) are in remarkable agreement and demonstrates the consistency between our two methods of analyzing dissipation. For the parameters of our experiment ($`d10\mu `$m, $`n_0=1.3\times 10^{14}`$ cm-3, $`l_z=66\mu `$m, $`N=1.810^7`$) the overall heating rate predicted by Eq.(2) is 4.5 times larger than that obtained from calorimetry. This discrepancy is not surprising, because Eq.(2) does not account for dynamical pressure, the inhomogeneous density profile, and the 3D nature of the flow field.

In one model , the drag force arises from the periodic emission of vortex lines at a rate that increases with velocity. The vortices reduce the pressure gradient across the object and the predicted heating rate is

$$\frac{dE}{dt}\frac{1}{N}\frac{\mu }{h}\frac{v}{c_s^2}(vv_c)ϵ_{vortex}$$
(3)
where $`ϵ_{vortex}=2\pi \stackrel{~}{n}\mathrm{}^2/M\mathrm{ln}(d/\xi )`$ is the vortex energy and $`\xi =(8\pi na)^{1/2}`$ is the healing length. Equating this to the heating rate in Eq.(2), predicts that the slope of the asymmetry $`\alpha =A(v)/(vv_c)`$ should scale as $`\alpha 1/c_s^2`$. For the two densities in Fig. 4, the model predicts the ratio $`\alpha _{low}/\alpha _{high}`$ = 2.1, while our measurement yields $`(1.9\pm 0.2)`$. This, and the absolute value of the observed heating rate (see Ref. ), are further evidence that the drag force is closely linked to vortex shedding.

Our new measurements are consistent with our earlier experiment  where dissipation at subsonic velocities was observed. However, the scatter in the observed heating rates in  prevented a clear observation of the threshold, and the critical velocity was determined by linear extrapolation from high heating rates to be $`v_c/c_s=0.25\pm 0.04`$. Other possible extrapolations tend to lower the value of $`v_c`$, for instance a function $`v(vv_c)`$ as suggested by Eq.(2), yielded $`v_c/c_s=0.17\pm 0.09`$. With improved calorimetry we have taken data for parameters similar to  in the range 0.1 $`<v/c_s<`$ 0.25 and observed a small heating rate ($``$ 20 nK/s at 1 mm/s) that was previously indiscernible.

Since our first report on critical velocities in a Bose Einstein condensate, several theoretical papers have added further insight . The observed critical Mach number $`v_c/c_s0.1`$ is lower than the predictions for homogeneous  and inhomogeneous  2D systems. This discrepancy is most likely due to the 3D geometry of our experiment, where the laser beam pierced lower density regions of the condensate. It was noted in Ref.  that the critical velocity for phonon excitation is lowered by the inhomogeneous density distribution. Ref.  shows that vortex stretching and half-ring vortices can lower the 3D critical velocity below the 2D value. Jackson et al.  performed 3D simulations for a geometry similar to our experiment and obtained a critical velocity as low as $`0.13c_s`$, in good agreement with our results. The relevant critical velocity in our experiment is most likely related to vortex nucleation , which is usually smaller than the Landau and Feynman critical velocities  at which phonons  or vortices  become energetically favorable.

In conclusion, we have demonstrated a new technique to study the onset of dissipation in a stirred condensate. This method directly measures density distributions and pressure differences. The results obtained are consistent with calorimetric measurements but are more immune to shot-to-shot noise because multiple non-perturbative images of the flow field are taken. This technique is well suited to shed further light on contradictory predictions about the onset of dissipation in a gaseous superfluid  and to study in detail how the drag force depends on condensate size and density, as well as the object size.

We thank A. Görlitz, C. E. Kuklewicz and Z. Hadzibabic for contributing to the early stages of this experiment, and D. Kokorowski for a critical reading of the manuscript. This work has been supported by ONR, NSF, ARO, NASA, and the David and Lucile Packard Foundation. A.P.C. acknowledges additional support from NSF.


text:
# ENTROPY GENERATION IN COSMOLOGICAL PARTICLE CREATION

## Acknowledgements

This work was partially supported by Universidad de Buenos Aires.