text:
# Coherent acceleration and Landau-Zener tunneling of Bose-Einstein condensates in 1-D optical lattices

## Abstract

We have loaded Bose-Einstein condensates into one-dimensional, off-resonant optical lattices and accelerated them by chirping the frequency difference between the two lattice beams. For small values of the lattice well-depth, Bloch oscillations were observed. Landau-Zener tunneling out of the lowest lattice band, leading to a breakdown of the oscillations, was also studied. In order to allow in-trap measurements of these phenomena, we dynamically compensated for the intrinsic micromotion of the atoms in our time-orbiting potential trap.

The properties of ultra-cold atoms in periodic light-shift potentials in one, two and three dimensions have been investigated extensively in the past ten years . In near-resonant and, more recently, far-detuned optical lattices, a variety of phenomena have been studied, such as the magnetic properties of atoms in optical lattices, revivals of wave-packet oscillations, and Bloch oscillations in accelerated lattices . While in most of the original optical lattice experiments the atomic clouds had temperatures in the the micro-Kelvin range, corresponding to a few recoil energies of the atoms, atomic samples with sub-recoil energies are now routinely produced in Bose-Einstein condensation experiments. Since the first experimental realizations in 1995, many aspects of Bose-Einstein condensed atomic clouds (BECs) have been studied , ranging from collective excitations to superfluid properties and quantized vortices. So far, the majority of these experiments have been carried out essentially in harmonic-oscillator potentials provided by magnetic traps or optical dipole traps. The properties of BECs in periodic potentials constitute a vast new field of research, for

instance . Several experiments in the pulsed standing wave regime  as well as studies of the tunneling of BECs out of the potential wells of a shallow optical lattice in the presence of gravity , and recent search for the superfluid dynamics have taken the first steps in that direction. In this paper, we present the results of experiments on BECs of $`{}_{}{}^{87}\mathrm{Rb}`$ atoms in accelerated optical lattices. In particular, we demonstrate coherent acceleration of BECs adiabatically loaded into optical lattices. For small values of the lattice depth we observed Bloch oscillations, which exhibited Landau-Zener breakdown when the lattice depth was further reduced and/or the acceleration increased. We loaded the condensate into optical lattices with different spatial periods, generating the periodic optical lattice either from two counterpropagating laser beams or two laser beam enclosing an angle $`\theta `$ different from $`180`$ degrees.

The properties of a Bose-Einstein condensate located in a periodic optical lattice with depth $`U_0`$ are described through the Gross-Pitaevskii equation valid for the single-particle wavefunction. Because the nonlinear term in the Gross-Pitaevskii equation reproduces the spatial periodicity of the wavefunction, the condensate ground state is periodic. In agreement with the Bloch approach, the condensate excitation spectrum exhibits a band structure and in presence of an acceleration of the optical lattice, Bloch oscillations of the condensate should occur . We present experimental results for Bloch oscillations preserving the condensate wavefunction. The role of the nonlinear interaction term of the Gross-Pitaevskii equation may be described through an effective potential in a noninteracting gas model . In the perturbative regime of ref.  the effective potential is $`U_{\mathrm{eff}}=U_0/(1+4C)`$ with $`C=g/E_\mathrm{B}`$ the interaction ratio between the nonlinear interaction term $`g=4\pi n_0\mathrm{}^2a/M`$ and the lattice Bloch energy $`E_\mathrm{B}=\mathrm{}^2(2\pi )^2/Md^2`$. The parameter $`C`$ contains the condensate density $`n_0`$, the s-wave scattering length $`a`$, the atomic mass $`M`$, the lattice constant $`d=\pi /\mathrm{sin}(\theta /2)k`$, with $`k`$ the laser wavenumber, and $`\theta `$ the angle between the two laser beams creating the 1-D optical lattice. From this it follows that a small angle $`\theta `$ should result in a large interaction term $`C`$. In the following, the parameters $`d`$, $`E_\mathrm{B}`$ and $`C`$ always refer to the respective lattice geometries with angle $`\theta `$.

Our apparatus used to achieve Bose-Einstein condensation of $`{}_{}{}^{87}\mathrm{Rb}`$ is described in detail in . Essentially, $`5\times 10^7`$ atoms captured in a magneto-optical trap (MOT) were transferred into a triaxial time-orbiting potential trap (TOP) . Subsequently, the atoms were evaporatively cooled down to the transition temperature for Bose-Einstein condensation, and after further cooling we obtained condensates of $`10^4`$ atoms without a discernible thermal component in a magnetic trap with frequencies around $`1530\mathrm{Hz}`$. In one set of experiments, the magnetic trap was then switched off and a horizontal 1-D optical lattice was switched on, while in the other case the interaction between the condensate and the lattice took place inside the magnetic trap, which was subsequently switched off to allow time-of-flight imaging. The lattice beams were created by a $`50\mathrm{mW}`$ diode slave-laser injected by a grating-stabilized master-laser blue-detuned by $`\mathrm{\Delta }2835\mathrm{GHz}`$ from the $`{}_{}{}^{87}\mathrm{Rb}`$ resonance line. After passage through an optical fibre, the laser light was split and passed through two acousto-optic modulators (AOMs) that were separately controlled by two phase-locked RF function generators operating at frequencies around $`80\mathrm{MHz}`$, with a frequency difference $`\delta `$. The first-order output beams of the AOMs generated the optical lattice. An acceleration of the lattice was effected by applying a linear ramp to $`\delta `$. For the values of the detuning and laser intensity used in our experiment, the spontaneous photon scattering rate ($`10\mathrm{s}^1`$) was negligible during the interaction times of a few milliseconds. In our experimental setup, we realized a counter-propagating lattice geometry with $`\theta =180\mathrm{deg}`$ and an angle geometry with $`\theta =29\mathrm{deg}`$, leading to lattice constants $`d`$ of $`0.39`$ and 1.56 $`\mu \mathrm{m}`$, respectively. The typical condensate density of $`n_05\times 10^{13}\mathrm{cm}^3`$ for our trap parameters in the Thomas-Fermi limit leads to $`C=0.01`$ for the counterpropagating configuration and $`C=0.17`$ for $`\theta =29\mathrm{deg}`$ . We therefore expect mean-field effects to be negligible in the counter-propagating lattice geometry, whereas in the geometry with a larger lattice constant the effective potential acting on the condensate should be significantly reduced. The two lattice beams with $`3\mathrm{mW}`$ each were expanded to a waist of $`1.8\mathrm{mm}`$, giving a theoretical maximum lattice depth (for $`\mathrm{\Delta }28\mathrm{GHz}`$) of $`U_04E_\mathrm{B}`$ for the counter-propagating lattice geometry.

In a preliminary experiment aimed at determining the depth of the periodic optical potential, we flashed on the counterpropagating lattice with $`\delta =\mathrm{}k^2/(2\pi M)=15.08\mathrm{kHz}`$ for $`10400\mu \mathrm{s}`$. This detuning corresponds to the first Bragg resonance, causing the condensate to undergo Rabi oscillations between the momentum states $`|p=0`$ and $`|p=2\mathrm{}k`$ (see Fig. 1). From the measured Rabi frequency we could then determine the lattice depth . The results of those measurements fell short by a factor of about 2 of the theoretical prediction, which was mainly due to the $`1015\%`$ uncertainty in our intensity measurements and imperfections in the alignment and polarization of the lattice beams. As shown above, mean-field effects are expected to be negligible in this geometry.

In order to accelerate the condensate, we adiabatically loaded it into the lattice by switching one of the lattice beams on suddenly and ramping the intensity of the other beam from $`0`$ to its final value in $`200\mu \mathrm{s}`$ . Thereafter, the linear increase of the detuning $`\delta `$ provided a constant acceleration $`a=\frac{\lambda }{2\mathrm{sin}(\theta /2)}\frac{d\delta }{dt}`$ of the optical lattice, leading to a final lattice velocity $`v_{lat}=\frac{\lambda }{2\mathrm{sin}(\theta /2)}\delta `$, where $`\delta `$ is the final detuning between the beams. After a few milliseconds of acceleration, the lattice beams were switched off and the condensate was imaged after another $`1015\mathrm{ms}`$ of free fall. As the lattice can only transfer momentum to the condensate in units of the Bloch momentum $`p_\mathrm{B}=\mathrm{}(2\pi /d)`$, the acceleration of the condensate showed up as diffraction peaks corresponding to higher momentum classes as time increased (Fig. 2). Since for our magnetic trap parameters the initial momentum spread of the condensate (which is transferred into a spread of the lattice quasimomentum during an adiabatic switch-on) was much less than a recoil momentum of the optical lattice, the different momentum classes $`|p=\pm np_\mathrm{B}`$ (where $`n=0,1,2,\mathrm{}`$) occupied by the condensate wavefunction could be resolved directly after the time-of-flight. The average velocity of the condensate was derived from the occupation numbers of the different momentum states. Figure 3 shows the results of the acceleration of a condensate in the counterpropagating lattice with $`U_{\mathrm{eff}}=0.29E_\mathrm{B}`$ and $`a=9.81\mathrm{m}\mathrm{s}^2`$. In the rest-frame of the lattice (Fig. 3 (b)), one clearly sees Bloch oscillations of the condensate velocity corresponding to a Bloch-period $`\tau _B=\frac{h}{M_{Rb}ad}=1.2\mathrm{ms}`$. The shape of these oscillations agrees well with the theoretical curve calculated from the lowest energy band of the lattice. As described in , the acceleration process within a periodic potential can also be viewed as a succession of adiabatic rapid passages between momentum states $`|p=2n\mathrm{}k`$ and $`|p=2(n+1)\mathrm{}k`$. We observed a momentum transfer of up to $`6p_B`$ without a detectable reduction of the phase-space density of the condensate .

