text:
# Parametric Resonant Phenomena in Bose-Einstein Condensates: Breaking of Macroscopic Quantum Self-Trapping

## I Introduction

Macroscopic quantum tunneling with dilute Bose-Einstein condensates of alkali-metal atoms has been the subject of many theoretical \[1-4\] and experimental \[5-7\] studies. Here we consider the problem of trapped Bosons in a double-well potential and investigate the effect of a periodically varying barrier in the tunneling of the Bose condensate. We study the problem using the two-mode classical-like equations  and find that remarkable effects are related to parametric resonance  both in the periodic macroscopic quantum tunneling regime and in the macroscopic quantum self-trapping regime.

## II Tunneling in a double-well potential

In a dilute gas of $`N`$ Bosons at zero temperature practically all particles are in the same single-particle state of the density matrix $`\rho (𝐫,𝐫^{};t)`$ . This macroscopically occupied single-particle state is called Bose-Einstein condensate and its wavefunction $`\psi (𝐫,t)`$ is well described by the Gross-Pitaevskii equation

$$i\mathrm{}\frac{}{t}\psi (𝐫,t)=\left[\frac{\mathrm{}^2}{2m}^2+U(𝐫)+g|\psi (𝐫,t)|^2\right]\psi (𝐫,t),$$
(1)
where $`U(𝐫)`$ is the external potential, $`g=4\pi \mathrm{}^2a_s/m`$ is the inter-atomic strength with $`a_s`$ the s-wave scattering length, and the wavefunction is normalized to $`N`$ .

Let us consider a double-well external potential $`U(𝐫)`$. In the tunneling regime, the wavefunction of the Bose-condensate can be approximated in the following way

$$\psi (𝐫,t)=\sqrt{N_1(t)}e^{i\varphi _1(t)}\psi _1(𝐫)+\sqrt{N_2(t)}e^{i\varphi _2(t)}\psi _2(𝐫),$$
(2)
where $`N_1(t)`$ is the number of particles in the first well and $`N_2(t)`$ is the number of particles in the second well, such as $`N=N_1(t)+N_2(t)`$ . If the double-well is simmetric, it is not difficult to show that the time-dependent behavior of the condensate in the tunneling energy range can be described, with a suitable rescaling of time, by the two-mode equations

$$\dot{\zeta }=\sqrt{1\zeta ^2}\mathrm{sin}\varphi ,\dot{\varphi }=\mathrm{\Lambda }\zeta +\frac{\zeta }{\sqrt{1\zeta ^2}}\mathrm{cos}\varphi ,$$
(3)
where $`\zeta =(N_1N_2)/N`$ is the fractional population imbalance of the condensate in the two wells, $`\varphi =\varphi _1\varphi _2`$ is the relative phase (which can be initially zero), and $`\mathrm{\Lambda }=2E_I/E_T`$ with $`E_I`$ the inter-atomic energy and $`E_T`$ the tunneling energy, i.e. the kinetic+potential energy splitting between the ground state and the quasi-degenerate odd first excited state of the GP equation . Note that $`E_T=\mathrm{}\omega _0`$, where $`\omega _0`$ is the oscillation frequency of the Bose condensate between the two wells when the inter-atomic interaction is zero ($`E_I=0`$) \[1-4\].

The two-mode equations have been studied by many authors and they are the Hamilton equations of the following Hamiltonian

$$H=\frac{\mathrm{\Lambda }}{2}\zeta ^2\sqrt{1\zeta ^2}\mathrm{cos}\varphi ,$$
(4)
where $`\zeta `$ is the conjugate momentum of the generalized coordinate $`\varphi `$. Actually, it is possible to introduce  an equivalent Hamiltonian, given by

$$H^{}=\frac{1}{2}p_\zeta ^2+a\zeta ^2+b\zeta ^4,$$
(5)
where

$$a=\frac{1}{2}\sqrt{1\mathrm{\Lambda }H_0},b=\frac{\mathrm{\Lambda }^2}{8},$$
(6)
with $`H_0=\mathrm{\Lambda }\zeta (0)^2/2\sqrt{1\zeta (0)^2}\mathrm{cos}\varphi (0)`$. For the Hamiltonian $`H^{}`$ the generalized coordinate is $`\zeta `$ and its conjugate momentum is $`p_\zeta `$. The equations of motion derived by $`H^{}`$ are

$$\dot{\zeta }=p_\zeta ,\dot{p}_\zeta =2a\zeta 4b\zeta ^3.$$
(7)
The analysis of the two-mode equations (3) or (7) has shown that for $`\mathrm{\Lambda }<2`$ there is a periodic macroscopic quantum tunneling (PMQT) of the Bose-condensate with Josephson-like oscillations of frequency $`\omega _J=\omega _0\sqrt{1+\mathrm{\Lambda }}`$ . Instead, for $`\mathrm{\Lambda }>2`$, there exists a critical $`\zeta _c=2\sqrt{\mathrm{\Lambda }1}/\mathrm{\Lambda }`$ such that for $`0<\zeta (0)<<\zeta _c`$ there is PMQT of condensate with frequency $`\omega _J`$, but for $`\zeta _c<\zeta (0)1`$ there is macroscopic quantum self-trapping (MQST) of the condensate: even if the populations of the two wells are initially set in an asymmetric state ($`\zeta (0)0`$) they maintain, on the average, the original population imbalance with a very small periodic transfer of particles through the barrier with frequency $`\omega _{ST}=\omega _0\sqrt{2(\mathrm{\Lambda }H_01)}`$ .

## III An oscillating barrier in the double-well potential

In this section we analyze the effect of a periodic oscillating energy barrier in the double-well potential. The presence of a oscillating barrier can be modelled by a time-dependent $`\mathrm{\Lambda }`$. We choose the following form

$$\mathrm{\Lambda }(t)=\mathrm{\Lambda }_0(1+ϵ\mathrm{sin}(\omega t)),$$
(8)
where $`\mathrm{\Lambda }_0`$ is the static value, $`ϵ`$ is a small perturbation ($`ϵ1`$) and $`\omega `$ is the oscillation frequency of the barrier.

We numerically solve the two-mode equations (7) with (8) by using a forth-order Runge-Kutta algorithm. First we consider the case of pulsed macroscopic quantum tunneling (PMQT). We set $`\mathrm{\Lambda }_0=0.8`$, $`\zeta (0)=0.5`$ and $`\varphi (0)=0`$. In Figure 1 we plot the population imbalance $`(t,\zeta )`$ and its phase-space portrait $`(\zeta ,p_\zeta )`$ for some values of $`ϵ`$ and $`\omega `$. For $`ϵ=0`$ the population imbalance harmonically oscillates between the values $`0.5`$ and $`0.5`$ with frequency $`\omega _J=1.27`$, but for $`ϵ0`$ the motion is not fully harmonic. In particular, Figure 1 shows that under the parametric resonance condition $`\omega =2\omega _J`$ the dynamics of the population imbalance is modified: the amplitude of the oscillation is modulated and reaches the values $`1`$ and $`1`$.

Now we consider the case of macroscopic quantum self-trapping (MQST). We set $`\mathrm{\Lambda }_0=25`$, $`\zeta (0)=0.6`$ and $`\varphi (0)=0`$. Figure 2 confirms that for $`ϵ=0`$ the population imbalance does not change sign and it periodically oscillates between the values $`0.47`$ and $`0.6`$ with frequency $`\omega _{ST}=13.53`$. As shown in Figure 2, MQST is strongly affected by the parametric resonance between the frequency $`\omega _{ST}`$ of MQST oscillation and the frequency $`\omega `$ of oscillation of the energy barrier. Moreover, Figure 3 shows that at the parametric resonance condition $`\omega =2\omega _{ST}`$ and with a sufficiently large perturbation ($`ϵ=0.2`$) the system eventually escapes from the self-trapping configuration. Note that this remarkable effect has been confirmed  by using the nonpolynomial Schrödinger equation, an effective one-dimensional equation derived from the three-dimensional Gross-Pitaevskii equation .

## Conclusions

We have studied the periodic tunneling and the quantum self-trapping of a Bose-Einstein condensate in a double-well potential with an oscillating energy barrier. We have used the two-mode classical-like equations to find that the periodic macroscopic quantum tunneling can be enhanced and the macroscopic quantum self-trapping can be reduced or suppressed by changing the frequency of the oscillating barrier. In the latter case, we have proved that the system escapes from the self-trapping configuration if the the frequency of oscillation of the double-well energy barrier and the frequency of MQST oscillations of the condensate satisfy the parametric resonance condition with a sufficiently large perturbation of the energy barrier.

## References

 Smerzi, A., Fantoni, S., Giovannazzi, G. and Shenoy, S.R., 1997, Phys. Rev. Lett. 79 4950; Raghavan, S., Smerzi, A., Fantoni, S., and Shenoy, S.R., 1999, Phys. Rev. A 59 620.

 Milburn, G.J., Corney, J., Wright, E., and Walls, D.F., 1997, Phys. Rev. A 55 4318.

 Zapata, I., Sols, F., and Leggett, A.J., 1998, Phys. Rev. A 57 R28.

 Salasnich, L., Parola, A., and Reatto, L., 1999, Phys. Rev. A 60 4171; Pozzi, B., Salasnich, L., Parola, A., and Reatto, L., 2000, Eur. Phys. J. D 11 367.

 Anderson, B.P. and Kasevich, M., 1999, Science 282 1686.

 Cataliotti, F.S., Burger, S., Fort, C., Maddaloni, P., Trombettoni, A., Smerzi, A., and M. Inguscio, 2001, Science 293 843.

 Morsch, O., Muller, J.H., Cristiani, M., Ciampini, D., and Arimondo, E., 2001, Phys. Rev. Lett. 87 140402.

 Landau, L.D. and Lifsits, E.M., 1991, Mechanics, Course of Theoretical Physics, vol. 3 (Pergamon Press, Oxford); Arnold, V.I., 1990, Mathematical Methods of Classical Mechanics (Springer, Berlin).

 Courteille, P.W., Bagnato, V.S., and Yukalov, Y.I., 2001, Laser Phys. 11, 659.

 Gross, E.P., 1961, Nuovo Cimento 20 454; Pitaevskii, L.P., 1961, Zh. Eksp. Teor. Fiz. 40, 646.

 Salasnich, L., Parola, A., and Reatto, L., 2002, J. Phys. B 35, 3205.

 Salasnich, L., 2002, Laser Phys. 12 198; Salasnich, L., Parola, A., and Reatto, L., 2002, Phys. Rev. A 65 043614.


format: text

text:
# Slow relaxation in weakly open vertex-splitting rational polygons

## I INTRODUCTION

Polygonal classical billiards is an active subject of research in mathematics and physics (see for review Ref.). In view of the null Lyapunov exponent and the null Kolmogorov metric entropy the rational polygons, formed by the piece-line billiard boundary with the vertex angles are rational multiplies of $`\pi `$, are known to be nonchaotic systems. They are therefore well distinct from the Sinai billiard (SB) and the Bunimovich billiard (BB) where classical chaotic motion regimes are due to, respectively, dispersive effects caused by the circle disk and the squared boundary, and the interplay between boundary segments formed by the circle and the square. Meanwhile, the rational polygons of $`m`$ equal sides and equal vertices (hereafter, the $`m`$-gons) revealed positive Lyapunov exponents with increasing of $`m`$. Furthermore, the polygonal billiards exposed chaoticlike changes in the associated quantum-level spectra, which fluctuations are shown to be very close to the Gaussian orthogonal-ensemble-type statistics. In view of the splitting effects by the angle vertices, polygons do not satisfy the conditions of integrability