When we increased the acceleration of the lattice or decreased the lattice depth, not all of the condensate was coherently accelerated up to the final velocity of the lattice. This can be interpreted in terms of Landau-Zener tunneling of the condensate out of the lowest band when the edge of the Brillouin zone is reached. Each time the condensate is accelerated across this edge, a fraction $`r=\mathrm{exp}(a_c/a)`$ with $`a_c=\frac{\pi U_0^2}{8\mathrm{}p_B}`$ undergoes tunneling into the first excited band (and, therefore, effectively to the continuum, as the gaps between higher bands are negligible for the shallow potentials used here). In Fig. 4 (a), the average velocity of the condensate after acceleration of the counter-propagating lattice to the Bloch velocity $`v_\mathrm{B}=p_\mathrm{B}/M`$ is shown as a function of acceleration and lattice depth along with a theoretical prediction using the Landau-Zener tunneling probability, which gives a mean velocity $`v_m=(1r)v_\mathrm{B}`$ of the condensate for a final velocity $`v_\mathrm{B}`$ of the lattice. When the lattice is accelerated up to a velocity $`nv_\mathrm{B}`$, a straightforward generalization of this equation yields $`v_m=v_\mathrm{B}(1/r1)[1(1r)^n]`$ for the mean velocity of the condensate, which is in good agreement with our experimental findings (Fig. 4 (b)) . In this case, a fraction $`r`$ of the condensate undergoes Landau-Zener tunneling each time the Bragg resonance is crossed, with a remaining fraction $`1r`$ of the condensate being accelerated further. Agreement with theory is also good when, instead of changing the acceleration, the lattice depth is varied at fixed acceleration (Fig. 4(c)).

Similar experiments were performed in the geometry leading to the larger lattice constant of $`1.56\mu \mathrm{m}`$. Because of the reduced Bloch velocity $`v_\mathrm{B}`$ in this geometry, the acceleration process was extremely sensitive to any initial velocity of the condensate, which in our TOP trap is intrinsically given by the micromotion  at the frequency of the bias field. For the trap parameters used in our experiments, the velocity amplitude of the micromotion could be of the same order of magnitude as $`v_\mathrm{B}`$ and the condensates could, therefore, have quasimomenta close to the edge of the Brillouin zone. In fact, when the standing wave was flashed on as in the Rabi-oscillation measurement described above, Bragg diffraction could be observed for zero detuning between the lattice beams, with the diffraction efficiency depending on the initial velocity of the condensate. Moreover, an initial velocity close to the band edge would have made it impossible to switch on the lattice adiabatically. In order to counteract these effects, we performed the acceleration experiments inside the magnetic trap, eliminating the velocity of the condensate relative to the lattice by phase-modulating one of the lattice beams at the same frequency and in phase with the rotating bias field of the TOP trap. In this way, in the rest frame of the lattice the micromotion was compensated. Nevertheless, a residual sloshing of the condensate with amplitudes $`<3\mu \mathrm{m}`$ could not be ruled out, so that the uncertainty in the initial velocity of the condensate was still around $`0.5\mathrm{mm}\mathrm{s}^1`$, corresponding to $`0.3v_\mathrm{B}`$ in this geometry. In Fig. 3 (c), the results of the acceleration of the condensate with a nominal lattice depth of $`1.38E_\mathrm{B}`$ are shown together with the theoretical curves for $`U_0=1.38E_\mathrm{B}`$ and the (assumed) effective potential $`U_{\mathrm{eff}}0.65U_00.88E_\mathrm{B}`$. The Landau-Zener tunneling probabilities measured in this lattice geometry were compatible with the same effective potential assuming that the conversion factor between theoretical and actual lattice depth was the same as in the counter-propagating geometry in which $`U_{\mathrm{eff}}U_0`$. In order to demonstrate unequivocally the reduction of the effective lattice potential by the interaction term, however, it would be necessary to vary $`n_0`$ appreciably holding all other parameters constant, which in our experimental setup was not possible without creating systematic errors due to variations in the equilibrium position of the trap (and hence the local laser intensity in the Gaussian profile) when the trap frequency was changed. As the interaction term is expected to distort the band structure of the condensate in the lattice , it should affect all measurable quantities (Rabi frequency, amplitude of Bloch oscillations, and tunneling probablility ) in the same way, so that a differential measurement is necessary (as has been demonstrated in the pulsed Bragg-diffraction regime of Ref. ). On the theoretical side, the finite extent of the condensate leading to the occupation of only a few lattice sites and the three-dimensional nature of the condensate evolution as well as the role of the interaction term in the adiabaticity criterion for switching on the lattice will also have to be taken into account.

In summary, we have investigated the coherent acceleration of Bose-Einstein condensates adiabatically loaded into a 1-D optical lattice as well as Bloch oscillations and Landau-Zener tunneling out of the lowest Bloch band. The results obtained are in good agreement with the available theories and extend the corresponding work on ultra-cold atoms in optical lattices into the domain of Bose-Einstein condensates.

This work was supported by the MURST through the PRIN2000 Initiative, and by the EU through the Cold Quantum-Gases Network, contract HPRN-CT-2000-00125. O.M. gratefully acknowledges a Marie-Curie Fellowship from the EU within the IHP Programme.


format: text

text:
# Does Kolmogorov mean field theory become exact for turbulence above some critical dimension?

## Abstract

I give three different arguments for an upper critical dimension $`d_{max}>3`$ above which the 1941 Kolmogorov mean field theory becomes essentially exact, and anomalous scaling vanishes. The first argument concerns the number of degrees of freedom in a turbulent flow and indicates that $`d_{max}=4`$. The second argument is a naive estimate of dangerous fluctuations, and also suggests that $`d_{max}=4`$. The third argument is related to a known critical point of the GOY shell model when the amplitude of back energy transfer becomes small. This third argument does not give a numerical value for $`d_{max}`$. None of these arguments bears any known relationship to any of the others nor to the generally accepted qualitative physical picture of the dynamical origin of anomalous scaling in turbulence. Despite this, the three arguments together suggest that the suggestion of an upper critical dimension should be taken seriously.

Physics Department, New York University, New York, NY 10003,
and Levich Institute, CCNY, New York, NY 10031 USA
Electronic address: Mark.Nelkin@nyu.edu

The 1941 Kolmogorov theory (K41)  is a good first approximation to the universal statistics of the small scale velocity fluctuations in high Reynolds number turbulence. In this theory, the probability distribution function (pdf) of the velocity difference $`\mathrm{\Delta }u(r)`$ between two points separated by a distance $`r`$ is self-similar and nearly Gaussian. It is well established from experiment and direct numerical simulation (DNS) that this self-similarity is only approximate in three dimensional turbulence. As the scale separation $`r`$ is decreased, the tails of the pdf become enhanced. The deviations from self-similarity are known as anomalous scaling.

Most studies of anomalous scaling in the inertial range have concentrated on the longitudinal structure functions,

$$D_p(r)=[u(x+r)u(x)]^p=[\mathrm{\Delta }u(r)]^p,$$
(1)
which are the moments of the velocity difference $`\mathrm{\Delta }u(r)`$ between two points separated by a distance $`r`$. The term “longitudinal” refers to the component of the velocity difference along the line joining the two points. These structure functions are characterized by power-law behavior in the inertial range with anomalous scaling exponents $`\zeta _p`$ defined by

$$D_p(r)=Const.ϵ^{p/3}r^{p/3}(r/L)^{(\zeta _pp/3)}.$$
(2)
The “normal” scaling of K41 corresponds to $`\zeta _p=p/3`$. The observed values of $`\zeta _p`$ for $`p>4`$ are smaller than the K41 values corresponding to an intermittent form for the pdf where large values of $`\mathrm{\Delta }u(r)`$ are increasingly enhanced as $`r`$ is decreased. This phenomenon is often called the internal intermittency of turbulent flows.

It was natural to suggest in the 1970’s  that K41 theory played the role of a mean field theory for turbulence. By analogy with equilibrium critical phenomena, it might be that this mean field theory could become exact either below some critical spatial dimension $`d_{min}`$ or above some critical dimension $`d_{max}`$. There have been many speculations about this possibility, but until recently these speculations have not led to any serious physical insight. I suggest here that recent developments might change this situation.

A very important development has been the observation that the inverse cascade of energy in two dimensional turbulence shows normal K41 scaling. This has been observed in the laboratory by Paret and Tabeling . It was first seen in DNS by Smith and Yakhot , explored in more detail by Boffetta et al. , and has been given some theoretical understanding by Yakhot . Recently Yakhot  has revived the idea of continuously varying the dimension $`d`$ between $`d=3`$ and $`d=2`$, and has suggested that there is a lower critical dimension $`d_{min}`$ at which the direction of the energy flux changes sign. This idea has been given important and interesting support within the GOY shell model in a recent paper by Giuliani, Jensen and Yakhot . They find that intermittency does not vanish near this critical dimension even though it is absent in $`d=2`$. I return to a more quantitative discussion of this model later in the present paper.

It is generally agreed, as first suggested by Kraichnan , that anomalous scaling arises from a competition between a spatially local cascade, and a spatially non-local dynamical averaging of coupled local cascades. This spatial mixing is intimately associated with the non-local pressure term in the Navier-Stokes equation. This competition is expected to be a smooth function of spatial dimension with intermittency decreasing as the spatial dimension increases. There are indications, however, that a crossover dimension $`d_{max}`$ may exist with K41 mean field theory becoming exact for $`d>d_{max}`$. One such indication comes from a calculation of the number of degrees of freedom in a turbulent flow. Using the $`\beta `$ model of fractally homogeneous turbulence , Kraichnan  showed that anomalous scaling decreased the number of degrees of freedom for $`d<4`$, but increased the number of degrees of freedom for $`d>4`$. Meneveau and Nelkin  extended this result to a general multifractal model. The crossover occurs because the effects of anomalous scaling on the density of modes in the inertial range competes with the effects of anomalous scaling on the fluctuating viscous cutoff, and this competition depends on the spatial dimension $`d`$.

It is appealing to think that there is some physics in the above estimate. Anomalous scaling is often thought to be associated with coherent structures in the flow, and therefore with a reduction of randomness, and a reduction in the number of degrees of freedom. At first sight it does not seem physically plausible for anomalous scaling to increase the number of degrees of freedom. Furthermore it is important to remember that anomalous scaling is a property of a statistically stationary state, and need not have any direct dynamical origin. On the other hand, there is no known reason why the number of degrees of freedom in a turbulent flow should be of direct physical interest. Furthermore the calculation of the crossover dimension depends on the mechanism of viscous dissipation. If the viscous term proportional to $`\nu k^2`$ in the Navier-Stokes equation is replaced by a hypo-viscous term proportional to $`\lambda k^p`$ with $`p<2`$, Nelkin and Meneveau  showed that the same crossover effect calculated for ordinary viscosity at $`d=4`$ occurs for $`d=3`$ when $`p=5/3`$. This suggests looking at a DNS of three dimensional turbulence with varying hypoviscous index $`p`$ to see if any reduction of intermittency can be seen when $`p<2`$.

There is another argument which suggests a crossover at $`d=4`$. The anomalous scaling in turbulence is associated with anomalous fluctuations of the energy dissipation rate which are, in turn, related to the fluctuations of the nonlinear energy transfer through Kolmogorov’s refined similarity hypothesis . To identify “dangerous” fluctuations, we need a locally defined dynamical variable which is a reasonable surrogate for the nonlinear energy transfer. A possible choice is

$$T(𝐫)=\eta ^2(u/x)^3$$
(3)
where $`\eta `$ is the dissipation length scale, and $`u`$ is the local velocity. This variable is “regularized” in the sense that its average value is constant in the limit of zero viscosity. Since the scaling dimension of the velocity derivative is $`(2/3)`$, the scaling dimension of $`T(𝐫)`$ is $`6(2/3)=4`$. This suggests an upper critical dimension $`d_c=4`$. The preceding argument was given in an earlier paper . It is not clear if it is valid, but it is not obviously wrong. There is a quite different sense in which this argument might be “dangerous.” It admittedly arises from a search for a dynamical variable which might show the desired crossover properties.

Finally there are interesting crossover properties observed in the GOY   shell model of turbulence. Since this model is a one dimensional dynamical model, it does not make any explicit reference to the spatial dimension $`d`$, but there are reasonable arguments  which suggest how the parameters in the model might vary with $`d`$. The starting point is a set of wave numbers $`k_n=k_02^n`$ and an associated set of complex amplitudes $`u_n`$ of the velocity field. Each amplitude interacts with nearest and next-nearest neighbor shells, and the corresponding set of coupled ODE’s takes the form:

$$du_n/dt=\nu k_n^2u_n+f_n+C_n$$
(4)
where the nonlinear coupling term is given by

$$C_n=ik_n[a_nu_{n+1}^{}u_{n+2}^{}+(b_n/2)u_{n1}^{}u_{n+1}^{}+(c_n/4)u_{n1}^{}u_{n2}^{}].$$
(5)
I take my notation from .

The values of the coupling constants are fixed by imposing conserved quantities. To conserve the total energy when $`f_n=\nu =0`$, we require $`a_n+b_{n+1}+c_{n+2}=0`$. The time scale is fixed by the condition that $`a_n=1`$ leaving the free parameter $`\delta `$ by defining the coupling constants as

$$a_n=1,b_n=\delta ,c_n=(1\delta ).$$
(6)
It has been known for some time that the choice $`\delta =1/2`$ gives good agreement between the temporal intermittency in the GOY model and the observed spatial intermittency of thee dimensional turbulence. It was pointed out by Kadanoff et al  that this choice of $`\delta `$ gives a second conserved quadratic quantity which corresponds loosely to helicity. It is not at all understood, however, why the temporal intermittency in this one dimensional model should be related to the spatial intermittency in real three dimensional turbulent flows.

The GOY model can be generalized to allow the scale ratio $`\lambda `$ between adjacent shells to be different than $`2`$. It was pointed out by Kandanoff et al  that the agreement with experiment remains good for varying values of $`\lambda `$ as long as the second conserved quadratic quantity continues to correspond to helicity.

This observation strongly suggests that there could be physical interest in varying the parameter $`\delta `$. The recent paper by Giuliani, Jensen and Yakhot  shows that as $`\delta `$ is increased, there is a critical value at $`\delta =1`$ at which the direction of the energy cascade changes. The value $`\delta =5/4`$ corresponds to enstrophy as the second conserved quadratic quantity, and might be related to the inverse cascade in two dimensional turbulence. This supports the idea that increasing $`\delta `$ from $`\delta =1/2`$ corresponds in some sense to reducing the dimension in turbulence from $`d=3`$ towards $`d=2`$.

In this paper I would like to emphasize what happens when $`\delta `$ is reduced below $`\delta =1/2`$. In the same spirit as above, this might correspond to increasing the spatial dimension to $`d>3`$. It was shown by Biferale et al  that there is a critical value of $`\delta 0.38`$. For values of $`\delta `$ smaller than this, the K41 fixed point becomes stable and intermittency vanishes. Further studies have shown that the properties of this critical point are quite subtle and can be affected by viscosity , but I choose here to ignore these subtleties. Clearly this suggests that there is an upper critical dimension $`d_{max}>3`$ above which intermittency vanishes, and the K41 mean field theory becomes essentially correct.

I should emphasize that the dynamical mechanism in the GOY model for the vanishing of intermittency is quite different than the conventionally assumed competition between local cascade and spatial mixing . The one dimensional GOY model does not consider space at all. As $`\delta `$ is decreased, the relative amplitude of backwards energy transfer is decreased, and this reaches a critical value below which intermittency is no longer maintained.

An interesting physical argument concerning the origin of this critical point has been given by Biferale and Kerr . They suggest that the fluctuations in the second conserved quadratic quantity become “dangerous” when $`\delta `$ is larger than this critical value. The fluctuations are “dangerous” when they are dominated by small scale fluctuations in the K41 solution. This estimate of the transition value of $`\delta `$ agrees well with numerical experiment for varying values of the scale ratio $`\lambda `$.

In conclusion, I have suggested that there is a critical dimension $`d_{max}>3`$ above which anomalous scaling in turbulence vanishes and the K41 mean field theory becomes essentially exact. Counting degrees of freedom, and naive arguments about “dangerous” fluctuations suggest that $`d_{max}=4`$. The existence of a transition to a stable K41 fixed point in the one dimensional GOY shell model has been known for some time, and it suggests that $`d_{max}>3`$, but gives no indication of its numerical value. In all of these arguments, there is no connection made to the usual picture of the dynamical origin of intermittency nor is there any known connection among the three arguments presented here. Despite these limitations, it is suggestive that there is important physics in turbulence for $`d>3`$ which deserves serious further study.


text:
# Interfaces (and Regional Congruence?) in Spin Glasses

## Abstract

We present a general theorem restricting properties of interfaces between thermodynamic states and apply it to the spin glass excitations observed numerically by Krzakala-Martin and Palassini-Young in spatial dimensions $`d=3,4`$. We show that such excitations, with interface dimension $`d_s<d`$, cannot yield regionally congruent thermodynamic states. More generally, zero density interfaces of translation-covariant excitations cannot be pinned (by the disorder) in any $`d`$ but rather must deflect to infinity in the thermodynamic limit. Additional consequences concerning regional congruence in spin glasses and other systems are discussed.

Recent numerical studies  on finite-dimensional spin glasses observed unanticipated low-energy excitations that have generated considerable interest . The models examined were nearest-neighbor Ising spin glasses with the Edwards-Anderson  Hamiltonian

$$=\underset{<x,y>}{}J_{xy}\sigma _x\sigma _y,$$
(1)
where the couplings $`J_{xy}`$ are independently chosen from a Gaussian distribution with mean zero and variance one, the sum is over only nearest neighbors on the $`d`$-dimensional cubic lattice, and the spins $`\sigma _z=\pm 1`$.