These evidences for the chaoticlike features of the nonintegrable rational polygons were recently questioned by Mantica through the orbital complexity analysis. Unlike the case of the integrable billiards, a delicate problem of the interplay between the regular (piece-line) and the irregular (vertex-angle) boundary segments in polygons cannot be solved in terms of the first-order averaged polygonal characteristics, such as the average-orbit coding length, the mean wall-collision time, or the average collision number. As shown through analysis of the orbit-wall-collision statistics, the higher order correlation effects induced by vertices play a crucial role in the intrinsic dynamics of $`m`$-gons. Earlier this was experimentally corroborated through the quantum-level spectra statistics: ”unlike the case of the irrational polygons, the long-range level correlation effects are due to the rationality of vertex angles”. In the case of the curved-by-circle infinite-$`m`$ rational polygon ($`\mathrm{}`$-gon) the late-time memory effects arise from the sliding orbits, which have no analog in the ballistic-type dynamics of its counterpart given by the circle billiard (hereafter CB), where the orbit classification is well established. This implies that the quasi-classical approach has no justification for the ”quantized” vertex-spliting $`\mathrm{}`$-gon that is geometrically equivalent to the CB. In other words, the quantum-to-classical dynamics correspondence suggested between a given $`m`$-gon and its circumscribing counterpart is violated regardless of the fact that the geometric correspondence exists and can be achieved with any precision when $`m\mathrm{}`$ (with the help of the aforementioned first-order averaged characteristics). This is in agreement with a conclusion on inapplicability in polygons of the quantum-to-classic correspondence principle. The latter was elaborated within the scope of the conventional Wentsell-Kramer-Brillouin picture that failed to establish a one-to-one correspondence between classical orbits and their quantum counterparts.

The issue of the current paper is an investigation of the vertex-splitting effects revealed in the  decay dynamics of finite-$`m`$ rational polygons. An analysis is given through the survival orbit statistics in the weakly open billiards, which boundaries permit orbits to escape through a small opening. We will see that, similarly to the case of the intrinsic dynamics in the closed $`m`$-gons, the vertex-splitting effects are dual of vertex-ordering and vertex-disordering motion effects, that is manifested through, respectively, the orderlike and chaoticlike behavior observed in the late-time relaxation. The paper is organized as follows. The recent findings for the decay dynamics in nonchaotic and chaotic billiards are given in Sec.II and discussed within concepts of the primary and the secondary relaxation channels. The weakly open rational polygons are analyzed numerically and analytically in Sec.III for the cases of small and large numbers $`m`$. Discussion and conclusions are summarized in Sec.IV.

## II DECAY DYNAMICS OF THE CHAOTIC AND NONCHAOTIC BILLIARDS

The intrinsic dynamics of the closed classical billiards is commonly considered in terms of a temporal decay of the correlation functions for certain dynamical variables (see e.g. Ref.). Pure exponential lost of amount of memory on the initial states is not a unique channel of relaxation even in the chaotic systems (see e.g. Ref.). By studying the chaotic billiards, such as the SB, that is dynamically equivalent to the correspondent Lorentz gas (LG) model and the BB with a stadium geometry, it has been recognized that a crossover from the short-time exponential to the late-time algebraic decay is due to the long-term memory effects on a free regular motion. The algebraic tail of the correlation functions seems to be vanished only in the case of the fully hyperbolic systems correspondent to such geometries as the finite-horizon SB (equivalent to the high-density LG) or the diamond. Qualitatively the same can be referred to the  decay dynamics in weakly open billiards that describes a crossover from a bounded to unbounded orbit free motion. Such a decay dynamics is initially established by uniformly distributed $`N_0`$ point particles (of unit mass and unit velocity) moving inside of the closed planar billiard table and then allowing to escape through a small opening of width $`\mathrm{\Delta }`$. A temporal behavior of the dynamic observables can be scaled to the characteristic billiard times, namely

$$\tau _c=\frac{\pi A}{P}\text{ and }\tau _e=\frac{P}{\mathrm{\Delta }}\tau _c\text{.}$$
(1)
Here the mean collision time $`\tau _c`$ and by the mean escape time $`\tau _e`$ are introduced through the accessible area $`A`$ and the perimeter $`P`$ ($`\mathrm{\Delta }`$) for a given billiard table. The late-time ($`t\tau _e`$) algebraic-type evolution of non-escaped orbits (or particles) $`N_tN_0(\tau _e/t)^\delta `$, which time derivation is related to the survival probability, is definitely characterized by the decay dynamic exponent $`\delta `$.

In the integrable nonchaotic billiards, the first observation of the algebraic decay with exponents $`\delta 1`$ was given for the case of the square billiard. Recent study of the decay dynamics in the $`4`$–gon revealed two distinct channels of the algebraic-type slow relaxation. The first one is due to the regular-orbit motion with the decay dynamic exponent $`\delta =1`$, and the second channel is originated from the ”irregular” orbits induced by the singular vertex-spliting effects, which give rise to the subdiffusion motion regime indicated by $`\delta 0.85`$. Such a kind of diverse dynamical behavior was observed through the  survival orbit spectra defined by a number for the survived orbits and simulated in the integrable CB and in the almost-integrable square billiard (see, respectively, Figs. 5,6 and 3,4 in Ref.). In both the cases an irregularlike motion is due to orbit families known as, respectively, the ”whispering-gallery” (polygonlike) and the ”bouncing-ball” orbits. Meanwhile, in the late-time decay of the CB (see inserts in Figs. 3,4 and 5,6 in Ref.) the short-time living, frequently escaped ”whispering-gallery” orbits, unlike the nonintegrable case, did not contribute to the second channel of relaxation. Thus, the solely exponent $`\delta =1`$ was observed in the CB.

In the chaotic closed and weakly open (including Hamiltonian) classical systems, presented by the BB, the infinite-horizon SB, and by the low-density LG model, the algebraic-type decay was numerically revealed by the dynamic exponents $`\delta 1`$. Similarly to the nonchaotic case, it has been repeatedly recognized that the algebraic tail is caused by the ”arbitrary long segments” observed in the evolution of stochastic orbits, or by the regularlike orbit motion due to ”sticking particles” . This implies that in both the cases this relaxation is due to a free motion of the corresponding trajectories in the infinite (with respect to the observation time) distinct corridors which are open in the relevant phase space.

The algebraic-type relaxation channel with $`\delta =1`$ ($`\alpha `$) established in the chaotic as well as in nonchaotic billiards, seems to be generic for all non-fully hyperbolic systems with smooth convex boundaries. Its independence of the billiard space dimension, its insensitiveness to details of the boundary shape, including that to a position of the small opening, and to the initial conditions, suggests that the late-time $`\alpha `$-relaxation arise in classical systems as the universal primary relaxation. The latter is a part of the two-step relaxation scenario discussed in the chaotic and nonchaotic weekly open classical systems. This universal scenario was introduced by the short-time pure exponential and by the first-power-algebraic overall orbit decays given by $`N_tN_0e^{(t/\tau _e)^\gamma }`$ with $`\gamma =1`$, and by $`N_tN_0(\tau _e/t)^\delta `$ with $`\delta =1`$, respectively. Meanwhile the real relaxation does not exclude the existence of the additional intermediate transient regimes approximated by the stretched-exponential form with $`\gamma <1`$ (earlier discussed for chaotic billiards in analytical and numerical forms), and by other algebraic-decay forms with $`\delta 1`$. The escape mechanism of the primary relaxation in the chaotic and nonchaotic billiards was described in details within a coarse-grained approximation (see, respectively, Eq.(12) in Ref. and Eq.(15) in Ref.). Unlike the case of the primary relaxation, the temporal observation conditions (observation windows) for the secondary relaxation in the chaotic billiards ($`\delta =\beta >1`$) are shown to be very sensitive to the billiard geometry, to the dimension $`d`$ of a billiard table, and to the initial conditions. This can be exemplified by the dynamic-exponent constraint $`1<\beta d`$ proposed in Ref. and observed in the chaotic BB and SB billiards.

The algebraic-type decay of correlations in the low-density chaotic LG is due to evolution of trajectories within the infinite principal and/or ”hindered” open corridors. In the corresponding disk-dispersing SB (of side $`L`$ and disk radius $`R`$) with the infinite-horizon geometry ($`R<L/2`$) the principal and the ”hindered” corridors, respectively, can be governed by the limiting disk radii $`R_\alpha =L/2`$ and $`R_\beta =\sqrt{2}L/4`$ and observed through the established algebraic-tail temporal windows. Thus, the primary relaxation was observed within the domain $`R_\beta `$ $`<R<`$ $`R_\alpha `$, when the ”hindered” corridors are closed. Under the geometric conditions $`0<R`$ $`R_\beta `$ the transient $`\beta `$-relaxation channel was activated and indicated by the decay dynamic exponents $`\beta =1.2`$ and $`\beta =1.1`$, respectively, in Refs. and . One can see that the primary and the secondary relaxation channels can be geometrically associated with Bleher’s principal and ”hindered” corridors, respectively. The latter case was additionally characterized as a superdiffusive motion regime, observed earlier in the chaotic and very recently in the polygonal closed billiards.

We see that non-fully hyperbolic billiards, on one hand, are indiscernible within decay dynamics observed through the primary universal $`\alpha `$-relaxation channel. On the other hand, the chaotic and nonchaotic, the open and closed billiards are well distinguished with respect to the secondary $`\beta `$-relaxation given by the decay dynamic exponents $`\beta >1`$ and $`\beta <1`$, respectively. In what follows we give numerical and analytical analysis of the conditions of stabilization of both the primary and secondary relaxation channels in the weakly open rational polygonal billiards.

## III ORBIT DECAY IN POLYGONS

We deal with the rational polygons of $`m`$ equal sides, denominated as $`m`$-gons, circumscribed below a circle of radius $`R`$. The mean collision time $`\tau _{cm}=(\pi R/2)\mathrm{cos}(\pi /m)`$ and the mean escape time $`\tau _{em}=(\pi R^2m/2\mathrm{\Delta })\mathrm{sin}(2\pi /m)`$ are given with the help of Eq.(1) through area $`A_m=(mR^2/2)\mathrm{sin}(2\pi /m)`$ and perimeter $`P_m=2mR\mathrm{sin}(\pi /m)`$. In the limit $`m\mathrm{}`$ one naturally arrives at the circle geometry of the $`\mathrm{}`$-gon with the mean times $`\tau _c\mathrm{}=\tau _{cR}^{(CB)}=\pi R/2`$ and $`\tau _e\mathrm{}=\tau _{eR}^{(CB)}=\pi ^2R^2/\mathrm{\Delta }`$ , both characteristic of the CB. This demonstrates how the first-order dynamic characteristics can be introduced through the aforementioned geometrical correspondence that takes place between the $`\mathrm{}`$-gon and the CB. Meanwhile the dynamical correspondence does not exist because the vertex-memory effects violate commutation between of the temporal ($`t\mathrm{}`$) and spatial ($`m\mathrm{}`$) limits.

### A Small Number of Vertices

Similarly to the closed $`m`$-gons, let us consider the case of small number of vertices, with $`m<10`$ , within the scope of the deterministic approach. This is given by generalization of the regular-orbit description introduced for the particular case of $`m=4`$ and is straightforwardly based on the fact that the wall-collision angles $`\phi `$ (counted off the normal to the boundary and preserved by elastic refections) are integrals of motion. This is true for the integrable billiards where dispersing or splitting effects are absent. A dynamic description of the regular-orbit motion can be introduced with accounting of the fact that $`m`$ (or $`m/2`$) sides of a given $`m`$-gon, with odd (or even) number of vertices, are dynamically equivalent. A description of the wall-collision statistics can be therefore reduced to the collision-angle domain $`\phi =[0,\phi _m]`$ with