These studies examined the ground states $`\pm \sigma ^L`$ and their low-energy excitations corresponding to the coupling realization $`𝒥^L`$ inside cubes $`\mathrm{\Lambda }_L`$ of side $`L`$, centered at the origin, with periodic boundary conditions. Krzakala and Martin (KM) , working in $`d=3`$, forced a random pair of spins, $`\sigma _z,\sigma _z^{}`$, to assume a relative orientation opposite to the one in $`\sigma ^L`$, and then relaxed the rest of the spins to a new lowest-energy configuration. This ensures that at least some bonds in the excited spin configuration, $`\sigma ^L`$, must be changed (i.e., satisfied $``$ unsatisfied) from $`\sigma ^L`$. It also ensures that the energy of $`\sigma ^L`$ is no more than $`O(1)`$ above that of $`\sigma _L`$, regardless of $`L`$. Interestingly, rather than simply generating local droplet flips, the KM procedure was observed to yield large-scale (i.e., of lengthscale $`l=O(L)`$) excitations $`\sigma ^L`$, with $`O(L^d)`$ spins flipped from $`\pm \sigma ^L`$, whose interface with $`\sigma ^L`$ (i.e., the set of bonds satisfied in $`\sigma ^L`$ but not in $`\sigma ^L`$, or vice-versa) had the property that the number of interface bonds scaled as $`L^{d_s}`$, with $`d_s<d`$.

Palassini and Young (PY) , working in both $`d=3`$ and $`4`$, excited $`\sigma ^L`$ differently by adding a novel coupling-dependent bulk perturbation to the Hamiltonian (1) in $`\mathrm{\Lambda }_L`$. PY also interpreted their results as evidence of an excited $`\sigma ^L`$ with $`d_s<d`$. At the same time the interface energy scaled as $`l^\theta ^{}`$ ($`L^\theta ^{}`$), with $`\theta ^{}=0`$ ; i.e., it remains finite as $`l\mathrm{}`$ (presuming that $`l`$ continues to grow as $`O(L)`$).

It is generally agreed  that the excitations of KM and of PY correspond to the same physical objects. We hereafter refer to them as KMPY excitations. Given a ground state $`\sigma ^L`$ (or ground state pair $`\pm \sigma ^L`$) in $`\mathrm{\Lambda }_L`$ with coupling-independent boundary conditions (such as periodic), KMPY excitations are global spin excitations whose structural and energetic properties scale with their length in a unique manner. More precisely, they are excitations $`\sigma ^L`$ that are characterized by three properties: 1) they have $`O(L^d)`$ spins flipped from $`\pm \sigma ^L`$ (and thus also have lengthscale $`l=O(L)`$); 2) the dimension $`d_s`$ of their interface with $`\sigma ^L`$ satisfies $`d_s<d`$; and 3) their energy difference with $`\sigma ^L`$ (i.e., the interface energy) scales as $`l^\theta ^{}`$ with $`\theta ^{}=0`$; i.e., the excitation energy remains of order one independently of $`l`$ (and $`L`$).

Questions have been raised over the correct interpretation of these numerical results , and correspondingly whether KMPY excitations really exist in the spin glass phase . We do not address these questions here. Rather, we take the point of view that if KMPY excitations do exist, then their physical meaning and implications for the low-temperature spin glass phase could be fundamentally important in the physics of disordered systems.

However, given the recentness of the discovery of KMPY excitations, their physical meaning and relevance remain unclear. By showing here rigorously that their interfaces cannot be pinned by the quenched disorder, we conclude that they cannot yield new ground or pure states. Such restrictions on their large-scale structure clarify the physical role they might play at low temperature $`T`$ and provide general results on a type of ground and/or pure state multiplicity not heretofore investigated: the possibility in spin glasses of regional congruence  (see below).

Pinning. A crucial question about a new type of excitation is whether its boundary is pinned. This has not yet been addressed for KMPY interfaces, to our knowledge.

To understand pinning in this context, consider at $`T=0`$ an increasing sequence of $`\mathrm{\Lambda }_L`$’s. For each $`L`$, use the procedures of  or  to create KMPY excitations; e.g. , first generate the periodic boundary condition ground state $`\sigma ^L`$, and then add a bulk perturbation (Eq. (2) of ) to the Hamiltonian to generate $`\sigma ^L`$, a perturbed ground state (i.e., an excited state for the original Hamiltonian). Then study the bonds $`x,y`$ that obey

$$\sigma _x^L\sigma _y^L=\sigma _x^L\sigma _y^L,$$
(2)
i.e., those that are satisfied in one state but not the other. The corresponding set of bonds in the dual lattice comprises the finite-size “domain wall” (for that $`L`$) or interface of the excitation (i.e., the boundary of the set of spins that are flipped to go from $`\sigma ^L`$ to $`\sigma ^L`$).

By pinning we mean the following. Consider a fixed $`L_0`$ (which can be arbitrarily large). Apply the KM or PY procedure on cubes $`\mathrm{\Lambda }_L`$, with $`LL_0`$. Observe $`\sigma ^{(L,L_0)}`$, the excited spin configuration $`\sigma ^L`$ restricted to $`\mathrm{\Lambda }_{L_0}`$. If the excitation interface remains inside $`\mathrm{\Lambda }_{L_0}`$ as $`L\mathrm{}`$ , then the interface is pinned. Pinning of lower-dimensional interfaces by quenched disorder is known to occur, e.g., in disordered ferromagnets for sufficiently large $`d`$ ; this example, and its relevance to spin glasses, will be discussed further below.

If the interface is not pinned, we say it “deflects to infinity”. Here, for any fixed $`L_0`$, the interface, for all $`L`$ above some $`L^{}`$, will not enter $`\mathrm{\Lambda }_{L_0}`$. (This is what occurs with interface ground states in disordered ferromagnets for small $`d`$ .) See Fig. 1 for a schematic illustration.

Scenarios. Assuming now that KMPY excitations do exist on at least small to moderate lengthscales, there are four possibilities for larger lengthscales. One is that the excitations disappear altogether, noted as a possibility in . Of the three remaining (more interesting) scenarios, two have the KMPY interface pinned, and the third has it deflect to infinity.

Scenario (1): KMPY excitations are pinned and give rise to new ground states at $`T=0`$ and new pure states at $`T>0`$ (below some $`T_c`$). I.e., their interfaces with the original ground state, or with each other, become relative domain walls between distinct ground (or pure) states not related by a global spin flip. (For more details on relative domain walls, see .) This interesting picture differs from previous proposals for the low-$`T`$ spin glass phase, and would be the first example of regional congruence  in spin glasses. Two distinct thermodynamic states that are not global flips of each other are regionally congruent if the total relative domain wall density vanishes; i.e., essentially the condition above that $`d_s<d`$. (If the density is nonzero so that $`d_s=d`$, as is more usually supposed, the states are said to be incongruent.)

Scenario (2): The excitations are pinned, but do not give rise to new ground states. For this to occur, any $`L\mathrm{}`$ limit of $`\sigma ^L`$ must be energetically unstable (in the Hamiltonian (1)) to the flip of some fixed finite droplet. Even if this occurs, it could be that KMPY excitations, at $`T>0`$, still give rise to new pure states.

Scenario (3): KMPY excitations persist on all lengthscales, but their interfaces deflect to infinity, as in Fig. 1. If this occurs, they cannot give rise to new thermodynamic ground or pure states, but could still be relevant to the excitation spectrum in finite volumes.

Which of these scenarios actually occurs was not addressed in , although  as well as  implicitly or explicitly assume Scenario (1). Determining which occurs is crucial to understanding the role of these excitations in the spin glass phase.

We now prove that KMPY excitations cannot be pinned by the quenched disorder, ruling out Scenarios (1) and (2). (A heuristic argument against Scenario (1) was presented in .) The remaining possibility (in which KMPY excitations persist on all lengthscales) is Scenario (3), where they deflect to infinity. Before presenting our theorem, we introduce the concept of metastate, which is implicit in the theorem and useful for understanding its applications.

Metastates. A ($`T=0`$) metastate gives the probability of finding various ground state pairs in typical large volumes. There are different constructions (technical details are in ); we give here only a simple physical description.

Consider the cube $`\mathrm{\Lambda }_{L_0}`$ with (large) fixed $`L_0`$, and then examine a sequence of $`\mathrm{\Lambda }_L`$’s with $`L\mathrm{}`$, all with (for example) periodic boundary conditions. Look, for each $`L`$, at the part of the ground state $`\sigma ^L`$ inside the smaller fixed $`\mathrm{\Lambda }_{L_0}`$, and keep a record of the fractions of $`L`$’s in which different spin configurations appear inside $`\mathrm{\Lambda }_{L_0}`$. (If there’s only a single pair of thermodynamic ground states, then one eventually sees one fixed configuration inside $`\mathrm{\Lambda }_{L_0}`$, and thus asymptotically the ground state pair appears inside $`\mathrm{\Lambda }_{L_0}`$ for a fraction one of the $`L`$’s). The resulting ($`T=0`$, periodic boundary condition) metastate is a probability measure on infinite-volume ground state pairs that provides information on the fraction of $`L`$’s for which each of them appears inside $`\mathrm{\Lambda }_{L_0}`$.

For the coming theorem, we extend the notion of metastate to uniform perturbation metastate . For each $`L`$, consider, in addition to $`\sigma ^L`$, a second $`\sigma ^L`$ (which in practice will be the ground state of a perturbed Hamiltonian) chosen so as to respect the finite-volume (torus) translation-invariance already present for $`(𝒥^L,\sigma ^L)`$ due to the periodic boundary conditions. (We show later that this holds for $`\sigma ^L`$ in both KM and PY constructions.) Then do for the pair $`(\sigma ^L,\sigma ^L)`$ what was done for $`\sigma ^L`$ in the original metastate. The resulting uniform perturbation metastate gives for both (infinite-volume) ground and excited states their relative frequency of appearance inside large volumes. The theorem presents implications of translation-invariance on the resulting interfaces, and shows that if there is a (pinned) interface at all, it must have strictly positive density (so $`d_s=d`$).

Theorem. On the cube $`\mathrm{\Lambda }_L`$ with periodic boundary conditions, for a given coupling configuration $`𝒥^L`$ let $`\sigma ^L`$ and $`\sigma ^L`$ be a pair of spin configurations such that the joint distribution of $`(𝒥^L,\sigma ^L,\sigma ^L)`$ is invariant under (torus) translations of $`\mathrm{\Lambda }_L`$. Let $`(𝒥,\sigma ,\sigma ^{})`$ be any limit in distribution as $`L\mathrm{}`$ of $`(𝒥^L,\sigma ^L,\sigma ^L)`$. Then for almost every $`(𝒥,\sigma ,\sigma ^{})`$ either $`\sigma ^{}=\pm \sigma `$ or else $`\sigma `$ and $`\sigma ^{}`$ have a relative interface of strictly positive density.

Proof. Because the joint distribution $`\kappa ^L`$ of $`(𝒥^L,\sigma ^L,\sigma ^L)`$ is for every $`L`$ invariant under torus translations, any limiting distribution $`\kappa `$ of $`(𝒥,\sigma ,\sigma ^{})`$ is invariant under all translations of the infinite-volume cubic lattice. The translation-invariance of $`\kappa `$ allows its decomposition into components in which translation-ergodicity holds (see, e.g., ). For each bond $`x,y`$ consider the event $`A_{x,y}`$ that $`x,y`$ is satisfied in one but not the other of $`\sigma ,\sigma ^{}`$. In each ergodic component, either the probability of $`A_{x,y}`$ equals zero and then there is no interface (i.e., $`\sigma ^{}=\pm \sigma `$) or else it equals some $`\rho >0`$. In the latter case, by the spatial ergodic theorem (see, e.g., ) the spatial density of $`x,y`$’s such that $`A_{x,y}`$ occurs must equal $`\rho `$, i.e., the $`(\sigma ,\sigma ^{})`$ interface has a nonzero density. Thus there is zero probability (with respect to $`\kappa `$) of a $`(𝒥,\sigma ,\sigma ^{})`$ such that there is a $`\sigma ,\sigma ^{}`$ interface, but with zero density.

Remarks. 1) If one takes the $`\kappa (𝒥,\sigma ,\sigma ^{})`$ of the theorem and conditions it on the coupling realization $`𝒥`$, the resulting conditional distribution $`\kappa _𝒥(\sigma ,\sigma ^{})`$ is the uniform perturbation metastate discussed above. Equivalent to the theorem’s conclusion is that for almost every $`𝒥`$, if $`\sigma `$ and $`\sigma ^{}`$ are chosen from $`\kappa _𝒥(\sigma ,\sigma ^{})`$, there is zero $`\kappa _𝒥`$-probability of an (infinite-volume) interface with zero density. This rules out $`(\sigma ^L,\sigma ^L)`$ interfaces with $`d_s<d`$ that are also pinned.

2) Although the theorem as formulated here addresses ground and excited states at $`T=0`$, it should be extendable to pure states at (low) $`T>0`$ by “pruning” small thermally induced droplets.

3) Although our construction used periodic boundary conditions (both for simplicity and because KM and PY used them), the theorem can be applied to other boundary conditions chosen independently of the couplings (as in , to which we refer the reader for details).

4) Although our focus here is on spin glasses, the theorem applies equally to many other systems, including disordered and homogeneous ferromagnets.

Application to KMPY excitations. An immediate application of the theorem is that KMPY excitations cannot be pinned, and so cannot give rise to (regionally congruent) ground states. To see that the theorem applies to the PY construction, note first that the periodic boundary condition clearly implies torus translation invariance for the distribution of $`(𝒥^L,\sigma ^L)`$. The Hamiltonian perturbation (Eq. (2) of ) is constructed from $`\sigma ^L`$ in a translation-covariant way, and thus the distribution of $`(𝒥^L,\sigma ^L,\sigma ^L)`$ is also invariant.

A priori, the situation for the KM construction is more subtle. Here if $`\sigma ^L`$ results from a fixed pair of spins at sites $`z_0,z_0^{}`$, translation-invariance would be lost. So instead, for fixed $`𝒥^L`$ and $`\sigma ^L`$ we regard $`\sigma ^L`$ as the random excitation resulting when $`z,z^{}`$ are chosen for each $`\mathrm{\Lambda }_L`$ from the uniform distribution on its sites; this restores translation-invariance. A first result is then that for a fraction one (as $`L\mathrm{}`$) of such choices of $`z,z^{}`$, the interface cannot be pinned. But there’s a second, more noteworthy result: it follows that a fixed $`z_0,z_0^{}`$ cannot yield a KMPY (or similar) excitation (which now would necessarily be pinned) at all, since if it did there would be a positive fraction of random $`(z,z^{})`$’s yielding the same excitation. Thus a fixed $`z_0,z_0^{}`$ must either yield a droplet excitation with volume $`o(L^d)`$ (e.g., bounded) as $`L\mathrm{}`$ or else one with $`d_s=d`$.

Regional Congruence. As discussed above, regionally congruent thermodynamic states  are those whose relative interfaces have zero density. Examples are the interface states in homogeneous and disordered Ising ferromagnets for sufficiently large $`d`$. But prior to  and  there had been no evidence for the existence of regional congruence in spin glasses.

In  and , we proved that thermodynamic states generated by coupling-independent boundary conditions (e.g., periodic, antiperiodic, free, fixed, etc.) are (with probability one) either the same (modulo a global spin flip) or else incongruent (so $`d_s=d`$). Thus regionally congruent states can only arise, if at all, through a sequence of $`\mathrm{\Lambda }_L`$’s with coupling-dependent boundary conditions; i.e., those that are conditioned on the $`𝒥^L`$’s. This is consistent with the above theorem, because such boundary conditions typically violate translation-invariance. E.g., interface states in ferromagnets arise through boundary conditions, such as Dobrushin  (i.e., plus boundary spins above the “equator” and minus below), that are not translation-covariant and implicitly use the knowledge that all couplings are positive. But in a spin glass, they are now gauge-equivalent to boundary conditions that are both translation-covariant and coupling-independent.

But now conclusions beyond those of  follow from the theorem, because in addition regional congruence cannot arise through any translation-covariant construction. A consequence is that any algorithm of the types currently known cannot yield regionally congruent states. The KM and PY procedures typify two main approaches in the search for new excitations and/or states, given that a priori the interface location is unknown: either a uniformly random sampling procedure (KM) or else a global perturbation (PY) (or similar coupling to a carefully chosen external field). Although the latter approach is coupling-dependent, it is also translation-covariant and so cannot generate regional congruence. One method that could find regional congruence, if it occurs, is an exhaustive search through all $`2^{O(L^{d1})}`$ fixed boundary conditions on $`\mathrm{\Lambda }_L`$ for each $`L`$; but that is hardly an option.

We conclude by noting that the order one energy difference of $`\sigma ,\sigma ^{}`$ plays no role in our theorem, which in fact has far wider applicability. Only one of its applications is to KMPY excitations, but these have generated interest because both $`\theta ^{}=0`$ and $`d_s<d`$. However, the order one energy does play a crucial role in our earlier heuristic argument , which we briefly summarize (and slightly extend) here. It was conjectured  and subsequently proved  (also ) that free energy fluctuations in spin glasses with coupling-independent boundary conditions on $`\mathrm{\Lambda }_L`$ scale as $`L^\theta `$, with $`\theta (d1)/2`$. It should then follow that the exponent $`\theta _r`$ characterizing interface energies for regionally congruent states satisfies the bound $`\theta _r>\theta `$ (as argued in ) because otherwise regional congruence could be seen with coupling-independent boundary conditions. So pinned interfaces with $`d_s<d`$ (hence regionally congruent states) must have $`\theta _r>\theta `$; if the latter inequality is violated, as in KMPY excitations, then the interface should deflect to infinity.

Acknowledgments. This research was supported in part by NSF Grants DMS-98-02310 (CMN) and DMS-98-02153 (DLS). The authors thank Peter Young and Matteo Palassini for useful correspondence. They also thank the Santa Fe Institute, where part of this work was done, for its hospitality.


text:
# Phenomenological Theory of Survival

## Abstract

Theoretical analysis proves that human survivability is dominated by an unusual physical, rather than biological, mechanism, which yields an exact law. The law agrees with all experimental data, but, contrary to existing theories, it is the same for an entire species, i.e., it is independent of the population, its phenotypes, environment and history. The law implies that the survivability changes with environment via phase transitions, which are simultaneous for all generations. They allow for a rapid (within few percent of the life span) and significant increase in the life expectancy even above its value at a much earlier age.

Mortality is one of the most universal and important phenomena in biology. Human mortality is extensively studied over two centuries  and is arguably the best statistically quantified biological phenomenon. There exist several evolutionary theories of aging  (the first one  is only 50 years old). Yet, C. Franceschi (in ref.) notes: ”Longevity is a trait with some peculiarities because of the unnecessary nature of aging”. This paper proves that, in contrast to all existing theories, human mortality is dominated by an unusual physical (rather than biological) mechanism. Presumably, the same mechanism dominates mortality of laboratory fly populations.

I study a large and diverse amount of over 150,000 data points from 1528 period life tables, for 3 races on 4 continents in 16 countries during over a century of their history, 36 cases total (here and on ”a case” denotes the population of a given sex in a given country at any time in its history; in the USA white, black and total populations are considered as separate cases). Each period life table presents in any given case the probability $`l_x`$ to survive to any (yearly) age x in a given calendar year. To escape any arbitrary scaling or adjustment, I use life table variables only.