$$\phi _m\text{ }=\{\begin{array}{cc}\pi /2m& \text{ for odd }m\text{,}\\ \pi /m& \text{for even }m.\end{array}$$
(2)
In turn, the $`\phi `$-family regular–orbit sets can be introduced (for details see Appendix) through the characteristic collision times, namely

$`t_{cm}(\phi )`$ $`=`$ $`{\displaystyle \frac{\pi R}{2\phi _m}}{\displaystyle \frac{\mathrm{cos}\pi /m\text{ }\mathrm{sin}\phi _m}{\mathrm{cos}(\phi \psi _m)}}\text{ , with}`$ (3)
$`\psi _m`$ $`=`$ $`\{\begin{array}{cc}0& \text{for odd }m\text{ and }m/2,\\ \pi /m& \text{for even }m/2.\end{array}`$ (6)
The collision time $`t_{cm}(\phi )`$ is related to the billiard mean collision time $`\tau _{cm}`$ through the ”mean-collision-time” equation, i.e., $`<t_{cm}(\phi )>_c_0^{\phi _m}t_{cm}(\phi )f_{0m}(\phi )𝑑\phi =`$ $`\tau _{cm}`$. The latter is equivalent to the known ”mean-free-path” equation considered in the uniformly populated two-dimensional ($`2D`$) collision subspace $`\mathrm{\Omega }_{cm}`$ (see also Eqs.(3,4) in Ref.) the $`3D`$ phase space $`\mathrm{\Omega }_m`$. The $`\phi `$-set-orbit distribution function $`f_{0m}`$ is defined here as

$$f_{0m}(\phi )=\{\begin{array}{cc}\frac{\mathrm{cos}(\phi \psi _m)}{\mathrm{sin}\phi _m}\text{}& \text{for subspace }\mathrm{\Omega }_{cm},\\ \frac{t_{cm}(\phi )}{\tau _{cm}}\frac{\mathrm{cos}(\phi \psi _m)}{\mathrm{sin}\phi _m}\frac{1}{\phi _m}\text{,}& \text{for space }\mathrm{\Omega }_m\end{array}$$
(7)
by generalization of Eqs.(6,7) in Ref..

We discuss the late-time ($`t\tau _{em}`$) survival dynamics in a given $`m`$-gon through the $`\phi `$-set-orbit decay spectra defined by numbers of the survived orbits $`N_m(t,\phi )`$ and by corresponding overall numbers $`N_{tm}=<N_m(t,\phi )>_c`$. Here a procedure of averaging over collisions, denoted by $`<\mathrm{}>_c`$, is introduced through the aforesaid ”mean-collision-time” equation. The universal relaxation channel, associated solely with the regular orbits, is given through numbers of the nonescaped orbits, namely

$$\text{ }\frac{N_m(t,\phi )}{N_{0m}}=C_m(\phi )\frac{t_{cm}(\phi )}{\tau _{cm}}f_{0m}(\phi )\frac{\tau _{em}}{t}\text{ and }\frac{N_{tm}}{N_{0m}}=D_m\frac{\tau _{em}}{t}\text{,}$$
(8)
which are the late-time solutions of the relevant decay-kinetics master equation (see Eqs. (16,17) in Ref.). As seen from Eq.(8), the fundamental characteristics $`t_{cm}(\phi ),`$ $`f_{0m}(\phi )`$ and $`\tau _{cm}`$ are common for both the decay and intrinsic dynamics. The $`\phi `$-set orbit-partial weight $`C_m(\phi )`$ and the orbit-overall weight $`D_m`$ of the corresponding algebraic tails can be established in explicit form within a certain coarse-grained scheme and directly observed. However an interesting analysis that indicates a departure of $`m`$-gons from the true integrable systems due to vertices can be given through the algebraic-tail weights without detail calculations.

By taking into account the mentioned equation $`N_{tm}=<N_m(t,\phi )>_c`$, the overall regular-orbit-decay weight $`D_{cm}=C_{cm}^{(\text{reg})}=<C_m(\phi )>_c^{(\text{reg})}`$ immediately follows from Eq.(8). The validity of this relation can be straightforwardly tested for the integrable CB where the irregular orbits do not survive in the late-time relaxation. Indeed, in this case $`C_{cR}^{(\mathrm{exp})}=0.206`$ and $`D_{cR}^{(\mathrm{exp})}=0.214`$ that follows from the observation (see Tab.2 in Ref.) of the primary relaxation of the collision space $`\mathrm{\Omega }_{cR\text{ }}`$. One can see that the algebraic-tail parametrization of the regular-orbit decay effects in the CB is self-consistent that is experimentally justified with a high precision, i.e., $`C_{cR}^{(\mathrm{exp})}=D_{cR}^{(\mathrm{exp})}=0.210\pm 2\%`$. In the case of polygons, a violation of the relations $`C_{cm}^{(\text{reg})}=C_{cm}^{(\mathrm{exp})}`$ and $`C_{cm}^{(\mathrm{exp})}=D_{cm}^{(\mathrm{exp})}`$ is expected in view of the long-living irregular-orbit motion induced by vertex-splitting effects.

We have performed numerical experiments on decay dynamics in $`m`$-gons with small number of vertices: $`m=3,4\mathrm{}8`$. The initial particles ($`N_{0m}=10^6`$) have been distributed randomly within the two distinct phase spaces described in Eq.(7) and then allowed to escape through a small opening $`\mathrm{\Delta }R`$. The condition $`\tau _{em}`$ ($`=300`$) is chosen to be common for all $`m`$ and that has been provided with the accuracy of $`\pm 5\%`$ with the help of Eq.(1). In all cases of $`m8`$ the typical algebraic decay is observed within a certain temporal windows given by, approximately, $`10^1\tau _{em}<t10^3\tau _{em}`$. The particular cases of the observed decay spectra for pentagon and heptagon are exemplified in Fig.1. In general, the overall-orbit late-time decay in $`m`$-gons with small number of vertices do not show noticeable deviations from the linear relaxation (see the left insert in Fig.1). Thus the partial weights $`C_m^{(\mathrm{exp})}(\phi )`$ are derived from the observed numbers $`N_m^{(\mathrm{exp})}(t,\phi )`$ through Eq.(8) with accounting of the estimated distribution function $`f_{0m}(\phi )`$ and the collision times $`t_{cm}(\phi )`$ given in, respectively, Eqs.(2) and (3). These equations were additionally experimentally tested (e.g. see the right insert in Fig.1).

As seen from Fig.1, the observed partial weights $`C_m^{(\mathrm{exp})}(\phi )`$ exhibit regular (small) and irregular (large) deviations from the mean magnitude $`C_{cm}^{(\mathrm{exp})}=C_{cm}^{(\text{tot})}`$ shown by a solid (horizontal) line. The latter and the regular-orbit weights $`C_{cm}^{(\text{reg})}`$ (evaluated with regardless of the large isolated peaks) are accumulated in Table.1. Through analysis of a difference between the overall and partial weights $`\mathrm{\Delta }C_{cm}=C_{cm}^{(\text{tot})}C_{cm}^{(\text{reg})}`$, one can see with the help of Table 1 that, similarly to the case of the intrinsic dynamics (see Fig.2 in Ref.), the vertex-splitting effects in the even-gons are more pronounced than those in the odd-gons. In all the $`m`$-gon cases a deviation of the total weights, treated as $`\mathrm{\Delta }D_{cm}=D_{cm}^{(tot)}D_{cm}^{(\text{reg})}`$ (with $`D_{cm}^{(\text{tot})}=D_{cm}^{(\mathrm{exp})}`$ and $`D_{cm}^{(\text{reg})}=C_{cm}^{(\text{reg})}`$ given in Table 1), exceeds experimental error ($`\pm 2\%`$) established above for the integrable SB, and we therefore infer that $`D_{cm}^{(tot)}>D_{cm}^{(\text{reg})}`$. This implies that the irregular-orbit motion is involved into the observed relaxation. On the other hand, no noticeable deviation from the primary relaxation is indicated in the observed decay-orbit dynamics when $`m8`$. Similarly to the case of the intrinsic dynamics observed in the $`m`$-gons with small number of vertices, we deduce that the regular-orbit motion dominates in the late-time decay dynamics.

### B Large Number of Vertices

The universal two-step relaxation in open billiards is shared by $`m`$-gons with arbitrary number of vertices. In Figs. 2 and 3 we analyze our numerical results for the late-time overall-orbit decay dynamics in the $`m`$-gons with large number ($`m=2^n`$ with $`n=3,\mathrm{}6`$) of vertices for the cases of relatively small and large opening widths. In general, one can observe that the decay dynamics of the rational polygons with increasing of number of sides moves away from that given by the geometrically corresponding CB: the chaotic effects, manifested by the secondary relaxation channel with $`\beta >1`$, become more pronounced with number of vertices. As seen from the observed relaxation of the initially equivalent states (given in phase spaces $`\mathrm{\Omega }_{\mathrm{}}`$ and $`\mathrm{\Omega }_R`$) is qualitatively distinct, and the open $`\mathrm{}`$-gon and the CB are not therefore dynamically equivalent. Moreover, the $`m`$-gons with $`m>8`$ do not expose the algebraic decay with $`\beta <1`$, characteristic for subdiffusive motion regimes in nonchaotic open systems observed in $`4`$-gon.

In the particular case of $`\mathrm{\Delta }=0.05R`$ shown in Fig. , the universal relaxation remains stable until $`m=64`$, but when $`mm_\alpha ^{(\mathrm{exp})}=128`$ the primary relaxation channel turns up to be closed. In other words, the regular-orbit-motion relaxation affected by vertices is transformed into the irregularlike-motion relaxation indicated by dynamic exponent $`\beta =1.2`$. Qualitatively the same follows from Fig. 3, but the upper limit for the $`\alpha `$-channel observation window shows its dependence on $`\mathrm{\Delta }`$, i.e., $`m_\alpha ^{(\mathrm{exp})}=32`$ for $`\mathrm{\Delta }=0.20R`$. Thus Bleher’s principal corridors of the regular-motion relaxations are open for $`3m<m_\alpha ^{(\mathrm{exp})}`$ (similarly to the case of the SB given by $`R_\beta `$ $`R<`$ $`R_\alpha `$ that is discussed in Sec.II) when the vertex-splitting (or disk-dispersing effects) are weak. We see that the $`\alpha `$-relaxation ( $`3m<m_\alpha ^{(\mathrm{exp})}`$) occurs as a precursor of the $`\beta `$-relaxation regime, realized for $`mm_\alpha ^{(\mathrm{exp})}`$, just as in the chaotic SB case the $`\beta `$-relaxation ($`0<R`$ $`R_\beta `$) was observed before the $`\alpha `$-relaxation ($`R_\beta `$ $`R<`$ $`R_\alpha `$). Conversely, the observed $`\alpha `$-to-$`\beta `$-relaxation crossover in the $`m`$-gons, induced by the interplay between the piece-line regular and the vertex-angle singular boundaries, is similar to that between the semi-square and semi-circle parts of the stadium boundary. The latter was deduced in the case of the BB with small ($`\mathrm{\Delta }=0.01`$) and large ($`\mathrm{\Delta }=0.25`$) opening widths, with the observed exponents $`1\beta 2`$ (see Table 1 in Ref.).

Qualitatively, the effect of closing of the principal corridors in a given $`m`$-gon can be understood by difficulties, increasing with $`m`$, to draw the long segments of free motion, which intersect the polygonal sides in the correspondent LG lattice but avoid the vertex angles (see also discussion in Appendix). By contrast, the $`\beta `$-relaxation revealed in Figs. 2, 3 is associated with stabilization of the ”irregular”-type trajectories, which are effectively deviated by vertex angles. More precisely, the observed order-to-chaoticlike crossover can be understood as a regular-to-irregular orbit transformation of the aforesaid sliding orbits (formed by $`\phi \pi /2`$ -sets with characteristic free-path times $`t_{cm}^{(reg)}(\phi )=\tau _c\mathrm{}\mathrm{cos}^1\phi `$ following from Eq.(3)), which survive in the open $`m`$-gons with $`m<m_\alpha `$, into those, renormalized substantially by vertex angles, called by vortexlike orbits (with the large, but finite mean characteristic time $`\tau _{cm}^{(irreg)}=\tau _c\mathrm{}m/\pi `$), which expected to be stable for $`mm_\alpha 1`$. Within this context, the survival conditions of the regular (sliding) and the irregular (vortexlike) orbits driven, respectively, by piece-line and vertex-angle parts of the open polygonal boundary, can be introduced as follows. On one hand, the favorable observation conditions for the $`\alpha `$-relaxation (or the $`\beta `$-relaxation) should ensure to exclude (or to include) the vertex-angle effects under the constraint $`m<m_\alpha `$ (or $`mm_\alpha `$). On the other hand, in the weakly open ($`\mathrm{\Delta }P_m`$) $`m`$-gon of a side length $`L_m=P_m/m`$, the survival conditions for the regular sliding or irregular vortexlike orbits are satisfied by geometric constraints, respectively, $`\mathrm{\Delta }L_m`$ or $`\mathrm{\Delta }L_m`$. Hence, the $`\alpha `$-$`\beta `$-relaxation crossover, observed at $`m=m_\alpha `$, is ensured by the condition $`\mathrm{\Delta }=L_m`$. With taking into account that perimeter in the polygons with large number of sides is well approximated by $`P_m2\pi R`$, one arrives at the desirable criterium

$$m_\alpha =\frac{2\pi R}{\mathrm{\Delta }}.$$
(9)
This finding provides the estimates $`m_\alpha ^{(\text{theor})}=126`$ and $`31`$ for the particular cases of $`m_\alpha ^{(\mathrm{exp})}=128`$ and $`32`$ realized in Figs. 2 and 3, respectively. We infer that unlike the case of the closed polygons, the vertex-splitting effects in the open $`m`$-gons with large number of vertices give rise to stabilization of the vortexlike-orbit motion.

## IV DISCUSSION AND CONCLUSIONS

The mild discontinues caused by vertex angles and relative lengths of the edges is the central problem of the intrinsic dynamics of the ”almost integrable” polygonal billiards commonly discussed in terms of the orbit ergodicity, mixing, entropy, coding, complexity, pseudo-integrability, spectral level statistics, and the orbit collision statistics. The problem is now addressed to the decay dynamics in the $`m`$-gons and is discussed through the orbit survival probability $`\mathrm{\Psi }_m(t)=|d(N_{tm}/N_{0m})/dt|`$ related to the number of the survived orbits $`N_{tm}`$. The decay spectra given by the $`\phi `$-set regular-orbit numbers $`N_m(t,\phi )`$, are also studied for the case of small number of vertices $`m`$.

A general approach to the decay problem based on a simple decay kinetic equation naturally arrives at the primary slow relaxation of the regular-orbit sets given by $`\mathrm{\Psi }_m^{(\alpha )}t^2`$ (see Eq.(8)). We have demonstrated that the universal $`\alpha `$-channel, attributed for both the chaotic and nonchaotic billiards, is also characteristic of nonintegrable rational polygons. The primary relaxation-motion regime originated from the piece-line parts of the polygonal table is associated with long-living sliding orbits with large collision angles $`\phi `$. In the corresponding phase space these orbits are unbounded trajectories (see Fig.4) that move without splitting at angle vertices along Bleher’s principal corridors. Following to the simplest polygonal orbit classification by Gutkin, the regular sliding-orbit sets can be presented by the ”infinite-past-to-infinite-future” trajectories. They ”never” hit vertices, preserve the initial linear momenta, and show a regular behavior in the orbit decay spectra $`N_m(t,\phi )`$ (see Fig.1). Conversely, the singular orbit sets caused by the ”infinite-past-to-vertex”, the ”vertex-to-infinite-future”, and the diagonal ”vertex-to-vertex” trajectories exhibit pronounced weights $`C_m(\phi )`$ in the orbit-decay process (shown by high peaks in Fig.1). Eventually, they do not play any significant role in the wall-collision statistics in $`m`$-gons with small number of vertices limited by, approximately, $`3m8`$, and thus the primary relaxation dominates. This corroborates by our numerical study (analyzed in Table 1and Fig.1) and, in general, is in accord with studies of the closed polygons by the orbit-wall collision statistics and by the orbit complexity.

When the number of vertices is large, the secondary relaxation with the survival probability $`\mathrm{\Psi }_m^{(\beta )}t^{\beta 1}`$ predominates over the primary relaxation (see Figs. 2,3). The established domain for the decay exponent $`1<\beta <2`$ corresponds to that known for the chaotic SB. Qualitatively, the survival probability function $`\mathrm{\Psi }_m^{(\beta )}`$ can be associated with the distribution function for trajectories trapped by the strange attractors, discussed in the theory of the open classical chaotic systems, or with the corresponding waiting-time probability function. With accounting of findings for the decay dynamics on the SB by Fendrik’s group, one can expect that the secondary relaxation is due to the singular trapped orbit sets that move freely along Bleher’s ”hidden corridors”. Similarly to the case of the chaotic billiards, the observation conditions for the secondary relaxation in rational polygons are sensitive to the initial conditions and to the geometrical constrains. Indeed, the $`\beta `$-relaxation channel turns up to be closed if the initial particle distribution is simulated in the 2D collision subspace $`\mathrm{\Omega }_{cm}`$. In the case of the 3D $`\mathrm{\Omega }_m`$ space the secondary relaxation appears to be dynamically stable under the geometrical constraint $`m>m_\alpha `$, where $`m_\alpha =2\pi R/\mathrm{\Delta }`$ is given by the $`\alpha `$-to-$`\beta `$-relaxation criterium estimated in Eq.(9). As shown, this criterium meets the favorable survival conditions for the regular-motion regime with those induced by rationality of the vortices. The latter are generated by the sliding orbits through the vertex-”ordering” effects and are associated with the vortexlike orbits. As follows from Eq.(9), the observation window of such a motion disappears in the limit $`\mathrm{\Delta }0`$, when the vortexlike orbits do not survive in the closed polygons (see Fig.2 in Ref.). Finally, we have demonstrated that the vortexlike orbits become stable in the open rational polygons and visible through the secondary slow relaxation common for the chaotic systems.

ACKNOWLEDGMENTS

The authors are grateful to Josef Klafter for drawing their interest to the escape problem in chaotic systems. Special thanks are due to Mario Jorge Dias Carneiro for illuminating discussions. The financial support of the Brazilian agency CNPq is also acknowledged.

## V APPENDIX. ORBIT-SET COLLISION TIMES

In a given $`m`$-gon a number of geometrically equivalent walls $`k`$ is bounded above by

$$q_m=\{\begin{array}{cc}m& \text{for odd }m;\\ m/2& \text{for even }m.\end{array}$$
(10)
The current collision angle $`\phi _{km}`$ with a wall $`k`$ $`(=1,2\mathrm{}q_m)`$ of a $`\phi `$-orbit with $`\phi =[0,\pi /2q_m]`$ is reduced through the relation $`\phi _{km}=\phi \mathrm{\Theta }_{km}`$ with the help of $`\mathrm{\Theta }_{km}=[\pi /2,\pi /2]`$ introduced as the lowest angle between the $`k`$-wall and axis $`x`$ , namely

$$\mathrm{\Theta }_{km}=\frac{\pi }{2q_m}\{\begin{array}{cc}q_m2k+1& \text{for odd }q_m;\\ q_m2k& \text{for even }q_m.\end{array}$$
(10)
As shown in Fig.4 for the particular case $`m=3`$, the estimates for the wall-collision times $`t_{cm}(\phi )`$ are found through summation of numbers of intersections $`n(t,\phi _{km})`$ for a trajectory, induced by a given $`\phi `$-set orbit, considered in the correspondent infinite LG lattice, namely

$$n_{cm}(t,\phi )\frac{t}{t_{cm}(\phi )}=\underset{k=1}{\overset{q_m}{}}n(t,\phi _{km})$$
(11)
The estimation procedure can be exemplified by a relation $`t\mathrm{cos}(\phi _{13})=n(t,\phi _{13})3a_3`$. The latter employes the fact that a distance between the equivalent walls is $`3a_3`$, where $`a_m=R\mathrm{cos}(\pi /m)`$ stands for the apothem in a given $`3`$-gon. This yields

$$t_{cm}(\phi )=a_mq_m\left[\underset{k=1}{\overset{q_m}{}}\mathrm{cos}(\phi \mathrm{\Theta }_{km})\right]^1$$
(12)
where $`q_m`$ and $`\mathrm{\Theta }_{km}`$ are given in Eqs.(A1) and (A2), respectively. Straightforward estimation of Eq.(A4) results in the collision times $`t_{cm}(\phi )`$ given in Eq.(3).

FIGURE CAPTIONS

.

Fig. 1. Analysis of the algebraic decay simulated in pentagon ($`5`$-gon) and heptagon ($`7`$-gon) within the collision $`\mathrm{\Omega }_{cm}`$ space. Symbols \- numerical data on the partial weights $`C_5^{(\mathrm{exp})}(\phi )`$ and $`C_7^{(\mathrm{exp})}(\phi )`$ deduced from the observed spectra of the survived orbits $`N_{t5}^{(\mathrm{exp})}(\phi )`$ and $`N_{t7}^{(\mathrm{exp})}(\phi )`$ with the help of Eq.(8) and simulated for $`\phi `$-set orbits with $`0\phi \phi _m`$ at distinct times $`t=20,30\tau _{em}`$ . Line: the overall-collision-anlle weight $`C_{cm}^{(\mathrm{exp})}`$.

Insert left: Points \- data $`N_{t5}`$ for the overall survived orbits at late times and their regular-orbit analysis with the help of Eq.(8).

Insert right: Points \- data for $`\phi `$-set collision time $`t_{c5}(\phi )`$ simulated within the reduced domain $`0\phi \pi /10`$. Line \- the same predicted in Eq.(3).

.

Fig. 2. Temporal evolution of the survived orbits in rational polygons with small opening width ($`\mathrm{\Delta }=0.05R`$) against the reduced time in log-log coordinates. Reduction is given by the help of Eq.(1) for the escape time $`\tau _e=300`$, chosen common for all cases. Points: numerical data for decay of the $`\mathrm{\Omega }_m`$ space phase simulated by $`N_0=10^6`$ particles in the $`m`$-gons (squares) and the correspondent CB (circles).

.

Fig. 3. Temporal evolution of the survived orbits in rational polygons with large opening width ($`\mathrm{\Delta }=0.20R`$) against the reduced time in log-log coordinates. Notations are the same as in Fig.2.

.

Fig.4. Estimation of the $`\phi `$-orbit collision time $`t_{cm}(\phi )`$ on the bases of Eq.(A3) for the case of $`m=3`$. The regular piece-line orbit $`a,b,c,d,e\mathrm{}`$ is represented by the infinite straight-line trajectory in the triangle LG lattice with the intersection-point sequences $`1,2,3,4\mathrm{}n(t,\phi _{km})`$. The equivalent walls $`k`$, the unreduced collision angles $`\phi _{km}`$, and the axillar angles $`\mathrm{\Theta }_{km}`$ are shown.

.

TABLE 1

.

| $`m`$ | $`C_{cm}^{(\text{tot})}`$ | $`C_{cm}^{(\text{reg})}`$ | $`D_{cm}^{(\text{tot})}`$ |
| --- | --- | --- | --- |
| 3 | 0.140 | 0.116 | 0.135 |
| 4 | 0.220 | 0.219 | 0.210 |
| 5 | 0.094 | 0.086 | 0.090 |
| 6 | 0.149 | 0.139 | 0.150 |
| 7 | 0.092 | 0.069 | 0.090 |
| 8 | 0.099 | 0.096 | 0.100 |

.

Table 1. Fitting parameters of the temporal algebraic decay of the collision space $`\mathrm{\Omega }_{cm}`$ simulated in the weakly open $`m`$-gons with $`\mathrm{\Delta }=0.05R`$. Notations: $`C_{cm}^{(\text{tot})}`$ and $`C_{cm}^{(\text{reg})}`$ correspond to the data on $`\phi `$-sets observed in the decay spectra $`C_{cm}^{(\text{exp})}(\phi )`$ (see Fig. 1) and averaged over, respectively, all collision angles and with excluding singular-orbit angles manifested by the high peaks; $`D_{cm}^{(\text{tot})}=D_{cm}^{(\text{exp})}`$\- the overall-set weights of the algebraic tail given in Eq.(8) and derived within the primary relaxation window (see the left insert in Fig.1).

.


text:
# The Gaiasphere and the limits of knowledge

## 1. Introduction

Eddington once remarked that ‘the contemplation in natural science of a wider domain than the actual leads to a far better understanding of the actual.’ Our intuition is that any dynamical phase early on in the history of the Galaxy which is dominated by relaxation or dissipation probably removes more information than it retains. This remains largely true for controlled experiments within terrestrial laboratories for the reason that it is exceedingly difficult to track fluid particles, for example, in order to unravel the process under study. However, at least in principle, the situation may be more tractable in astrophysics in that individual stars carry information about the birth site at the time of their formation. Nature offers us a vast untapped reservoir of information in the detailed chemical abundances stored within stellar atmospheres.

The fossil evidence of galaxy formation is not restricted to chemical abundances alone. In our review, “The New Galaxy – Signatures of its Formation,” we stress the importance of new observations across a wide range of parameter space (Freeman & Bland-Hawthorn 2002). We are coming into a new era of galactic investigation, in which one can study the fossil remnants of the early days of the Galaxy in a broader and more focussed way, not only in the halo but throughout the major luminous components of the Galaxy. The goal of these studies is to reconstruct as much as possible of the early galactic history.

What do we mean by the reconstruction of early galactic history? We seek a detailed physical understanding of the sequence of events which led to the Milky Way. Ideally, we would want to tag (i.e. associate) components of the Galaxy to elements of the protocloud – the baryon reservoir which fueled the stars in the Galaxy.

Our approach is anchored to observations of the Galaxy, interpreted within the broad scope of the CDM hierarchy. Many of the observables in the Galaxy relate to events which occurred long ago, at high redshift. Fig. 1 shows the relationship between look-back time and redshift in the context of the $`\mathrm{\Lambda }`$CDM model: the redshift range (z $`\stackrel{<}{}`$ 6) of discrete sources in contemporary observational cosmology matches closely the known ages of the oldest components in the Galaxy. The Galaxy (near-field cosmology) provides a link to the distant universe (far-field cosmology).

Here we address the prospect of chemically ‘tagging’ stars to common sites of formation, in order to reconstruct star clusters which have long since dispersed. High resolution spectroscopy at high signal-to-noise ratio of many stellar types reveals an extraordinarily complex pattern of spectral lines. The spectral lines carry key information on element abundances that make up the stellar atmosphere. Some of these elements cannot arise through normal stellar evolution, and therefore must reflect conditions in the progenitor cloud at the time of its formation. A cornerstone of near-field cosmology must be to explore how much of the Galaxy’s past can be reconstructed from the chemical signatures in stellar atmospheres.

## 2. Chemical signatures

A major goal of near-field cosmology is to tag individual stars with elements of the protocloud. For many halo stars, and some outer bulge stars, this may be possible with phase space information provided by Gaia. But for much of the bulge and the disk, secular processes cause the populations to become relaxed (i.e. the integrals of motion are partially randomized). In order to have any chance of unravelling disk formation, we must explore chemical signatures in the stellar spectrum. Ideally, we would like to tag a large sample of representative stars with a precise time and a precise site of formation.

Over the last four decades, evidence has gradually accumulated (Fig. 2) for a large dispersion in metal abundances \[Xi/Fe\] (particularly n-capture elements) in low metallicity stars relative to solar abundances (Wallerstein et al. 1963; Pagel 1965; Spite & Spite 1978; Truran 1981; Luck & Bond 1985; Clayton 1988; Gilroy et al. 1988; McWilliam et al. 1995; Norris et al. 1996; Burris et al. 2000). Elements like Sr, Ba and Eu show a 300-fold dispersion, although \[$`\alpha `$/Fe\] dispersions are typically an order of magnitude smaller.

In their celebrated paper, Burbidge et al. (1957 $``$ B2FH) demonstrated the likely sites for the synthesis of slow (s) and rapid (r) n-capture elements. The s-process elements (e.g. Sr, Zr, Ba, Ce, La, Pb) are thought to arise from the He-burning phase of intermediate to low mass (AGB) stars (M $`<`$ 10M), although at the lowest metallicities, trace amounts are likely to arise from high mass stars (Burris et al. 2000; Rauscher et al. 2001).

In contrast to the s-process elements, the r-process elements (e.g. Sm, Eu, Gd, Tb, Dy, Ho) cannot be formed during quiescent stellar evolution. While some doubts remain, the most likely site for the r-process appears to be SN II, as originally suggested by B2FH (see also Wallerstein et al. 1997). Therefore, r-process elements measured from stellar atmospheres reflect conditions in the progenitor cloud. In support of Gilroy et al. (1988), McWilliam et al. (1995) state that ‘the very large scatter means that n-capture element abundances in ultra-metal poor stars are products of one or very few prior nucleosynthesis events that occurred in the very early, poorly mixed galactic halo’, a theme that has been developed by many authors (e.g. Audouze & Silk 1995; Shigeyama & Tsujimoto 1998; Argast et al. 2000; Tsujimoto et al. 2000).

Supernova models produce different yields as a function of progenitor mass, progenitor metallicity, mass cut (what gets ejected compared to what falls back towards the compact central object), and detonation details. The $`\alpha `$ elements are mainly produced in the hydrostatic burning phase within the pre-supernova star. Thus $`\alpha `$ yields are not dependent on the mass cut or details of the fallback/explosion mechanism which leads to a much smaller dispersion at low metallicity.

There is no known age-metallicity relation that operates over a useful dynamic range in age and/or metallicity. (This effect is only seen in a small subset of hot metal-rich stars). Such a relation would require the metals to be well mixed over large volumes of the ISM. For the forseeable future, it seems that only a small fraction of stars can be dated directly. Arguably, the most promising stellar chronometers make use of radioactive dating or asteroseismology. Nucleo-cosmochronology is not yet a precise science but has potential for future development (Sneden et al. 2001$`a`$; Goriely & Arnould 2001).

## 3. Reconstructing ancient star groups

We now conjecture that the heavy element metallicity dispersion may provide a way forward for tagging groups of stars to common sites of formation. With sufficiently detailed spectral line information, it is feasible that the ‘chemical tagging’ will allow temporal sequencing of a large fraction of stars in a manner analogous to building a family tree through DNA sequencing.

Most stars are born within rich clusters of many hundreds to many thousands of stars (Clarke et al. 2000; Carpenter 2000). McKee & Tan (2002) propose that high-mass stars form in the cores of strongly self-gravitating and turbulent gas clouds. The low mass stars form within the cloud outside the core, presumably at about the same time or shortly after the high mass stars have formed. The precise sequence of events which give rise to a high mass star is a topic of great interest and heated debate in contemporary astrophysics (e.g. Stahler et al. 2000).

A necessary condition for ‘chemical tagging’ is that the progenitor cloud is uniformly mixed in key chemical elements before the first stars are formed. Another possibility is that a few high mass stars form shortly after the cloud assembles, and enrich the cloud fairly uniformly. Both scenarios would help to ensure that long-lived stars have identical abundances in certain key elements before the onset of low-mass star formation.

For either statement to be true, an important requirement is that open clusters of any age have essentially zero dispersion in some key metals with respect to Fe. There has been very little work on heavy element abundances in open clusters. The target clusters must have reliable astrometry so as to minimize ‘pollution’ from stars not associated with the cluster (Quillen 2002).

If our requirement is found not to be true, then either the progenitor clouds are not well mixed or high mass stars are formed after most low mass stars. A more fundamental consequence is that a direct unravelling of the disk into its constituent star groups would be impossible, in other words, the epoch of dissipation cannot be unravelled after all.

Consider the (extraordinary) possibility that we could put many coeval star groups back together over the entire age of the Galaxy. This would provide an accurate age for the star groups either through the color-magnitude diagram, or through association with those stars within each group that have \[n-capture/Fe\] $``$ 0, and can therefore be radioactively dated. This would provide key information on the chemical evolution history for each of the main components of the Galaxy.

But what about the formation site? The kinematic signatures will identify which component of the Galaxy the reconstructed star group belongs, but not specifically where in the Galactic component (e.g. radius) the star group came into existence. For stars in the thin disk and bulge, the stellar kinematics will have been much affected by the bar and spiral waves; it will no longer be possible to estimate their birthplace from their kinematics. Our expectation is that the derived family tree will severely restrict the possible scenarios involved in the dissipation process. In this respect, a sufficiently detailed model may be able to locate each star group within the simulated time sequence.

Freeman & Bland-Hawthorn (2002) argue that, in addition to open clusters, the thick disk is an extremely important fossil of the processes behind disk formation. The thick disk is thought to be a snap-frozen relic of the early disk, heated vertically by the the infall of an intermediate mass satellite. Chemical tagging of stars that make up the thick disk would provide clues on the formation of the first star clusters in the early disk.

## 4. Chemical abundance space

An intriguing prospect is that reconstructed star clusters can be placed into an evolutionary sequence, i.e. a family tree, based on their chemical signatures. Let us suppose that a star cluster has accurate chemical abundances determined for a large number $`n`$ of elements (including isotopes). This gives it a unique location in an $`n`$-dimensional space compared to $`m`$ other star clusters within that space. We write the chemical abundance space as $`𝒞`$(Fe/H, X1/Fe, X2/Fe, …) where X1, X2 … are the independent chemical elements that define the space (i.e. elements whose abundances are not rigidly coupled to other elements).

The size of $`n`$ is unlikely to exceed about 50 for the foreseeable future. Hill et al. (2002) present exquisite data for the metal-poor star CS 31082-001, where abundance estimates are obtained for a total of 44 elements, almost half the entire periodic table (see also Cayrel et al. 2001). In Fig.3, we show what is now possible for another metal-poor star, CS 22892-052 (Sneden et al. 2001$`a`$). The $`\alpha `$ elements and r-process elements, and maybe a few canonical s-process elements at low \[Fe/H\], provide information on the cloud abundances prior to star formation, although combinations of these are likely to be coupled (Heger & Woosley 2001; Sneden et al. 2001$`a`$). There are 24 r-process elements that have been clearly identified in stellar spectra (Wallerstein et al. 1997).

The size of $`m`$ is likely to be exceedingly large for the thin disk where most of the baryons reside. For a rough estimate, we take the age of the disk to be 10 Ga. If there is a unique SN II enrichment event every 100 years, we expect of order 108 formation sites. Typically, a SN II event sweeps up a constant mass of $`5\times 10^4`$M (Ryan et al. 1996; Shigeyama & Tsujimoto 1998). Simple chemical evolution models indicate that this must be of the right order to explain the metallicity dispersion at low \[Fe/H\] (Argast et al. 2000). Roughly speaking, there have been 103 generations of clouds since the disk formed, with about 105 clouds in each star-forming generation, such that cloud formation and dispersal cycles on a 107 yr timescale (Elmegreen et al. 2000).

Whereas the total number of star clusters over the lifetime of the thin disk is very large, the size of $`m`$ for the stellar halo (Harding et al. 2001), and maybe the thick disk (Kroupa 2002), is likely to be significantly smaller. Our primary interest is the oldest star clusters. Reconstructing star clusters within the thick disk is a particularly interesting prospect since the disk is likely to have formed within 1$``$1.5 Ga of the main epoch of baryon dissipation (Prochaska et al. 2000).

The task of establishing up to 108 unique chemical signatures may appear to be a hopeless proposition with current technology. But it is worth noting that more than 60 of the chemical elements ($`Z>30`$) arise from n-capture processes. Let us suppose that half of these are detectable for a given star. We would only need to be able to measure two distinct abundances for each of these elements in order to achieve 109 independent cells in $`𝒞`$-space. If many of the element abundances are found to be rigidly coupled, of course the parameter space would be much smaller.

It may not be necessary to measure as many as 30 elements if some can be found which are highly decoupled and exhibit large relative dispersions from star to star. Burris et al. (2000) demonstrate one such element pair, i.e. \[Ba/Fe\] and \[Sr/Fe\]. It is difficult at this stage to suggest which elements are most suited to chemical tagging. In part, this depends on the precise details and mechanism of formation of the n-capture elements at low \[Fe/H\].

The element abundances \[Xi/Fe\] show three main peaks at Z $``$ 26, Z $``$ 52, and Z $``$ 78; the last two peaks are evident in Fig. 3. There have been suggestions that the r-process gives rise to random abundance patterns (e.g. Goriely & Arnould 1996) although this is not supported by new observations of a few metal poor stars. Heavy r-process elements around the second peak compared to the Sun appear to show a universal abundance pattern (Sneden et al. 2000; Cayrel et al. 2001; Hill et al. 2002). However, Hill et al. find that the third peak and actinide elements (Z $``$90) are decoupled from elements in the second peak. We suspect that there may be a substantial number of suitable elements ($``$10) which could define a sufficiently large parameter space.

Our ability to detect structure in $`𝒞`$-space depends on how precisely we can measure abundance differences between stars. It may be possible to construct a large database of differential abundances from echelle spectra, with a precision of 0.05 dex or better; differential abundances are preferred here to reduce the effects of systematic error.

## 5. Chemical trajectories

Our simple picture assumes that a cloud forms with a unique chemical signature, or that shortly after the cloud collapses, one or two massive SN II enrich the cloud with unique yields which add to the existing chemical signature. The low-mass population forms with this unique chemical signature. If the star-formation efficiency is high ($`\stackrel{>}{}30`$%), the star group stays bound although the remaining gas is blown away. If the star-formation efficiency is low ($`\stackrel{<}{}10`$%), the star cluster disperses along with the gas. In a closed box model, the dispersed gas reforms a cloud at a later stage.

In the closed box model, each successive generation of supernovae produce stellar populations with progressive enrichments. These will lie along a trajectory in $`𝒞`$-space which can be identified in principle using minimum spanning tree methods (Sedgewick 1992). The overall distribution of the trajectories will be affected by fundamental processes like the star formation efficiency, the star formation timescale, the mixing efficiency, the mixing timescale, and the satellite galaxy infall rate.

As we approach solar levels of metallicity in \[Fe/H\], the vast number of trajectories will converge. By \[Fe/H\] $``$ -2.5, AGB stars will have substantially raised the s-process element abundances; by \[Fe/H\] $``$ -1, Type Ia supernovae will have raised the Fe-group abundances. Star clusters that appear to originate at the same location in this $`𝒞`$-space may simply reflect a common formation site, i.e. the resolution limit we can expect to achieve in configuration space. The ability to identify common formation sites rests on accurate differential abundance analyses (Edvardsson et al. 1993; Prochaska et al. 2000).

Even with a well established family tree based on chemical trajectories in the chemical $`𝒞`$-space, this information may not give a clear indication of the original location within the protocloud or Galactic component. This will come in the future from realistic baryon dissipation models. Forward evolution of any proposed model must be able to produce the chemical tree.

However, the $`𝒞`$-space will provide a vast amount of information on chemical evolution history. It should be possible to detect the evolution of the cluster mass function with cosmic time (Kroupa 2002), the epoch of a starburst phase and/or associated mass ejection of metals to the halo (Renzini 2001), and/or satellite infall (Noguchi 1998).

As we go back in time to the formation of the disk, we approach the chemical state laid down by population III stars. The lack of stars below \[Fe/H\] $``$ -5 suggests that the protocloud was initially enriched by the first generation of stars (Argast et al. 2000). However, the apparent absence of any remnants of population III remains a puzzle: its stars may have had a top-heavy initial mass function, or have dispersed into the intra-group medium of the Local Group. If one could unravel the abundances of heavy elements at the time of disk formation, this would greatly improve the precision of nucleo-cosmochronology.

## 6. Candidates for chemical tagging

Chemical tagging will not be possible for all stars. In hot stars, our ability to measure abundances is reduced by the stellar rotation and lack of transitions for many ions in the optical. The ideal candidates are the evolved FGK stars that are numerous (10%) and intrinsically bright. These can be observed at echelle resolutions ($`R>30,000`$) over the full Gaiasphere. Moreover, giants have deep, low density atmospheres that produce strong low-ionization absorption lines compared to higher gravity atmospheres. Even in the presence of significant line blending, with sufficient signal, it should be possible to derive abundance information by comparing the fine structure information with accurate stellar synthesis models. Detailed abundances of large numbers of F and G subgiants would be particularly useful, if it becomes possible to make such studies, because direct relative ages can be derived for these stars from their observed luminosities.

It is not clear at what \[Fe/H\] the r-process elements become swamped by the ubiquitous Fe-group and s-process elements. At a resolution of $`R10^5`$, many r-process elements can be seen in the solar spectrum, although the signal-to-noise ratio of about 1000 is needed, and even then the spectral lines are often badly blended (Kurucz 1991; 1995). Travaglio et al. (1999) suggest that the s-process does not become significant until \[Fe/H\] $``$ -1 because of the need for pre-existing seed nuclei (Spite & Spite 1978; Truran 1981), although Pagel & Tautvaisiene (1997) argue for some s-process production at \[Fe/H\] $``$ -2.5. Prochaska et al. (2000) detected Ba, Y and Eu in a snapshot survey of thick disk G dwarfs in the solar neighbourhood with -1.1 $`\stackrel{<}{}`$ \[Fe/H\] $`\stackrel{<}{}`$ -0.5. This survey only managed to detect a few transitions in each element although their spectral coverage was redward of 440nm with SNR $``$ 100 per pixel at $`R50,000`$. Longer exposures with $`R10^5`$ and spectral coverage down to 300nm would have detected more heavy elements.

## 7. Future progress

In our view, observations of nucleosynthetic signatures of metal-poor stars provide a cornerstone of near-field cosmology. Success in this arena requires major progress across a wide front, including better atomic parameters (Truran et al. 2001), improved supernova models, better stellar synthesis codes and more realistic galaxy formation models. There are no stellar evolutionary models that lead to a self-consistent detonation and deflagration in a core-collapse supernova event or, for that matter, detonation in a thermonuclear explosive event. Realistic chemical production at the onset of the supernova stage requires a proper accounting of a large number of isotope networks (400$``$2500) that cannot be adequately simulated yet. Modern computers have only recently conquered relatively simple $`\alpha `$ networks involving 13 isotopes. The inexorable march of computer power will greatly assist here.

There is also a key experimental front both in terms of laboratory simulations of nucleosynthesis, and the need for major developments in astronomical instrumentation. Many authors (e.g. Sneden et al. 2001$`b`$) have stressed the importance of greatly improving the accuracy of transition probabilities and reaction rates for both heavy and light ion interactions. This will be possible with the new generation of high-intensity accelerators and radioactive-beam instruments (Käppeler et al. 1998; q.v. Manuel 2000).

Progress on all fronts will require iteration between the different strands. Already, relative r-process and $`\alpha `$ element abundances for metal-poor stars have begun to constrain the yields for different stellar masses and associated mass cuts of progenitor supernovae (Mathews et al. 1992; Travaglio et al. 1998; Ishimaru & Wanajo 2000).

Detailed high resolution abundance studies of large samples of galactic stars will be crucial for the future of fossil astronomy. But we stress that in order to access a representative sample of the Gaiasphere, this will require a new generation of ground-based instruments, in particular, a multi-object echelle spectrograph with good blue response on a large aperture telescope. There is a real need for a high-resolution spectrograph which can reach hundreds or even thousands stars in a square degree or more. The Gemini Wide Field proposal currently under discussion provides an opportunity for this kind of instrument (S. Barden, personal communication). The instrument will be expensive and technically challenging, but we believe this must be tackled if we are to ever unravel the formation of the Galaxy.

Could some of the residual inhomogeneities from prehistory have escaped the dissipative process at an early stage? We may not know the answer to this question with absolute certainty for many years. But it is an intriguing thought that one day we may be able to identify hundreds or even thousands of stars throughout the Gaiasphere that were born within the same cloud as the Sun.

Acknowledgment. I am deeply indebted to Ken Freeman for his inspiration and dialogue over many years. This contribution is in honour of him, and in fact arose out of our work on the 2002 Annual Reviews article.

## References

Argast D, Samland M, Gerhard OE, Thielemann F.-K. 2000. Astron.Astrophys. 356:873–87

Audouze J, Silk J. 1995. Ap.J. 451:L49–52

Burbidge EM, Burbidge GR, Fowler WA, Hoyle F. 1957. Rev.Mod.Phys. 29:547–650

Burris DL, Pilachowski CA, Armandroff TE, Sneden C, Cowan JJ, Roe H. 2000, Ap.J. 544:302–319

Carpenter JM. 2000. Astron.J. 120:3139–61

Cayrel R, Hill V, Beers TC, Barbuy B, Spite M et al. 2001. Nature 409:691–92

Clayton D. 1988. MNRAS 234:1–36

Edvardsson B, Andersen J, Gustafsson B, Lambert DL, Nissen PE, Tomkin J. 1993. Astron.Astrophys. 275:101–52

Elmegreen BG, Efremov Y, Pudritz RE, Zinnecker H. 2000. In Proto stars and Planets IV, eds V Mannings, AP Boss, SS Russell, pp 179–215 (Arizona : Univ of Arizona Press)

Freeman KC, Bland-Hawthorn J. 2002. Ann. Rev. Astron. Astrophys., 40:487–537

Gilroy K, Sneden C, Pilachowski CA, Cowan JJ. 1988. Ap.J. 327:298–320

Goriely S, Arnould M. 1996. Astron.Astrophys. 312:327–337

Goriely S, Arnould M. 2001. Astron.Astrophys. 379:1113–22

Harding P, Morrison HL, Olszewski EW, Arabadjis J, Mateo M, et al. 20 01. Astron.J. 122: 1397–1419

Hill, V et al. 2002. astro-ph/0203462

Ishimaru Y, Wanajo S. 2000 In The First Stars ed A Weiss, TG Abel, V Hill, p 189– (Berlin: Springer)

Käppeler F, Thielemann FK & Wiescher M. 1998. Ann. Rev. Nucl. Part. Sci. 48:175–251

Kroupa P. 2002. MNRAS 330:707–18

Kurucz RL. 1991. In The Solar Interior and Atmosphere eds. AN Cox, WC Livingston, M Matthews. p 663 (Tucson: U. Arizona Press)

Kurucz RL. 1995. In Laboratory and Astronomical High Resolution Spectra eds. AJ Sauval, R Blomme, N Grevesse p 17 (San Francisco: Astron.Soc.Pacific)

Lineweaver, CH. 1999. Science 284:1503–1507

Luck RE, Bond HE. 1985. Ap.J. 292:559–77

Manuel O. 2000. Origins of Elements in the Solar System, (Kluwer: NewYork)

Mathews GJ, Bazan G, Cowan JJ. 1992. Ap.J. 391:719–735

McKee CF, Tan JC. 2002. Nature 416:59–61

McWilliam A, Preston GW, Sneden C, Searle L. 1995. Astron.J. 109:2736–2756

Noguchi M. 1998. Nature 392:253–55

Norris JE, Ryan SG, Beers TC. 1996. Ap.J.Suppl. 107:391–421

Pagel BEJ. 1965. Royal Obs Bull 104:127–51

Pagel BEJ, Tautvaisiene G. 1997. MNRAS 288:108–116

Prochaska JX, Naumov SO, Carney BW, McWilliam A, Wolfe AM. 2000. Astron.J. 120:2513–2549

Quillen AC. 2002. astro-ph/0202253

Rauscher, T. et al. 2001. astro-ph/0112478

Renzini A. 2001. In Chemical Enrichment of Intracluster and Intergalactic Medium ed F Matteucci and R Fusco-Femiano, ASP. Conf. Ser. (San Francisco: ASP)

Ryan SG, Norris JE, Beers TC. 1996. Ap.J. 471:254–78

Sedgewick R. 1992. Algorithms in C++ (Addison–Wesley)

Shigeyama T, Tsujimoto T. 1998, Ap.J. , 507, L135–39

Sneden C, Cowan JJ, Ivans II, Fuller GM, Burles S, et al.. 2000. Ap.J. 533:L139–L142

Sneden C, Cowan JJ, Beers TC, Truran JW, Lawler JE, Fuller G. 2001$`a`$. In Astrophysical Ages and Time Scales, eds T von Hippel, C Simpson & N Manset, ASP Conf. Ser. 245:235–243 (San Francisco: Astron. Soc. Pac.)

Sneden C, Lawler JE, Cowan JJ. 2001$`b`$. astro-ph/0109194

Spite M, Spite F. 1978. Astron.Astrophys. 67:23–31

Stahler SW, Palla F, Ho PTP. 2000. In Protostars and Planets IV , eds V Mannings, AP Boss, SS Russell, pp 327–351 (Arizona: Univ of Arizona Pre ss)

Travaglio C, Galli D, Gallino R, Busso M, Ferrini F, et al. 1999. Ap.J. 521:691–702

Travaglio C, Gallino R, Zinner E, Amari S, Woosley S. 1998. Meteoritics & Planetary Science 33:A155-156

Tsujimoto T, Shigeyama T, Yoshii Y. 2000. Ap.J. 531:L33–36

Truran JW. 1981. Astron.Astrophys. 97:391–93

Truran JW, Burles S, Cowan J, Sneden C. 2001. In Astrophysical Ages and Time Scales, eds T von Hippel, C Simpson & N Manset, ASP Conf. Ser. 245:226–34 (San Francisco: Astron. Soc. Pac.)

Wallerstein G, Greenstein JL, Parker R, Helfer HL, Aller LH. 1963. Ap.J. 137:280–300

Wallerstein G, Iben I, Parker P, Boesgaard AM, Hale GM, et al. 1997. Rev. Mod. Phys. 69:995–1084


text:
# Numerical Results for Ground States of Mean-Field Spin Glasses at low Connectivities

## Abstract

An extensive list of results for the ground state properties of spin glasses on random graphs is presented. These results provide a timely benchmark for currently developing theoretical techniques based on replica symmetry breaking that are being tested on mean-field models at low connectivity. Comparison with existing replica results for such models verifies the strength of those techniques. Yet, we find that spin glasses on fixed-connectivity graphs (Bethe lattices) exhibit a richer phenomenology than has been anticipated by theory. Our data prove to be sufficiently accurate to speculate about some exact results.

A theoretical understanding of the intricate dynamics of disordered systems has been a major goal of statistical physics at least since the introduction of the Edwards-Anderson spin glass model EA . Already the study of the equilibrium at low temperature, a state real disordered materials rarely achieve FH , reveals a stunning range of new phenomena, even in the simplest models such as the Sherrington and Kirkpatrick model (SK) where all spins are mutually connected SK ; TAP ; ParisiRSB ; MPV . As an intermediate step in extending the mean-field techniques toward finite-dimensional models, spin glasses on random graphs are an area of active research. Those systems have been of interest from early on because they combine infinite-range connections (like SK) with a finite, decidedly low connectivity. But those earlier studies have focused either on temperatures at the glass transition VB , on purely replica symmetric (RS) solutions KaSo , or on perturbative approaches in the (SK-)limit of large connectivity PT ; DG ; LG .

Simultaneously, the formal similarity between spin-glass Hamiltonians and the objective function of combinatorial optimization problems has been realized SA  and exploited to make RS predictions, for instance, for the bipartitioning problem on random graphs MPGB ; WS ; Banavar . This connection, with the discovery of phase transitions in combinatorial optimization problems HH  and the application of replica techniques to their study Monasson , has recently rejuvenated interest in spin glasses on random graphs. But to obtain quantitatively valuable predictions for NP-hard problems required the application of replica symmetry breaking (RSB) ParisiRSB  to those problems at finite connectivities and low temperatures which was accomplished recently MP1 ; Franz ; MPZ . Finally, these RSB methods are now being applied to spin glasses on random graphs, producing quantitatively valuable results MP1 ; MP2  at accuracies below 0.1%. At this level of accuracy, a comparison between theoretical and simulation results becomes valuable at least in two respects: Convergence of the numerical with the RSB result can verify the assumptions underlying RSB as well as the quality of the numerical method used to approximate an NP-hard problem.

In this letter we apply the extremal optimization (EO) heuristic BoPe1 ; eo\_prl  to investigate the ground state properties of spin glasses on random graphs. With this method we have sampled system sizes up to $`n=4096`$ on low-connectivity graphs. We have obtained high-accuracy results for the ground-state energies of spin glasses on ordinary random graphs (ORG) with fluctuating connectivities, and for Bethe-lattice graphs (BL) with fixed connectivities. On a smaller sample of BL, we have also obtained results for the entropy of such graphs. The energies are in excellent agreement with RSB predictions for low-connectivity BL MP2 . Both, the energies and entropies reveal a sensitivity to the even-oddness of the BL, which may explain inconsistencies with results in the SK limit PT . No such inconsistency arises for ORG and the numerical extrapolation is in good agreement with analytic results for the large-connectivity limit in RSB PT . The BL energies seem to fall (within $`0.2\%`$) on a simple line ranging from the 2-connected graphs to the SK-limit, but only for even integer connectivities, without any obvious interpolation. The BL entropies decrease linearly with the inverse connectivity for odd connectivities, but are already consistent with zero (within accuracy) at small even connectivity. More details of the numerical procedure is given elsewhere eo\_bethe ; eo\_rg .

Of the two types of random graphs are considered for this study, the BL are regular random graphs Bollobas . These graphs consist of $`n`$ vertices where each vertex possesses a fixed number $`k+1`$ of bonds MPGB ; MP1 ; MP2  with randomly selected other vertices. Alternatively, ORG are obtained by randomly connecting any pair of vertices with a specified probability $`p=c/(n1)`$, leading to a graph of average connectivity $`c`$ but where the connectivities of individual vertices are Poissonian distributed Bollobas . Note that each vertex’ connectivity, and thus $`k+1`$, is inherently discrete, while $`c`$ can take on any real value.

Once a graph of connectivity $`c`$ is generated, randomly chosen quenched couplings $`J_{i,j}\{1,+1\}`$ are assigned to existing bonds between neighboring vertices $`i`$ and $`j`$. Each vertex $`i`$ is occupied by a spin variable $`x_i\{1,+1\}`$. The energy of the system is defined as the difference in number between violated bonds and satisfied bonds, $`H=_{\{bonds\}}J_{i,j}x_ix_j`$, and we will focus on the energy and entropy per spin, resp.,

$`e_c={\displaystyle \frac{1}{n}}H,s_c={\displaystyle \frac{1}{n}}\mathrm{ln}\mathrm{\Omega },`$ (1)
where $`\mathrm{\Omega }`$ is the degeneracy of the configurations exhibiting the ground state energy.

For our numerical procedure we used the following implementation of EO BoPe1 ; eo\_prl : For a given spin configuration on a graph, assign to each spin $`x_i`$ a “fitness” $`\lambda _i=\mathrm{\#}violatedbonds=0,1,2,\mathrm{},c_i`$, so that $`e_c=_i\lambda _i/(2n)`$ is satisfied. Here, $`c_i`$ is the integer connectivity of vertex $`i`$, and $`c_ik+1`$ for every vertex in BL. If $`c_{\mathrm{max}}=\mathrm{max}_ic_i`$, each spin falls into one of only $`c_{\mathrm{max}}+1`$ possible states. Say, currently there are $`n_{c_{\mathrm{max}}}`$ spins with the worst fitness, $`\lambda =c_{\mathrm{max}}`$, $`n_{c_{\mathrm{max}}1}`$ with $`\lambda =c_{\mathrm{max}}+1`$, and so on up to $`n_0`$ spins with the best fitness $`\lambda =0`$, where $`_jn_j=n`$. Now draw a “rank” $`l`$ according to the distribution $`P(l)l^\tau `$. Then, determine $`0jc_{\mathrm{max}}`$ such that $`_{i=j+1}^{c_{\mathrm{max}}}n_i<l_{i=j}^{c_{\mathrm{max}}}n_i`$. Finally, select any one of the $`n_j`$ spins in state $`j`$ and reverse its orientation unconditionally. As a result, it and its neighboring spins change their fitness. After all the effected $`\lambda `$’s and $`n`$’s are reevaluated, a new spin is chosen for an update.

The arguments given in eo\_jam  and a few experiments indicate that $`\tau =1.3`$ is a satisfactory choice to find ground states efficiently on either type of graph. Our implementation restarts for each instance at least $`r_{\mathrm{max}}=4`$ times with new random initial spin assignments, executing $`0.1n^3`$ updates per run. If a new, lower-than-previous energy state is encountered in run $`r`$, we adjust $`r_{\mathrm{max}}=2+2r`$ for that instance so that EO runs at least twice as many restarts as were necessary to find the lowest state in the first place. Especially for small $`n`$, $`r_{\mathrm{max}}`$ hardly ever exceeds 4; for larger $`n`$ a few graphs require up to 30 restarts before termination. Since EO perpetually explores new configurations it is well suited to explore also the degeneracy $`\mathrm{\Omega }`$ of low-energy states. In these runs, we used a similar approach to the above, except for setting $`r_{\mathrm{max}}=8+2r`$, where $`r`$ is the latest run in which a new configuration of the lowest energy was located.

We have simulated spin glasses on BL with this algorithm for $`k+1`$ between 3 and 26, and graph sizes $`n=2^l`$ for $`l=5,6,\mathrm{},10`$ to obtain results for ground state energies eo\_bethe . In particular, for $`k+1=3`$ we have used the methods described in Ref. eo\_rg  to reach system sizes of $`n=4096`$. In a separate simulation, using $`\tau =1.4`$, we have explored BL of size $`n[16\mathrm{}256]`$ to determine their entropy. We have used the same algorithm, preceded by a graph reduction procedure eo\_rg , to study ORG ranging from $`n2^{15}`$ for $`c=2`$ to $`n2^9`$ at $`c=25`$. Amazingly, as is shown in Refs. eo\_bethe ; eo\_rg , in all these cases our data can be extrapolated for $`n\mathrm{}`$ via

$`e_c(n)e_c+{\displaystyle \frac{A}{n^{2/3}}}(n\mathrm{}).`$ (2)
Deviations from these scaling corrections are generally small eo\_bethe  and we assume Eq. (2) to be exact here. There does not appear to be a theoretical justification for Eq. (2), but in Fig. 1 we show representatively the large-$`n`$ extrapolation of our data for $`k+1=3`$. Both, finite-$`n`$ energies and entropies, appear to scale linearly when plotted for $`1/n^{2/3}`$. The extrapolation results for the energies of BL and ORG are given in Tab. 1. The extrapolation results for the BL entropies are given in Tab. 2. Variation in the (estimated) errors reflect differences in computational effort (number of instances, largest $`n`$) and in the quality of the extrapolation eo\_bethe .

We can compare our results with existing theoretical predictions at the RS and the 1RSB level at least for the case of BL at $`k+1=3`$. For this case, Ref. MP2  reproduced $`e_3=1.2777`$ at the RS level, and yielded $`e_3=1.2717`$ at the 1RSB level (further replica corrections are expected to be small). These values and our extrapolation result of $`e_3=1.2716(1)`$ are indicated in Fig. 1. Clearly, the extrapolation result is extremely close to the 1RSB results, but inconsistent with the RS result.

We have also used EO to sample the degeneracy $`\mathrm{\Omega }`$ of the lowest-energy states found for BL. In these simulations we focused on smaller system sizes of $`n256`$ for $`k+1=3,\mathrm{},9`$ and $`10,14,\mathrm{},26`$ only. As a test for the accuracy of our implementation, we have run the simulation for $`k+1=3`$ a second time on identical instances, using different initial conditions and $`n/5`$ more updates, obtaining identical results for each instance.

When plotted for $`1/(k+1)`$ in Fig. 2, the BL entropies for even connectivities decrease about linearly toward zero. The entropies for odd connectivities, clearly non-zero at $`k+1=3`$ (see Fig. 1), drop more rapidly and are essentially indistinguishable from zero already at $`k+1=9`$. While any rapid, smooth decay could easily escape our limited accuracy, the plot still raises the question whether there may be a finite connectivity beyond which odd entropies vanish. The qualitative difference between even and odd BL entropies can be understood in the presents or absence, resp., of “free spins,” a finite fraction of spins in the ground state which violate exactly half of their bonds and may flip at no cost eo\_bethe .

We have also plotted all BL and ORG energies as $`e_c/\sqrt{c}`$ vs. $`1/c`$ (where $`c=k+1`$ for BL) in Fig. 3. We expect that $`lim_c\mathrm{}e_c/\sqrt{c}=E_{SK}=0.76321`$ for RSB CR ; MPV . All energies for ORG appear to fall on a single smooth curve. We can fit those energies with a parabola, projecting $`E_{SK}^{ORG}0.761`$. In contrast, the BL energies split into even and odd values, each located apparently on a simple line. Each line separately extrapolates very close to the exact value: $`E_{SK}^{even}0.763`$ and $`E_{SK}^{odd}0.765`$. Amazingly, the trivial value of $`e_2=1`$ is very close to the linear fit for the even results. Clearly, a function that would interpolate continuously all the BL energies will have to be very complicated (oscillatory). But we may speculate that its envelope for the even $`e_{k+1}`$ is a simple line, passing $`e_2=1`$ and the SK result:

$`_{k+1}=E_{SK}\sqrt{k+1}{\displaystyle \frac{2E_{SK}+\sqrt{2}}{\sqrt{k+1}}}.`$ (3)
In Fig. 4 we plot the deviation $`ϵ=e_{k+1}/_{k+1}1`$ of the extrapolated energies from Eq. (3) for even $`k+1`$. While the extrapolated values do not fall exactly onto the proposed function, they are all within about 0.2% of it. In fact, all points are slightly too high, which may indicate a more complex functional correction to Eq. (3), or a systematic error, say, in the extrapolation due to higher-order corrections.

It has been pointed out Mezard\_pc  that Eq. (3) would imply that a first-order perturbation around the trivial $`k+1=2`$ solution would be exact and give the RSB result for the SK model. But the obvious continuation of BL off the even integers fails to interpolate the data smoothly. If we “interpolate” BL for each $`k+1=2,4,6,\mathrm{}`$ with a mix of $`(1p)n`$ $`(k+1)`$-vertices and $`pn`$ $`(k+3)`$-vertices for $`0p1`$, the resulting energies only provide a set of secants, $`e_{k+1+2p}=pe_{k+3}+(1p)e_{k+1}`$, to the even-integer data. We have also plotted our (somewhat less accurate) extrapolation results for those interpolating graphs in Fig. 4. On this scale, the singular behavior of this continuation at the even integers becomes obvious.

We can further compare our extrapolated energies with perturbative calculations in the SK-limit of infinite connectivity PT ; DG ; LG . A recent RSB calculation indicates a $`1/c`$-correction for the ground state energy of $`f_1^{BL}0.32`$ for BL and $`f_1^{ORG}0.17`$ for ORG (see Figs. 3 and 4 of Ref. PT ). While our (crude) fit in Fig. 3 predicts a slope at the origin of $`0.16`$ for ORG, the slopes for either even or odd BL data would predict $`0.110.10`$, or $`0.1122`$ from Eq. (3), far from the perturbative result. It appears that the oscillation between even and odd $`k+1`$ complicates also the analytic continuation of the BL problem for large connectivities.

It will be most interesting to see how well upcoming RSB calculations at $`T=0`$ for even $`k+1`$ will correspond to the proposed function in Eq. (3), or to the extrapolated energies in Tabs. 1 in general. While the EO algorithm in itself can not provide information about the physics at $`T>0`$, the results presented in this Letter are sufficiently promising to apply EO also to sample other models eo\_prl  and more complicated properties of the ground states, such as overlap distributions and excitations.


text:
# Towards strangeness saturation in central heavy-ion collisions at high energies

## 1 Introduction

It is now well established that particle abundances can be described by statistical-thermal models. In such a way, a large number of observables can be reproduced with a small number of parameters, namely the temperature, baryo-chemical potential and a factor measuring the degree of strangeness saturation. Here we focus on the centrality dependence of the hadron multiplicities and adjust the thermal parameters so as to reproduce the experimental data.

The recent study in  (cf. table III therein) impressively demonstrated that, with increasing system size at SPS energies, the strangeness saturation increases. In  we have shown that at a beam energy of 158 AGeV, in collisions of lead-on-lead nuclei, the strangeness saturation continuously increases with centrality. However, strangeness is well below saturation. A preliminary analysis  of the centrality dependence at RHIC energy of $`\sqrt{s}_{NN}=130`$ GeV points to a further increase of strangeness towards saturation for central collisions of gold nuclei. An independent analysis  confirms this finding.

In this paper, we compare the analyses of two data sets: (i) NA49 $`4\pi `$ multiplicities of $`\pi `$, $`K^\pm `$, $`\varphi `$, $`N_{\mathrm{part}}`$ (taken as the sum over all baryons) and $`\overline{p}`$ in 6 centrality bins in the reaction Pb (158 AGeV) + Pb  (our previous study  did not include the $`\varphi `$ multiplicities). It should be emphasized that protons are not included in our analysis since, particularly in non-central collisions, there may be a large spectator component in the quoted experimental yields. (ii) PHENIX mid-rapidity densities of $`\pi ^\pm `$, $`K^\pm `$ and $`p^\pm `$ in the reaction Au + Au at $`\sqrt{s}=130`$ AGeV in 5 centrality bins . These yields were not corrected for weak decays. PHENIX estimates the probability for reconstructing protons from $`\mathrm{\Lambda }`$ decays as prompt protons at 32% at $`p_T`$ = 1 GeV/c . In the analysis of the PHENIX data, two fits were performed: with no feeding from weak decays, and with a 32% feed-down from $`\mathrm{\Lambda }`$ decay.

The use of $`4\pi `$ data at SPS energy is in the spirit of the fireball model, since many dynamical effects cancel out in ratios of fully-integrated hadron yields. In particular, effects due to flow disappear if the freeze-out surface is characterized by a single temperature and chemical potential .

## 2 Analyses of Hadron Multiplicities

Hadron multiplicities can be described  by the grand-canonical partition function $`𝒵(V,T,\stackrel{}{\mu }_i)=\text{Tr}\{\text{e}^{\frac{\widehat{H}\stackrel{}{\mu }_i\stackrel{}{Q}_i}{T}}\}`$, where $`\widehat{H}`$ is the statistical operator of the system, $`T`$ denotes the temperature, and $`\mu _i`$ and $`Q_i`$ represent the chemical potentials and corresponding conserved charges respectively. In the analysis of $`4\pi `$ data, the net-zero strangeness and the baryon-to-electric charge ratio of the colliding nuclei constrain the components of $`\stackrel{}{\mu }_i=(\mu _B,\mu _S,\mu _Q)`$. These constraints have to be relaxed when considering data in a limited rapidity window, increasing the number of free parameters. The particle numbers are given by

$$N_i^{\mathrm{prim}}=V(2J_i+1)\frac{d^3p}{(2\pi )^3}𝑑m_i\frac{1}{\gamma _s^{\left|S_i\right|}\text{e}^{\frac{E_i\stackrel{}{\mu }_i\stackrel{}{Q}_i}{T}}\pm 1}\text{BW}(m_i),$$
(1)
where we include phenomenologically a strangeness saturation factor $`\gamma _s`$ with $`\left|S_i\right|`$ the number of valence strange quarks and anti-quarks in species $`i`$  (e.g. $`\gamma _s`$ suppression for the kaons and $`\gamma _s^2`$ for $`\varphi `$) to account for incomplete equilibration in this sector, $`E_i=\sqrt{\stackrel{}{p}^{\mathrm{\hspace{0.17em}2}}+m_i^2}`$, and BW is the Breit-Wigner distribution. The particle numbers to be compared with experiment are $`N_i=N_i^{\mathrm{prim}}+_j\text{Br}^{ji}N_j^{\mathrm{prim}}`$, due to decays of unstable particles with branching ratios $`\text{Br}^{ji}`$.

The results of our fits are displayed in figs. 1 and 2. The inclusion of weak decays in the PHENIX analysis lowers the temperature but does not appreciably affect either $`\gamma _s`$ or $`\mu _B`$. Both at SPS and RHIC, the number of kaons and antikaons increase with centrality relative to the charged pion and antiproton multiplicities. This is responsible for the increase of the strangeness saturation factor $`\gamma _s`$. At RHIC energy the temperature also appears to increase with centrality, as shown in fig. 1.

The findings for central collisions at RHIC agree fairly well with the results in , as do our extracted values of $`T`$ and $`\mu _B`$ in central collisions at SPS energy. In contrast to , however, we find a pronounced strangeness under-saturation, as in  at SPS energy (cf. also  for other analyses with $`\gamma _s<1`$ at SPS). What is striking is the increase of $`\gamma _s`$ with energy and centrality, and the tendency to approach full strangeness saturation at RHIC, i.e. $`\gamma _s1`$. The observed increase of $`\gamma _s`$ with increasing centrality is in agreement with parton kinetics as described in , for example.

The baryo-chemical potential is fairly constant at both SPS and RHIC, as shown in fig. 2. The apparent centrality independence of the temperature and baryo-chemical potential at SPS leads to the model prediction that the ratio $`\varphi /K^+`$ should increase like $`\gamma _s`$ with centrality. This, however, does not agree with the experimentally observed centrality independence of this ratio. The trend can be reproduced by the model if the $`\varphi `$ is suppressed by a factor of $`\gamma _s`$ (instead of $`\gamma _s^2`$ ), but this leads to higher values of $`\chi ^2`$ in all but the most peripheral bin. This is shown quantitatively in fig. 2 where the $`\varphi /K^+`$ ratio is plotted against the number of participants. As can be seen, if the $`\varphi `$ mesons are suppressed only by a factor $`\gamma _s`$, one overestimates their number but the dependence as a function of $`N_{\mathrm{part}}`$ can be reproduced. If however, the $`\varphi `$ has a quadratic suppression factor, $`\gamma _s^2`$, then the normalization is reproduced, at least in the three most central bins, but the shape of the curve does not match the observed dependence. It would be of interest to establish this point more clearly, since it touches on the precise way in which strangeness suppression occurs. A more accurate determination of the $`\varphi `$ yield for small values of $`N_{\mathrm{part}}`$, especially for 150 - 200 participants, would be extremely useful. Also at AGS energies, the $`\varphi /K^+`$ ratio is centrality independent . The situation of the $`\varphi `$ yield as a function of $`N_{\mathrm{part}}`$ remains inconclusive.

## 3 Summary

In summary, we have shown that the analyses of full phase-space multiplicities at SPS energy and mid-rapidity densities at RHIC energy point to strangeness saturation in central collisions at still higher energies.