The survivability $`l_x`$ depends on age $`x`$, calendar year $`t`$ and a complete set $`A`$ of all parameters, comprehensively describing a considered population, i.e. its phenotypes and their heterogeneous social, medical, dietary, etc. conditions during the previous $`x`$ years (from $`t`$ \- $`x`$ till $`t`$). So,$`l_x=l_x(t,A)`$. Mortality rate, and thus its statistics and accuracy, are low in early age. The number of survivors, and thus accuracy, decrease in old age. So, middle age survivability, e.g., $`l_{40}`$, is determined more accurately. Present $`l_x(t,A)`$ and $`l_{40}(t,A)`$ from the same life table as the ordinate and abscissa of a point. According to Fig. 1, for different races, countries, continents, histories $`l_x`$ is dominated by the function $`l_x(l_{40})`$, which depends only on $`l_{40}`$, i.e. invariant to $`A`$ and is therefore the same for an entire species. The deviation of $`l_x`$ from this function depends on $`A`$ (and thus on population phenotypes, living conditions, history). It is relatively small. From now on I denote $`l_x(l_{40})`$ as ”the invariant survivability” and, unless specifically stated otherwise, consider this (dominant) fraction of survivability only.

By Fig. 1, at any age and in every case, $`l_{40}`$ significantly changes with living conditions. Thus, since living conditions of different groups in a given case are different, these groups may have different probabilities $`l_x^{}`$ to survive to $`x`$ years (the value of $`l_x`$ in the population is the average $`l_x^{}`$; $`0l_x,l_x^{}1`$). Suppose the probability density of a given $`l_{40}^{}`$ in the population is c($`l_{40}`$,$`l_{40}^{}`$), i.e.

$$\underset{0}{\overset{1}{}}c(l_{40},l_{40}^{})𝑑l_{40}^{}=1,l_x(l_{40})=\underset{0}{\overset{1}{}}l_x^{}c(l_{40},l_{40}^{})𝑑l_{40}^{}$$
(1)
By Eq. 1, $`0\mathrm{min}(l_{40}^{})l_{40}\mathrm{max}l_{40}^{}`$.So, $`l_{40}^{}=0`$ when $`l_{40}`$ = 0 and $`l_{40}^{}=1`$ when $`l_{40}=1`$. Thus, invariance, which allows one to introduce c and which yields Eq. 1, implies that any population is homogeneous in $`l_{40}^{}`$ at the distribution boundaries. In virtue of invariance, the dependence of $`l_x^{}`$ on $`l_{40}^{}`$ is the same as the dependence $`l_x(l_{40})`$, i.e., $`l_x^{}=l_x(l_{40}^{}).`$ This, by Eq. 1, yields a remarkable symmetry of the invariant survivability to the transformations, specified by the function $`c`$ of two variables:

$$l_x\left[\underset{0}{\overset{1}{}}l_{40}^{}c(l_{40},l_{40}^{})𝑑l_{40}^{}\right]=\underset{0}{\overset{1}{}}l_x(l_{40}^{})c(l_{40},l_{40}^{})𝑑l_{40}^{}$$
(2)
Such invariance is sufficiently restrictive to mathematically accurately determine $`l_x(l_{40})`$ (as an invariant of survivability dynamics with respect to any population and environmental changes). Indeed, in a general case the function $`l_x(l_{40})`$, which yields Eq. 2, must be linear. A special case of $`c(l_{40},l_{40}^{})=c^{}\delta (l_{40}^{}l^{})+c^\mathrm{"}\delta (l_{40}^{}l^\mathrm{"})`$yields $`c^{}+c^\mathrm{"}=1`$ by Eq.1, and $`l_x(l_{40})=l_x(c^{}l^{}+c^\mathrm{"}l^\mathrm{"})=c^{}l_x(l^{})+c^\mathrm{"}l_x(l^\mathrm{"})`$ by Eq. 2. The latter equation implies $`dl_x/dl_{40}=dl_x/dl_x^{}=`$const. On the other hand, a linear $`l_x(l_{40})`$ yields Eq. 2 with any arbitrary $`c`$ in virtue of Eq. 1. Similarly, if an entire population is distributed in a certain interval of $`l_{40}`$, then $`l_x(l_{40})`$ is linear in this interval, and the population is homogeneous at its ends (as it was in the previous case). Correspondingly, if any population is homogeneous at certain $`l_{40}`$ points only, then an entire population is distributed within one of the intervals $`\lambda _s<l_{40}<\lambda _{s+1}`$ (s = 0,1 is its ordinal number):

$$\text{}\underset{\lambda _s}{\overset{\lambda _{s+1}}{}}c(l_{40},l_{40}^{})𝑑l_{40}^{}=1;l_{40}=\underset{\lambda _s}{\overset{\lambda _{s+1}}{}}l_{40}^{}c(l_{40},l_{40}^{})𝑑l_{40}^{}$$
(3)
and $`l_x(l_{40})`$ is linear within these intervals. When $`l_{40}=\lambda _{s+1}`$, an entire population of all ages moves to the next interval. Clearly, such piecewise $`l_x(l_{40})`$ satisfies Eq.2. When $`\lambda _0=0,\lambda _1=1,`$ this case reduces to the previous one. If an infinite number of edge points in Eq.3 condenses and forms a continuous interval (where any population is homogeneous), then the invariant $`l_x(l_{40})`$ in this interval remains undetermined. Of course, the possibility to distinguish this case from the case of a large number of short linear segments is limited by the disregarded non-invariant contributions to survivability, but demographic studies prove that survivability is heterogeneous and exclude the former case (of a population which is homogeneous in a finite interval). Thus, in a general case

$$l_x=R_x^{(s)}l_{40}+a_x^{(s)},\lambda _sl_{40}\lambda _{s+1}.$$
(4)
So, except for the number of linear segments, the very existence of the invariant survivability allows one to establish its exact law, without any experiments, approximations, and assumptions.
Equation (4) is the implication of invariance only. Such invariance was also demonstrated  for medfly and fruitfly families whose different populations were extensively studied in different conditions. Fly statistics is rather low, thus their $`l_x`$ was studied as a function of the life expectancy at birth (which averages $`l_x`$ over different ages and is therefore more smooth). The study was based on fly populations hatched the same day. However, since in most cases fruitflies were kept in stationary conditions, it suggests that the invariant survivability law is general for humans and flies.

Piecewise linear law (4) and its age independent invariant intersections are the main predictions of this paper. They are verified with all experimental data - see, e.g., $`l_1`$, $`l_{60}`$,$`l_{80}`$ vs $`l_{40}`$ in Fig.2. Piecewise linear law agrees with (but has never been suggested in) demographic approximations . Slope jumps in Eq.(4) and Fig.2 are consistent with ref . However, in contrast to the qualitative observations in ref., accurate Eq.(4) yield quantitative predictions. Since all $`l_x`$ belong to the same calendar year, their slope jumps are simultaneous for all generations. Indeed, although the intersections in Fig.2 are reached at different times in different countries for different sexes, they are simultaneous (to the invariance accuracy) at different ages, and (in $`l_{40}`$) for different populations. The invariant law allows for the prediction of intersections. For instance, the extrapolation of $`l_{80}`$ beyond $`\lambda _1=0.95`$ in Fig. 1 yields $`l_{40}=1`$ (while $`l_{80}`$ is finite), i.e. implies no deceased until $`40y`$. This is hardly possible, and suggests a crossover, which is indeed seen in Fig.2 at $`\lambda _2=0.97`$ The extrapolations of $`l_{60}`$ and $`l_{80}`$ to $`l_{40}<0.94`$ yield $`l_{60}=0`$ and $`l_{80}=0`$ , i.e., no survivors beyond $`60y`$ (when $`l_{40}=0.27`$) and $`80y`$ (when $`l_{40}=0.42`$). This is never true, and suggests an intersection at $`0.4<l_{40}<0.7`$.

To elucidate the nature of the invariant law, present Eq. (4) in a different form:

$$\text{}l_x=c_xl_x^{(s)}+(1c_s)l_x^{(s+1)}whenl_x^{(s)}l_xl_x^{(s+1)}$$
(5)
Here $`l_{40}^{(s)}=\lambda _s;l_x^{(s)}=l_x(\lambda _s)`$ from Eq.(4); $`c_s`$ reduces to $`l_x^{(s)}`$ and $`l_{40}`$ \[thus, by Eq. (3), to c\]. Equation (5) accurately separates ”nature” and ”nurture” in survival. ”Nature” reduces to the fixed set of the intersection survivabilities . The set depends only on age and is invariant, i.e., independent of phenotypes (and thus of a specific DNA sequence at least in an entire species), their living conditions and life history. The dependence of the set on age is not determined by invariance, but invariance implies that it is the same for at least an entire species. (Moreover, it scales onto the same functions for species as remote as humans and flies ). Thus, it does not change at least as long as the species does not evolve into a different species (demographic data in Fig.2 verify it for 100-150 years). So, the set must be inheritable. Such set reminds of the body temperature (which in any living conditions is the same with few percent accuracy) of an entire class of birds and a subclass of placental mammals. Presumably, both the set and the body temperature are genetically determined, and independent of a specific DNA sequence. But the set, unlike the body temperature, strongly depends on age. ”Nurture” distributes the survivability $`l_x`$ at a given age between two adjacent intersection survivabilities. The concentration $`c_s(0c_s1)`$ is age independent (and may be related to, e.g., $`l_{40}`$). Thus, survivability follows environment (in particular, a new intersection survivability emerges) simultaneously for all generations.

Equation (5) relates $`c_s`$, and thus $`l_x,`$ to $`l_1`$. Since $`l_x`$ is the survival probability in the same calendar year as $`l_1`$, so, by Eq. (5), the survivability accurately and rapidly follows the change in environment according to the value of $`l_1`$ (i.e., the infant mortality $`q_0=1l_1`$), which is established in less than two years. Indeed, whatever the difference in environmental factors is, close values of infant mortality imply very close fractions of deceased at any given age in 1885 Swedish and 1947 Japanese females, despite of their different races, continents, countries, and 62 year gap in their different history. Since $`l_1`$ depends on the environment in a given year only, the invariant survivability $`l_x`$ is reversible (an entire survival curve comes back when $`q_0`$ changes non-monotonically and returns to its previous value), and statistically independent of the life history during its $`x`$ years, despite country specific, highly and irregularly changing, non-monotonic and non-stationary living conditions. World wars and epidemics, e.g., flu in 1918 Europe, significantly decrease $`l_x`$. (For instance, in 1915 the probability for a French male to survive to $`80y`$ was 5 times less than in 1913, twice less than in 1917, and 3.5 times less than for a 1861 Swedish male). Yet, they just slightly shift the plots (mostly vertically, and relatively little) in Fig.2. In a couple of years (which estimate the relaxation time at few percents of the life span) all plots restore their invariant dependence, i.e., the memory of the previous life history is erased. Accurate reversibility of the invariant survivability does not decrease with aging, even in old age. Unless such reversibility is related to some perfect biological rehabilitation (which is hardly possible), it implies an adiabatic change in a certain thermodynamic equilibrium. This is consistent with its relaxation time (rapid compared to the life span, but enormous on a microscopic scale). Equation (5) reduces $`l_x`$ to the fixed set of $`l_x^{(s)}`$. Since $`l_x^{(s)}`$ reversibly change into each other, they are related to different equilibrium thermodynamic states of the same system, i.e. to different phases. Thus, Eq. (5) relates the invariant survivability to a certain phase equilibrium, and the jumps in Fig.1 slopes to the emergence of a new phase. The phase concentration $`c_s`$, which Eq.(5) reduces to, e.g., $`l_1`$, is independent of age. This suggests that the age dependence of a $`l_x^{(s)}`$ is related to the difference in survivabilities provided by the same phase at different ages. Quasi-equilibrium phases in a living homeostatic being might be related to a meso- (e.g., a cell) or microscopic (e.g., DNA configuration) scale.

While the existence of survivability phases is accurately proven, their microscopic nature, as well as that of the parameter which determines their concentrations, remains unknown (as ”units of heredity” were to Mendel). However, any phase equilibrium may be reversibly manipulated. This means that human life expectancy may be rapidly (within few years) reversed to its value at a much earlier age. (Note that in just 8 years from 1947 till 1955 the life expectancy of Japanese females increased 26% at birth, 15% at $`60y`$, and 20% at $`80y`$).

Presumably, the change in phase concentrations may affect other than survivability characteristics, in particular, aging and disease resistance. Indeed, Eq.(5) is not violated even in Japan prior to 1949 (see Fig.1) and in Finland from 1890 till 1940 (see Fig.2), although their mortality has a strongly tubercular age pattern during this period.

The suggested phase equilibrium nature of survivability implies that non-equilibrium (in particular, sufficiently non-stationary and heterogeneous conditions which depend on age, sex, social mobility, immigration and other factors), may lead to more than two adjacent phases and to hysteresis in the adjustment. The contribution of non-equilibrium phases is most pronounced in old age, when the difference in $`l_x^{(s)}`$ is significantly higher. This agrees with Fig.2. Hysteresis yields small jumps in $`l_x`$ and shifts of the intersections, which are indeed present in certain cases. The accuracy of the invariant law estimates maximal concentrations of ”extra” phases. They are lower for females, presumably suggesting lower heterogeneity of female populations.

Survivability independence of life history implies no correlation between early and old age invariant mortalities in a given cohort (born the same year). This disagrees with evolutionary theories of aging . Stochastic mutation accumulation  theory is inconsistent with reversibility and rapid accurate survivability change with environment. Optimal allocation of metabolic resources (pleiotropy and disposable soma theories) implies strong correlation between survival in young and old age in populations born the same year. This is inconsistent with the survivability independence of life history. Reversibility is also inconsistent with mortality theories (thelomers, oxygen consumption, free radicals, somatic mutations). Although natural mortality in the wild is mostly due to extrinsic hazards, invariant mortality, which dominates in species as remote as humans and flies (in certain protected environment), and which rapidly, accurately and simultaneously in all generations changes with environment, calls for biological, evolutionary and microscopic physical theories.
These conclusions, as well as an unusual mechanism of mortality, are accurate implications of the exact law, which agrees with all demographic approximations and studies . For instance, the infant mortality is widely appreciated by demographers as a sensitive barometer of environmental conditions. However, demographic approximations are developed primarily as a useful tool of maximally accurate estimation and forecast. They are often country and time specific, and approximate specific demographic data better than Eq. (4), but they do not consider nor care about the underlying general law and its mechanism. In contrast, my goal is the exact law, albeit of the invariant mortality only (to discover inertia, one must disregard friction!) and its mechanism, which yield new biological insights.

This paper was stimulated by enlightening discussions with Profs. Y. Aharonov, S.M. Jazwinski, S. Horiuchi, D.S. Thaler, A. Libchaber,and B. Tsirelson. I am grateful to Prof. S. Horiuchi, S. Smallwood, A. Hanika, K. Anderson, S. Helgadottir, C. Capocci, E. Waeyetens, N. White, U. Wegmueller, V. Kannisto, M. Nieminen, K. Deaves for life tables; to I. Kolodnaya for crucial technical assistance; A. Kobi for computer simulations, J. & R. Meyerhoff chair for financial support.


text:
# Chapter 1 Binaries and Globular ClusterDynamics

## Chapter 1 Binaries and Globular Cluster <br>Dynamics

Dept of Physics, MIT, Cambridge, MA 02139, USA

### 1 Introduction

The dynamical evolution of dense star clusters is a problem of fundamental importance in theoretical astrophysics, but many aspects of the problem have remained unresolved in spite of years of numerical work and improved observational data (see Meylan & Heggie 1997 for a recent review). The realization over the last 10 years that primordial binaries are present in globular clusters in dynamically significant numbers has completely changed our theoretical perspective on these systems (see, e.g., Gao et al. 1991; Hut et al. 1992; Sigurdsson & Phinney 1995). Most importantly, dynamical interactions between hard primordial binaries and other single stars or binaries are now thought to be the primary mechanism for supporting a globular cluster against gravothermal contraction and avoiding core collapse. In addition, exchange interactions between primordial binaries and compact objects can explain very naturally the formation of large numbers of X-ray binaries and recycled pulsars in globular cluster cores (see, e.g., Camilo et al. 2000). Dynamical interactions involving primordial binaries can also result in dramatically increased collision rates in globular clusters. This is because the interactions are often resonant, with all the stars involved remaining together in a small volume for a long time ($`10^210^3`$ orbital times). For example, in the case of an interaction between two typical hard binaries with semi-major axes $`1`$AU containing $`1M_{}`$ main-sequence stars, the effective cross section for a direct collision between any two of the four stars involved is essentially equal to the entire geometric cross section of the binaries (Bacon, Sigurdsson & Davies 1996; Cheung, Portegies Zwart, & Rasio 2001). This implies a collision rate $`100`$ times larger than for single stars. Direct observational evidence for stellar collisions and mergers in globular clusters comes from the detection of large numbers of blue stragglers concentrated in the dense cluster cores (see, e.g., Bailyn 1995).

### 2 Monte Carlo Simulations of Cluster Dynamics

The first Monte Carlo methods for calculating the dynamical evolution of star clusters in the Fokker-Planck approximation were developed more than 30 years ago. They were first used to study the development of the gravothermal instability (Hénon 1971a,b; Spitzer & Hart 1971a,b). More recent implementations have established the Monte Carlo method as an important alternative to direct $`N`$-body integrations (see Spitzer 1987 for an overview). The main motivation for our recent work at MIT was our realization a few years ago that the latest generation of parallel supercomputers now make it possible to perform Monte Carlo simulations for a number of objects equal to the actual number of stars in a globular cluster (in contrast, earlier work was limited to using a small number of representative “superstars,” and was often plagued by high levels of numerical noise). Therefore, the Monte Carlo method allows us to do right now what remains an elusive goal for $`N`$-body simulations (see, e.g., Aarseth 1999): perform realistic, star-by-star computer simulations of globular cluster evolution. Using the correct number of stars in a dynamical simulation ensures that the relative rates of different dynamical processes (which all scale differently with the number of stars) are correct. This is particularly crucial if many different dynamical processes are to be incorporated, as must be done in realistic simulations.

Our implementation of the Monte Carlo method is described in detail in the papers by Joshi, Rasio, & Portegies Zwart (2000), Joshi, Nave, & Rasio (2001), and Joshi, Portegies Zwart, & Rasio (2001). We adopt the usual assumptions of spherical symmetry (with a 2D phase space distribution function $`f(E,J)`$, i.e., we do not assume isotropy) and standard two-body relaxation in the weak scattering limit (Fokker-Planck approximation). In its simplest version, our code computes the dynamical evolution of a self-gravitating spherical cluster of $`N`$ point masses whose orbits in the cluster are specified by an energy $`E`$ and angular momentum $`J`$, with perturbations $`\mathrm{\Delta }E`$ and $`\mathrm{\Delta }J`$ evaluated on a timestep that is a fraction of the local two-body relaxation time. The cluster is assumed to remain always very close to dynamical equilibrium (i.e., the relaxation time must remain much longer than the dynamical time). Our main improvements over Hénon’s original method are the parallelization of the basic algorithm and the development of a more sophisticated method for determining the timesteps and for computing the two-body relaxation from representative encounters between neighboring stars. Our new method allows the timesteps to be made much smaller in order to resolve the dynamics in the cluster core more accurately.

We have performed a large number of test calculations and comparisons with direct $`N`$-body integrations, as well as direct integrations of the Fokker-Planck equation in phase space, to establish the accuracy of our basic treatment of two-body relaxation (Joshi et al. 2000). Fig. 1 shows the results from a typical comparison between Monte Carlo and $`N`$-body simulations. Fig. 2 shows the results obtained with our code for Heggie’s Collaborative Experiment run (Heggie et al. 1998).

### 3 Summary of Recent Results

Our recent work has focused on the addition of more realistic stellar and binary processes to the basic Monte Carlo code, as well as a simple but accurate implementation of a static tidal boundary in the Galactic field (Joshi et al. 2001a). As a first application, we have studied the dependence on initial conditions of globular cluster lifetimes in the Galactic environment. As in previous Fokker-Planck studies (Chernoff & Weinberg 1990; Takahashi & Portegies Zwart 1998), we include the effects of a power-law initial mass function (IMF), mass loss through a tidal boundary, and single star evolution, and we consider initial King models with varying central concentrations. We find that the disruption and core-collapse times of our models are significantly longer than those obtained with previous 1D (isotropic) Fokker-Planck calculations, but agree well with more recent results from direct $`N`$-body simulations and 2D Fokker-Planck integrations. In agreement with previous studies, our results show that the direct mass loss due to stellar evolution causes most clusters with a low initial central concentration to disrupt quickly in the Galactic tidal field. The disruption is particularly rapid for clusters with a relatively flat IMF. Only clusters born with high central concentrations or with very steep IMFs are likely to survive to the present and undergo core collapse.

In another recent study, we have used our Monte Carlo code to examine the development of the Spitzer “mass stratification instability” in simple two-component clusters (Watters, Joshi, & Rasio 2000). We have performed a large number of dynamical simulations for star clusters containing two stellar populations with individual masses $`m_1`$ and $`m_2>m_1`$, and total masses $`M_1`$ and $`M_2<M_1`$. We use both King and Plummer model initial conditions and we perform simulations for a wide range of individual and total mass ratios, $`m_2/m_1`$ and $`M_2/M_1`$, in order to determine the precise location of the stability boundary in this 2D parameter space. As predicted originally by Spitzer (1969) using simple analytic arguments, we find that unstable systems never reach energy equipartition, and are driven to rapid core collapse by the heavier component. These results have important implications for the dynamical evolution of any population of primordial black holes or neutron stars in globular clusters. In particular, primordial black holes with $`m_2/m_110`$ are expected to undergo very rapid core collapse independent of the background cluster, and to be ejected from the cluster through dynamical interactions between single and binary black holes (see Portegies Zwart & McMillan 2000 and references therein). We have also used Monte Carlo simulations of simple two-component systems to study the evaporation (or retention) of low-mass objects in globular clusters, motivated by the surprising recent observations of planets and brown dwarfs in several clusters (Fregeau et al. 2001).

Much of our current work concerns the treatment of dynamical interactions with primordial binaries. We are in the process of completing a first study of globular cluster evolution with primordial binaries (Joshi et al. 2001b), based on the same set of approximate cross sections and recipes for dynamical interactions used in the Fokker-Planck simulations of Gao et al. (1991). Typical results are illustrated in Figs. 3–6. The heating of the cluster core generated by a small population of primordial binaries can support the cluster against core collapse for very long times, although the details of the evolution depend sensitively both on the number of stars $`N`$ and on the initial binary fraction $`f_\mathrm{b}`$. Clusters with smaller $`N`$ tend to evolve faster towards deeper core collapse (compare Figs. 3 and 4, and note that the main dependence on $`N`$ through $`t_{\mathrm{r}h}N/\mathrm{log}N`$ has been scaled out), and so do clusters with smaller initial binary fractions (Fig. 5). After the initial core collapse, gravothermal oscillations powered by primordial binaries can continue for very long times $`50100t_{\mathrm{r}h}`$ even with an initial binary fraction as small as 1%. With an initial binary fraction of 10%, we observe a single, very moderate phase of core collapse during which the core shrinks by a factor $`510`$ in radius and then re-expands rapidly on a timescale of a few central relaxation times (this occurs at $`t/t_{\mathrm{r}h}60`$ in Fig. 3 and at $`t/t_{\mathrm{r}h}20`$ in Fig. 4). The cluster then enters a second quasi-equilibrium phase of primordial binary burning with the core radius increasing slowly until well beyond $`100t_{\mathrm{r}h}`$.

Since $`t_{\mathrm{r}h}10^9`$yr for most clusters, if the type of evolution illustrated in Fig. 3 applied to all globular clusters, there would be no “core-collapsed” clusters in the Galaxy (10–15% of all Galactic globular clusters are classified observationally as “core collapsed”). However, the timescale on which real clusters will exhaust their primordial binary supply and undergo (deeper) core collapse depends on a number of factors not considered here, including the orbit of the cluster in the Galaxy (the simulations are for an isolated cluster, but mass loss and tidal shocking can accelerate the evolution dramatically), and the stellar IMF (the clusters shown here contain all equal-mass stars and binaries; a more realistic mass spectrum will also accelerate the evolution). In addition, some clusters may be formed with much fewer primordial binaries than even the 1% considered in Fig. 5. However, the simple picture that emerges from these simulations may well, to first approximation, describe the dynamical state of most Galactic globular clusters observed today. Indeed, for a cluster in the stable “binary burning” phase, the ratio of half-mass radius to core radius $`r_\mathrm{h}/r_\mathrm{c}210`$ (Fig. 6), which is precisely the range of values observed for the $`80\%`$ of globular clusters that have a well-resolved core and are well-fitted by King models (see, e.g., Djorgovski 1993). Some of these clusters may have gone in the past through a brief episode of “moderate core collapse” (as shown around $`t60t_{\mathrm{r}h}`$ in Fig. 3). Yet, they should not be called “core-collapsed” (nor would they be classified as such by observers). Unfortunately some theorists will even call “core-collapsed” clusters that have just reached the initial phase of binary burning ($`t1050t_{\mathrm{r}h}`$ in Fig. 3). Since the core has just barely contracted by a factor $`23`$ by the time it reaches this phase, it seems hardly justified to speak of a “collapsed” state.

The addition of binary stellar evolution processes to our simulations will allow us to study in detail the dynamical formation mechanisms for many exotic objects such as X-ray binaries and millisecond radio pulsars, which have been detected in large numbers in globular clusters. For example, exchange interactions between neutron stars and primordial binaries can lead to common-envelope systems and the formation of short-period neutron-star / white-dwarf binaries that can become visible both as ultracompact X-ray binaries and binary millisecond pulsars with low-mass companions (see, e.g., Camilo et al. 2000, on observations of 20 such millisecond radio pulsars in 47 Tuc). Rasio, Pfahl, & Rappaport (2000) present a preliminary study of this formation scenario, based on simplified dynamical Monte Carlo simulations. In a dense cluster such as 47 Tuc, most neutron stars acquire binary companions through exchange interactions with primordial binaries. The resulting systems have semimajor axes in the range $`0.11`$AU and neutron star companion masses $`13M_{}`$. For many of these systems it is found that, when the companion evolves off the main sequence and fills its Roche lobe, the subsequent mass transfer is dynamically unstable. This leads to a common envelope phase and the formation of short-period neutron-star / white-dwarf binaries. For a significant fraction of these binaries, the decay of the orbit due to gravitational radiation will be followed by a period of stable mass transfer driven by a combination of gravitational radiation and tidal heating of the companion. The properties of the resulting short-period binaries match well those of observed binary pulsars in 47 Tuc (Fig. 7). A similar dynamical scenario involving massive CO white dwarfs in place of neutron stars could explain the recent detections of double degenerate binaries containing He white dwarfs concentrated in the core of a dense globular cluster (Hansen et al. 2001; see Edmonds et al. 1999 and Taylor et al. 2001 on observations of He white dwarfs in the core of NGC 6397).

We are also currently working on incorporating into our Monte Carlo simulations a more realistic treatment of tidal interactions, and, in particular, tidal shocking through the Galactic disk (based on Gnedin, Lee, & Ostriker 1999). Tidal shocks can accelerate significantly both core collapse and the evaporation of globular clusters, reducing their lifetimes in the Galaxy (Gnedin & Ostriker 1997). Future work may include a fully dynamical treatment of all strong binary-single and binary-binary interactions (exploiting the parallelism of the code to perform separate numerical 3- or 4-body integrations for all dynamical interactions) as well as a fully dynamical treatment of tidal shocking (performing short, $`N`$-body integrations for each passage of the cluster through the Galactic disk or near the bulge).

### Acknowledgments

This work was supported in part by NSF Grant AST-9618116 and NASA ATP Grant NAG5-8460. Our computations were performed on the Cray/SGI Origin2000 supercomputer at Boston University under NCSA Grant AST970022N.

### References