text:
# Center of quantum group in roots of unity and the restriction of integrable models The talk given in RAQIS03 conference, Aneccy, 25-28/03/03

## 1 Introduction

F. Alcaraz et al.  discovered a remarkable fact: the $`XXZ`$ model with the special open boundary conditions (OBC) and a rational value of the anisotropy parameter admits a restriction. The model arising as a result of the restriction coincides in the thermodynamic limit with one of the Minimal Models of CFT. The algebraic reason for the restriction was explained in  and . It was shown in  that the $`XXZ`$ model with the OBC considered in  has not only integrability but also $`U_q(sl(2))`$ symmetry. In roots of unity, the state space of the model decomposes to the sum of “good” and “bad” representations of $`U_q(sl(2))`$. The restriction of Alcaraz et al. is equivalent throwing out “bad” parts and keeping only the highest vectors of “good” representations. In , a new monodromy matrix was constructed that is bilinear in terms of generators of the quantum group $`A(u)`$, $`B(u)`$, $`C(u)`$, $`D(u)`$ connected in the usual way with the Hamiltonian of the model in  and compatible with integrability. The twisted trace of this monodromy matrix (Sklyanin transfer matrix) also has $`U_q(sl(2))`$ symmetry and admits the restriction.

In this paper, we generalize the Pasqier–Saleur construction. We show that not only $`U_q(sl(2))`$-invariant but also a much wider class of OBC indeed admits the restriction (the Sklyanin construction of the transfer matrix also works for this wider class of OBC). The decisive condition for compatibility of the Hamiltonian, the transfer matrix, and other operators with the restiction is their “weak” commutativity with a special element of the quantum group. The notion of “weak” commutativity and its connection with the extension of the center of the quantum group in roots of unity is explained below.

In Sec. 2, we describe the conditions that the operators must have in order to admit the restriction. In Sec. 3, some such operators are found. We discuss some possible generalizations of the construction in the last section.

## 2 Center of the quantum group and the restriction in roots of unity

As usual, let the $`R`$-matrix $`R(u)`$ denote the solution of the Yang–Baxter equation. We consider the simplest and well-known $`R(u)`$ matrix corresponding to the six-vertex model, whose elements can be written as

$`R_{\alpha \alpha }^{\alpha \alpha }(u)=\rho \mathrm{sin}(u+\eta ),`$ (1)
$`R_{\alpha \beta }^{\alpha \beta }(u)=\rho \mathrm{sin}u,`$
$`R_{\beta \alpha }^{\alpha \beta }(u)=\rho \mathrm{sin}\eta ,`$
where $`\alpha ,\beta =1,2`$, $`\alpha \beta `$, and $`\eta `$ is the so-called anisotropy parameter.

The quantum group A connected with $`R(u)`$ is generated by $`A(u)`$, $`B(u)`$, $`C(u)`$, $`D(u)`$, entries of the monodromy matrix $`L(u)`$, which satisfies

$$R_{12}(uv)L_1(u)L_2(v)=L_2(v)L_1(u)R_{12}(uv).$$
(2)
As shown by V. Tarasov , the center of A in roots of unity, i.e., $`\eta =\pi m/N`$, is generated by the following elements of A:

$`A(u)=A(u)A(u+\eta )\mathrm{}A(u+(N1)\eta ),`$
$`B(u)=B(u)B(u+\eta )\mathrm{}B(u+(N1)\eta ),`$
$`C(u)=C(u)C(u+\eta )\mathrm{}C(u+(N1)\eta ),`$
$`D(u)=D(u)D(u+\eta )\mathrm{}D(u+(N1)\eta ).`$
For convenience, we let $`B(u)`$ denote the central element,

$$B(u)=B(u)B_1(u),$$

where

$$B_1(u)=B(u+\eta )\mathrm{}B(u+(N1)\eta ).$$

We now fix $`V=C^2\mathrm{}C^2`$ as the representation space of our quantum group. It is easy to see that for arbitrary $`v`$,

$$B(v)=B(v)B_1(v)=0$$

on this space. We can then define the state space of a restricted model as

$$W(v)=\text{Ker}B(v)/\text{Im}B_1(v).$$

In the limit $`v\mathrm{}`$, $`B(v)`$ coincides up to a scalar factor with $`X`$, one of the generators of $`U_q(sl(2))`$. As a result, $`W(\mathrm{})`$ coincides with the space of the “good” highest vectors of Pasqier–Saleur. It was shown in  that the Hamiltonian of the $`XXZ`$ chain with OBC of special type

$$H_{XXZ}=\underset{n=1}{\overset{L1}{}}\left[\sigma _n^+\sigma _{n+1}^{}+\sigma _n^{}\sigma _{n+1}^++\frac{\mathrm{cos}\eta }{2}\sigma _n^z\sigma _{n+1}^z+i\frac{\mathrm{sin}\eta }{2}(\sigma _n^z\sigma _{n+1}^z)\right]$$

is invariant under the quantum algebra $`U_q(sl(2))`$. Here, $`q=e^\eta `$. Because of this, $`H_{XXZ}`$ is properly defined on $`W(\mathrm{})=\text{Ker}X/\text{Im}X^{(N1)}`$.

In the thermodynamic limit, where $`L\mathrm{}`$, the spectrum of low-lying states coincides (in Cardy’s sense) with $`M(N1/N)`$, one of the Minimal Models of CFT. In the next section, we show that the construction in  can be generalized to arbitrary values of the parameter $`v`$. The corresponding Hamiltonian is

$`H_{XXZ}=`$ $`{\displaystyle \underset{n=1}{\overset{L1}{}}}\left[\sigma _n^+\sigma _{n+1}^{}+\sigma _n^{}\sigma _{n+1}^++{\displaystyle \frac{\mathrm{cos}\eta }{2}}\sigma _n^z\sigma _{n+1}^z\right]`$
$`+i{\displaystyle \frac{\mathrm{sin}\eta }{2}}(\mathrm{cot}v\sigma _1^z\mathrm{cot}(v+\eta )\sigma _L^z).`$ (3)
We now find the sufficient conditions for any operator $`Q`$ to be projectible on $`W(v)`$. It is easy to see that sufficient conditions are that there exist some operators $`\widehat{Q}`$ and $`\widehat{Q}_1`$ for a given $`Q`$ such that

$$B(v)Q=\widehat{Q}B(v)$$
(4)
and

$$QB_1(v)=B_1(v)\widehat{Q}_1.$$
(5)
Indeed, Eq. (4) guarantees that if a vector $`\psi \text{Ker}B(v)`$, then the vector $`Q\psi \text{Ker}B(v)`$. It follows from Eq. (5) that if the difference of two vectors $`\psi _1`$ and $`\psi _2`$ belong to $`\text{Im}B_1(v)`$, i.e., if $`\psi _1\psi _2=B_1(v)\chi `$, then the difference of $`Q\psi _1`$ and $`Q\psi _2`$ also belongs to $`\text{Im}B_1(v)`$.

## 3 Sklyanin transfer matrix and other operators respecting the restriction in roots of unity

In , E. Sklyanin explained the integrability of the $`XXZ`$ model with OBC of the form

$`H_{XXZ}=`$ $`{\displaystyle \underset{n=1}{\overset{L1}{}}}\left[\sigma _n^+\sigma _{n+1}^{}+\sigma _n^{}\sigma _{n+1}^++{\displaystyle \frac{\mathrm{cos}\eta }{2}}\sigma _n^z\sigma _{n+1}^z\right]`$
$`+i{\displaystyle \frac{\mathrm{sin}\eta }{2}}(\mathrm{cot}(\xi _++\eta /2)\sigma _1^z+\mathrm{cot}(\xi _{}\eta /2)\sigma _L^z).`$
by constructing a special monodromy matrix and using it to diagonalize Hamiltonian (3) and the corresponding transfer matrix by means of the algebraic Bethe ansatz.

Let $`K_+(u)=K(u+\eta /2,\xi _+)`$ and $`K_{}(u)=K(u\eta /2,\xi _{})`$, where

$$K(u,\xi )=\left[\begin{array}{cc}\mathrm{sin}(u+\xi )& 0\\ 0& \mathrm{sin}(u\xi )\end{array}\right].$$

Then $`K_\pm `$ satisfies the boundary Yang–Baxter equations .

The Sklyanin monodromy matrix is defined  as

$`\mathrm{\Lambda }(u)`$ $`=\sigma ^2L^t(u)\sigma ^2K_+(u)L(u)`$
$`=\left[\begin{array}{cc}\mathrm{\Lambda }_1^1& \mathrm{\Lambda }_1^2\\ \mathrm{\Lambda }_2^1& \mathrm{\Lambda }_2^2\end{array}\right].`$
Using (2) and the boundary Yang–Baxter equation, we can prove that $`\mathrm{\Lambda }(u)`$ satisfies the same equation as $`K_+(u)`$. It gives the commutation relations between $`\mathrm{\Lambda }_j^i`$. The Sklyanin transfer matrix is defined as

$`T_\text{S}(u)`$ $`=\text{tr}\mathrm{\Lambda }(u)K_{}(u)`$
$`=\mathrm{sin}(u\eta /2+\xi _{})\mathrm{\Lambda }_1^1\mathrm{sin}(u\eta /2\xi _{})\mathrm{\Lambda }_2^2.`$
Explicit expressions for the elements of the Sklyanin monodromy matrix are

$`\mathrm{\Lambda }_1^1(u)=\mathrm{sin}(u+\eta /2+\xi _+)A(u)D(u)+\mathrm{sin}(u+\eta /2\xi _+)C(u)B(u),`$
$`\mathrm{\Lambda }_2^2(u)=\mathrm{sin}(u+\eta /2+\xi _+)B(u)C(u)\mathrm{sin}(u+\eta /2\xi _+)D(u)A(u),`$
$`\mathrm{\Lambda }_1^2(u)=\mathrm{sin}(u+\eta /2+\xi _+)B(u)D(u)+\mathrm{sin}(u+\eta /2\xi _+)D(u)B(u),`$
$`\mathrm{\Lambda }_2^1(u)=\mathrm{sin}(u+\eta /2+\xi _+)A(u)C(u)\mathrm{sin}(u+\eta /2\xi _+)C(u)A(u).`$
Let $`\xi _+=v\eta /2`$ and $`\xi _{}=v\eta /2`$. Then the following relations are satisfied:

$$B(v)T_\text{S}(u)=\widehat{T}_\text{S}(u)B(v),$$
(6)
where the explicit expression for the Sklyanin transfer matrix for these values $`\xi _{pm}`$ is

$`T_\text{S}(u)=`$ $`\mathrm{sin}(uv\eta )\mathrm{sin}(u+v)A(u)D(u)`$
$`+\mathrm{sin}(uv\eta )\mathrm{sin}(uv+\eta )C(u)B(u)`$
$`+\mathrm{sin}(u+v)\mathrm{sin}(u+v)B(u)C(u)`$
$`+\mathrm{sin}(u+v)\mathrm{sin}(uv+\eta )D(u)A(u)`$ (7)
and

$`\widehat{T}_\text{S}(u)=`$ $`\mathrm{sin}(uv)\mathrm{sin}(u+v+\eta )A(u)D(u)`$
$`+\mathrm{sin}(uv)\mathrm{sin}(uv)C(u)B(u)`$
$`+\mathrm{sin}(u+v+\eta )\mathrm{sin}(u+v\eta )B(u)C(u)`$
$`+\mathrm{sin}(u+v+\eta )\mathrm{sin}(uv)D(u)A(u).`$
We have

$$B(v)\mathrm{\Lambda }_1^2(u)=\widehat{\mathrm{\Lambda }}_1^2(u)B(v),$$
(8)
where explicitly

$`\mathrm{\Lambda }_1^2(u)=\mathrm{sin}(u+v)B(u)D(u)+\mathrm{sin}(uv+\eta )D(u)B(u),`$
$`\widehat{\mathrm{\Lambda }}_1^2(u)=\mathrm{sin}(u+v+\eta )B(u)D(u)+\mathrm{sin}(uv)D(u)B(u).`$
Equations (6) and (8) were verified by direct calculation.

The operators $`T_\text{S}(u)`$ and $`\mathrm{\Lambda }_1^2(u)`$ thus satisfy the first condition, Eq. (4), for the restriction. We conjecture that they also satisfy the second one, Eq. (5), if $`q`$ is a root of unity. This conjecture was explicitly verified by direct calculation for third-degree roots for the case $`T_\text{S}(u)`$ (but not for the case $`\mathrm{\Lambda }_1^2(u)`$). It would nice to find an elegant general proof.

The operators $`T_\text{S}(u)`$ and $`\mathrm{\Lambda }_1^2(u)`$ depend on one parameter. There also exists a two-parameter family of operators satisfying (4) and (5).

By definition, let

$`T(x_{ij};u_1,u_2)=`$ $`x_{11}A(u_1)D(u_2)+x_{22}D(u_1)A(u_2)`$
$`+x_{12}B(u_1)C(u_2)+x_{21}C(u_1)B(u_2).`$
Then

$$B(v)T(x_{ij}(v);u_1,u_2)=T(\widehat{x}_{ij}(v);u_1,u_2)B(v),$$
(9)
where $`T(x_{ij}(v);u_1,u_2)`$ and $`T(\widehat{x}_{ij}(v);u_1,u_2)`$ are obtained from $`T(x_{ij};u_1,u_2)`$ by suitably substituting $`x_{ij}(v)`$ and $`\widehat{x}_{ij}(v)`$ for $`x_{ij}`$ and

$`x_{11}(v)=\mathrm{sin}(u_1\eta v)\mathrm{sin}(u_2v),`$
$`x_{22}(v)=\mathrm{sin}(u_1+\eta v)\mathrm{sin}(u_2v),`$
$`x_{12}(v)=\mathrm{sin}(u_2v)\mathrm{sin}(u_2v),`$
$`x_{21}(v)=\mathrm{sin}(u_1+\eta v)\mathrm{sin}(u_1\eta v),`$
$`\widehat{x}_{11}(v)=\mathrm{sin}(u_1v)\mathrm{sin}(u_2\eta v),`$
$`\widehat{x}_{22}(v)=\mathrm{sin}(u_1v)\mathrm{sin}(u_2+\eta v),`$
$`\widehat{x}_{12}(v)=\mathrm{sin}(u_2\eta v)\mathrm{sin}(u_2+\eta v),`$
$`\widehat{x}_{21}(v)=\mathrm{sin}(u_1v)\mathrm{sin}(u_1v).`$
We have

$`T(x_{ij}^{}(v);u_1,u_2)B(v+\eta )B(v+2\eta )`$
$`=B(v+\eta )B(v+2\eta )T(\widehat{x}_{ij}^{}(v);u_1,u_2),`$ (10)
where

$`x_{11}^{}(v)=\mathrm{sin}(u_1\eta v)\mathrm{sin}(u_2+3\eta v),`$
$`x_{22}^{}(v)=\mathrm{sin}(u_12\eta v)\mathrm{sin}(u_2v),`$
$`x_{12}^{}(v)=\mathrm{sin}(u_2v)\mathrm{sin}(u_23\eta v),`$
$`x_{21}^{}(v)=\mathrm{sin}(u_1\eta v)\mathrm{sin}(u_12\eta v),`$
$`\widehat{x}_{11}^{}(v)=\mathrm{sin}(u_13\eta v)\mathrm{sin}(u_2\eta v),`$
$`\widehat{x}_{22}^{}(v)=\mathrm{sin}(u_1v)\mathrm{sin}(u_22\eta v),`$
$`\widehat{x}_{12}^{}(v)=\mathrm{sin}(u_2\eta v)\mathrm{sin}(u_22\eta v),`$
$`\widehat{x}_{21}^{}(v)=\mathrm{sin}(u_1v)\mathrm{sin}(u_13\eta v).`$
If we require that the operators $`T(x_{ij}(v);u_1,u_2)`$ and $`T(x_{ij}^{}(v);u_1,u_2)`$ coincide, we can verify that this requirement is satisfied if $`\eta =\pi /3`$ or $`\eta =2\pi /3`$. As discussed above, it follows that the two-parameter family of operators $`T(x_{ij}(v);u_1,u_2)`$ can be restricted on $`W(v)`$.

We conjecture that $`T(x_{ij}(v);u_1,u_2)`$ satisfies restriction conditions (4) and (5) if $`\eta =m\pi /p`$, where $`m`$ and $`p`$ are coprime integers.

It is easy to verify that the relation

$$T_\text{S}(u)=T(x_{ij}(v);u,u)$$

holds, where $`T_\text{S}(u)`$ is Sklyanin transfer matrix in (7).

## 4 Discussion

It was shown in  that the Sklyanin transfer matrix for the Pasqier–Saleur case ($`v\mathrm{}`$) after the restriction satisfies the truncated system of fusion functional equations. This system defines the spectrum $`M(p/p+1)`$. This statement can also be generalized to finite $`v`$. It is remarkable that the spectrum of states surviving after the restriction is independent of $`v`$ . This fact was discovered numerically in .

The explicit construction for additional central elements of the elliptic Yang–Baxter algebra in roots of unity was given in . It would be interesting to generalize the approach in this paper to the elliptic case.

Another important problem is to generalize the Kitanine–Maillet–Terras construction  of the local operators of the $`XXZ`$ model in terms of elements of monodromy matrix to the restricted models. Namely, it would be interesting to build explicit operators that simultaneously respect the restriction and have mutual locality (i.e., commutativity). Constructing such operators would allow obtaining explicit formulas for the correlation functions in the restricted models.

### Acknowledgments

I am indebted to M. Jimbo, N. Kitanine, and E. Sklyanin for the useful discussions and also to W. Everett for the editorial assistance. This work was presented at the conference “On Recent Advances in the Theory of Quantum Integrable Systems.” I indebted to the organizers of the conference and especially to P. Sorba for the opportunity to participate in this nice and very interesting meeting. This work is supported in part by RFBR-01-02-16686, SSRF-20044.2003.2, and INTAS-00-00055.


format: text

text:
# Hysteresis in one-dimensional reaction-diffusion systems

## Abstract

We introduce a simple nonequilibrium model for a driven diffusive system with nonconservative reaction kinetics which exhibits ergodicity breaking and hysteresis in one dimension. These phenomena can be understood through a description of the dominant stochastic many-body dynamics in terms of an equilibrium single-particle problem, viz. the random motion of a shock in an effective potential. This picture also leads to the exact phase diagram of the system and suggests a new generic mechanism for “freezing by heating”.

The closely related questions of phase coexistence, ergodicity breaking and hysteresis in noisy one-dimensional systems with short range interactions and finite local state space (such as in spin systems or vertex models) are intriguing. In thermal equilibrium these phenomena cannot occur as there is no local mechanism that could limit the growth of islands of a minority phase inside a majority phase. Far from equilibrium one has found phase separation and spontaneous symmetry breaking in driven diffusive systems provided that either a bulk conservation law, viz. particle number conservation Evan95 ; Evan98a ; Arnd98b ; Barl98 ; Kafr02 , or vanishing local transition rates Helb99 ; Evan02  constrain the local dynamics. As already noted in Ref. Evan95  the only known exception to this rule, the error-correcting model by Gacs Gacs86 , is rather complicated and still not widely understood, see also Gray01 .

Recently it has been demonstrated that phase coexistence occurs in a one-dimensional driven diffusive system even in the presence of Langmuir kinetics $`A0`$ which break the bulk conservation law Parm03 . This mechanism is inspired by the process of motor proteins moving along actin filaments. Earlier this model was introduced as a toy model reproducing stylized facts in limit order markets Will02 . The formation of a localized shock in this system which separates a domain of low particle density from a domain of high density has been studied subsequently Popk03 ; Evan03 . However, the two different domains do not represent two possible global steady states. The process is ergodic even in the thermodynamic limit and no hysteresis is possible.

It is the purpose of this letter to present a simple nonequilibrium system with local non-conservative dynamics and finite local state space which exhibits ergodicity breaking and hysteresis in the thermodynamic limit, in the usual sense that in finite volume the sojourn time in two metastable steady states increases exponentially with system size. To be specific we investigate the totally asymmetric exclusion process (TASEP) augmented by nonconservative reaction kinetics. The TASEP is a stochastic model of diffusing particles on a one-dimensional lattice with a hopping bias in one direction Schu01 . Each site from 1 to $`L`$ is either empty or occupied by one particle. In the bulk particles (’$`A`$’) hop stochastically from site $`i`$ to $`i+1`$ with unit rate, provided that the target site is empty. The boundaries act as particle reservoirs with densities $`\rho _{}`$ on the left resp. $`\rho _+`$ on the right: On site $`1`$ particles are created with rate $`\rho _{}`$, provided the site is empty, which corresponds to a particle hopping from the left reservoir onto the first site. Particles on site $`L`$ are annihilated with rate $`1\rho _+`$, corresponding to a particle hopping from the last site into the right reservoir.

In our model particles also undergo the following reaction process: On a vacant site enclosed by two particles a particle may be attached with rate $`\omega _a`$, and a particle enclosed by two other particles may be detached with rate $`\omega _d`$. This process can be symbolically written as $`A\mathrm{}AAAA`$ and may be interpreted as activated Langmuir kinetics. Without the TASEP dynamics the stationary density of this process is either $`K=\omega _a/(\omega _a+\omega _d)`$ or zero, with no correlations Rako03b . As in previous work we consider the physically interesting case when $`L\mathrm{}`$ and these rates are proportional to $`1/L`$ Parm03 ; Will02 ; Popk03 ; Evan03 . Hence we define renormalized rates

$$\omega _a=\mathrm{\Omega }_a/L,\omega _d=\mathrm{\Omega }_d/L$$
(1)
where $`\mathrm{\Omega }_a`$ and $`\mathrm{\Omega }_d`$ are kept constant while $`L\mathrm{}`$. For other choices of the attachment/detachment (AD) rates the dynamics is either governed by the TASEP ($`\omega _{a,d}<𝒪(1/L)`$) or by the AD process ($`\omega _{a,d}>𝒪(1/L)`$).

We find a stationary phase diagram of the model with five distinct phases (Fig. 1). The stationary density profile $`\rho _i`$ is not constant as a function of lattice site $`i`$. Yet some of the phases are comparable to those of the usual TASEP with open boundaries Schu93b ; Derr93b : in the high density phases (HD1/2) one finds $`\rho _i>1/2`$ while in the low density phase (LD) $`\rho _i<1/2`$. In HD1 the bulk density profile is dependent on $`\rho _+`$, while it is independent of both boundaries in HD2 as in the maximal current phase of the TASEP. On the other hand two additional phases exist: (i) A coexistence phase which is characterized by a stable shock position, i.e., a jump in the density profile which is localized at a certain position in the bulk of the system. The shock connects a low density domain to its left with a high density domain to its right as known from related models studied previously Parm03 ; Popk03 . Notice that in the usual TASEP there is a coexistence line in the phase diagram with a nonlocalized shock. In a different parameter regime we find a novel phase with an unstable shock position in the bulk. In this phase both the LD and HD states are stable (if $`L\mathrm{}`$) which implies that ergodicity is broken in the thermodynamic limit. Although for finite systems a transition between the two states is possible, the mean life time of each steady state is exponentially large in the system size $`L`$ (see below). We note that this is not a spontaneous symmetry breaking since there is no symmetry relating the two metastable states. This phase has no analog in the TASEP with open boundaries.

Hysteresis in this nonequilibrium setting was observed by measuring the space-averaged density $`\overline{\rho }`$ along the curve of constant $`\rho _+=0.45`$ while changing $`\rho _{}`$ in such a way that the system starting from the LD phase passed through the nonergodic phase and ended up in the HD2 phase. Then the process of changing $`\rho _{}`$ was reversed. A relevant parameter in hysteresis phenomena is the speed of sweeping: in our simulations $`\rho _{}`$ was changed by $`10^4`$ in every $`k`$ MC steps ($`k=500,1500,5000`$). A time average was not taken, $`\overline{\rho }`$ was measured in every $`k`$ steps. On Fig. 2 one can see the resulting hysteresis loops. We found that the hysteresis loop inflates with increasing speed which is reminiscent of hysteresis in usual magnetic systems.

To rationalize these observations we first consider the hydrodynamic limit on the Euler scale, i.e., we take $`L\mathrm{}`$ while the lattice constant is scaled by $`a=1/L`$ and the time by $`t=t_{\text{lattice}}/L`$. Thus the spatial coordinate $`x=i/L`$ becomes continuous. Following the line of arguments of Ref. Popk03  the hydrodynamic equation for the density takes the form

$$\frac{}{t}\rho (x,t)+\frac{}{x}j(\rho )=S(\rho ),$$
(2)
with the exact current $`j(\rho )=\rho (1\rho )`$ of the TASEP and the cubic source term

$$S(\rho )=\mathrm{\Omega }_a\rho ^2(1\rho )\mathrm{\Omega }_d\rho ^3.$$
(3)
resulting from the activated Langmuir kinetics. For the stationary state $`_t\rho (x,t)=0`$ holds and using $`_xj=j/\rho \rho /x`$ we obtain

$$v_c(\rho )\frac{\rho (x)}{x}=S(\rho ),$$
(4)
with the collective velocity $`v_c=j/\rho `$. This nonlinear differential equation can be integrated analytically and yields the flow field

$$x(\rho )=\frac{1}{\mathrm{\Omega }_a\rho }+\frac{\mathrm{\Omega }_a\mathrm{\Omega }_d}{\mathrm{\Omega }_a^2}\mathrm{ln}\left|\frac{1}{K}\frac{1}{\rho }\right|+c$$
(5)
with an integration constant $`c`$.

As the differential equation is of first order and the boundary condition fixes the density at two positions, following a line of the flow field does not represent a solution of the boundary problem in general. In the original lattice model this inconsistency is resolved by the appearance of a shock and/or boundary layers as described in Popk03 . Apart from the discontinuities the stationary density profile follows the flow field of eq. (4).

In order to understand quantitatively the selection of the stationary shock position (which determines the phase diagram) and also to explain the phenomenon of hysteresis from a microscopic viewpoint we describe the dominant dynamical mode of the particle system in terms of the random motion of the shock. To this end we generalize the approach of Kolo98  and introduce space-dependent hopping rates

$`w_{xx+a}={\displaystyle \frac{j_R(x)}{\rho _R(x)\rho _L(x)}}`$
$`w_{x+ax}={\displaystyle \frac{j_L(x)}{\rho _R(x)\rho _L(x)}}.`$ (6)
for jumps of the shock over a lattice constant $`a`$. Here the indices $`L`$ and $`R`$ denote the solutions (lines of the flow field (5)) on the left resp. right of the shock. Similar hopping rates are used in Evan03 . The space-dependent hopping rates furnish us with the picture of a random walker in an effective energy landscape $`E(x)`$ inside a finite box. The energy landscape is generated by the interplay of the particle current with the reaction kinetics. In this way we relate the original nonequilibrium many-particle system to an equilibrium single-particle model. Let $`p(x)`$ be the equilibrium probability of the shock being at position $`x`$. Then due to detailed balance

$$\frac{w_{xx+a}}{w_{x+ax}}=\frac{p(x+a)}{p(x)}=\mathrm{exp}(E(x+a)+E(x)).$$
(7)
which defines the energy landscape.

The potential $`E(x)`$ is monotonically increasing (decreasing) function for the HD (LD) phase (Fig. 3). This implies that although there are fluctuations the shock is always driven to the left (right) boundary. In the coexistence phase there is a global minimum in the bulk resulting in a stable shock position (Fig. 3) at a macroscopic distance from the boundaries. Here the dynamics can be well approximated by a random walker in a harmonic potential which gives a Gaussian distribution for the shock position. Hence the width of the shock distribution is proportional to $`\sqrt{L}`$ Rako03b  which was also found in Parm03 ; Evan03  for the TASEP with Langmuir kinetics.

The nonergodic phase is characterized by a global energy maximum in the bulk (Fig. 3), leading to an unstable bulk fixed point of the shock. The two minima at the left and right boundary correspond to the two stable stationary states. Starting with an initial condition close to one of the minima, the random walker will drift most likely into this local minimum and stay in its vicinity for a time of the order of the mean first passage time $`\overline{\tau }`$ before it traverses to the other minimum. This leads to hysteresis. Using a formula for the mean first passage time derived by Murthy and Kehr MK  one expects that $`\overline{\tau }`$ grows exponentially with the system size $`L`$. Moreover, one expects the transition from one minimum to the other to be a random Poisson process with an average waiting time $`\overline{\tau }`$.

This simple one-particle picture is well borne out by MC simulations. For judiciously chosen parameters it is possible to perform simulations up to times much larger than $`\overline{\tau }`$. Using multispin coding Bark  for the MC algorithm rather good statistics become available for the waiting time $`\tau `$ (the time the system spends in one of the stationary states before switching to the other). For tracing the position of the shock we use the second class particle technique Ferr91 . We measured the position of the second class particle as a function of time: a typical realization is shown in Fig. 4.

As shown in Fig. 5 the numerically determined cumulative distribution function $`\mathrm{\Phi }(t)=P(\tau <t)`$ of the waiting time $`\tau `$ is hardly distinguishable from the expected exponential distribution

$$\mathrm{\Phi }(t)=1\mathrm{exp}(t/\overline{\tau }).$$
(8)
With this picture of a moving shock in mind and using the expression (5) it is also possible to derive the exact phase transition lines defining the phase diagram presented above. Adapting the arguments of Popk03  the analysis of the stability of boundary layers yields a high density phase for $`\rho _{}>1/2`$. Because of a boundary layer the bulk solution of the density profile is independent of $`\rho _+`$ if $`\rho _+<1/2`$. Thus in this region the phase diagram is only ruled by $`\rho _{}`$. For $`\rho _{}1/2`$ the two lines in the phase diagram bounding the coexistence phase and nonergodic phase resp. are determined by the stationary shock position. Crossing the line separating LD/coexistence phase and nonergodic/HD phase from left to right results in a change of the sign of $`_xE(1)`$ from $`+`$. Crossing the other line separating coexistence/HD phase and LD/nonergodic phase from left to right results in a change of the sign of $`_xE(0)`$ from $`+`$. The sign of the slope of the energy profile, i.e., the stability of the shock position can be analysed by considering the average shock velocity

$$v_s=\frac{j_R(x)j_L(x)}{\rho _R(x)\rho _L(x)}.$$
(9)
A shock position at the boundary is stable when it is driven toward the boundary, i.e., $`v_s(0)<0`$ at the left, $`v_s(1)>0`$ at the right boundary. Thus the lines separating the phases are calculated by comparing the values of $`\rho _L(x)`$ and $`\rho _R(x)`$ at the positions $`x=0,1`$.

To conclude we have demonstrated the existence of hysteresis and broken ergodicity (in the thermodynamic limit) in a driven diffusive system without bulk conservation law. We stress that the two different stationary distributions are not ordered states in which the activated Langmuir reaction kinetics would be dynamically suppressed. Surprisingly, adding noise which is on average spatially homogeneous (the nonconservative reaction process) to a conservative spatially homogeneous nonequilibrium system with a nonvanishing particle current leads to a space-dependent effective potential which determines the stationary position of the shock. In the absence of this noise, i.e., in the usual TASEP, the shock performs an unbiased random walk and hence is unlocalized, whereas suitably chosen reaction kinetics may create a variety of effective potentials which localize the shock. An increase in noise strength is usually associated with heating up a system whereas localization reduces the amount of disorder, corresponding to cooling. Thus we have identified a novel generic mechanism for the phenomenon of freezing by heating. The description of the nonequilibrium many-body dynamics in terms of a collective single-particle mode moving under equilibrium conditions yields the exact stationary phase diagram as well as the numerically verified flipping process between the metastable states of the finite system. Details of the flipping dynamics will be presented elsewhere Rako03b .

###### Acknowledgements.

A. Rákos thanks the Deutsche Forschungsgemeinschaft (DFG) for financial support.


text:
# Multiscaling of galactic cosmic ray flux

## I Introduction

Galactic cosmic rays (GCR) originate outside the solar system. They comprise protons (85$`\%`$), alpha particles (14$`\%`$) and heavy nuclei. Interaction with large-scale ordered magnetic field causes the gradient- and curvature-drift motion of GCR in the inner heliosphere, while the interaction with the irregular (stochastic) field component results in the pitch angle scattering of GCR. The scale sizes for the two effects, drifts which depend on variations in the mean field on the order of the heliocentric radial distance, and diffusion which depends on irregularities of sizes comparable to the particles gyro-radii, are quite distinct (see, for instance fs ,alania  and references therein).

Until recently space instruments have lacked the combination of large geometrical factor and good mass resolution required to address the question of the variation on short time scales of individual nuclides other than $`H`$ and $`He`$. In present Letter we report on an investigation of the flux fluctuations of GCR nuclei (Carbon - $`C`$: $`z=6`$) based on the data obtained by instrumentation carried aboard the Advanced Composition Explorer (ACE) spacecraft for 2000 year. The intensities of these low energy particles have been continuously monitored by the Cosmic Ray Isotope Spectrometr (CRIS).

In order to get away from the effects of the Earth’s magnetic field, the ACE spacecraft orbits at the L1 libration point which is a point of Earth-Sun gravitational equilibrium about 1.5 million km from Earth and 148.5 million km from the Sun.

The large collecting power and high resolution of the CRIS instrument allow us to investigate modulation of GCR flux on short time scales (beginning from 1 hour) and in different energy ranges.

Current theories of parallel and perpendicular diffusion are fairly well accepted for high energy particles. However, there are important outstanding issues pertaining to diffusion and transport of charged particles at medium to low energies. The particles at low energies are believed interact directly with the steepened dissipation range that sets in near the ion inertial scale. Dynamical scattering theory suggests that these particles may be sensitive to (1) the dynamical decorrelation mechanism that operates at these scales; (2) the ”geometry” of the fluctuations in the dissipation range; and (3) the spectral distribution of power at the small scales matt1 .

Multiscaling properties of cosmic rays have been already discussed in the light of the data obtained by neutron monitors with cut-off rigity 1.09 GeV (see, for instance, kv ,ber1  and references therein) as well as modulation of the cosmic rays by solar wind turbulence and magnetic field intermittency hor1 -hor2 . However, as far as we know, multiscaling analysis of the GCR heavy nuclei data obtained by satellite instruments is performed here for the first time. This analysis reveals a strong (turbulent-like) intermittency of the GCR dissipation rate and dependence of the intermittency type from the energy ranges (not to be confused with the HEP multiplicity intermittency bp ,ber2 ).

## II Differential flux dissipation rate

Differential flux of the GCR: $`J`$, as measured by CRIS, is proportional to the normal component of the velocity field of the GCR: $`u`$, and to their local concentarion: $`n`$,

$$Jun$$
$`(1)`$
Fluctuations of $`J`$ with time are then determined by corresponding fluctuations of $`u`$ and $`n`$.

Dissipation of passive admixture concentration in fluid turbulence is characterized by a ”gradient” measure my -sa :

$$\chi _r=\frac{_{v_r}(n)^2𝑑v}{v_r}$$
$`(2)`$
where $`v_r`$ is a subvolume with space-scale $`r`$ (for detail justification of this measure see handbook my : p. 381 and further). Scaling law of this measure moments,

$$\chi _r^pr^{\mu _p}$$
$`(3)`$
is an important characteristic of the dissipation rate field just in inertial interval of turbulence (see, for instance, my -sa ). Analogous measure is used to characterize also dissipation rate of turbulent velocity (or kinetic energy) in the inertial interval of scales my -sa .

For turbulent flows the Taylor hypothesis is generally used to interpret the data hor1 ,hor2 ,s . This hypothesis states that the intrinsic time dependence of the wavefields ($`u`$ and $`n`$) can be ignored when the turbulence is convected past the probes at nearly constant speed. With this hypothesis, the temporal dynamics should reflect the spatial one, i.e. the fluctuating velocity (concentration) field measured by a given probe as a function of time, $`u(t)`$ is the same as the velocity $`u(x/u)`$ where $`u`$ is the mean velocity and $`x`$ is the distance to a position ”upstream” where the velocity is measured at $`t=0`$.

With the Taylor hypothesis $`dJ/dx`$ is replaced by $`dJ/udt`$ and one can define GCR flux dissipation rate as:

$$\chi _\tau \frac{_0^\tau (\frac{dJ}{dt})^2𝑑t}{\tau }$$
$`(4)`$
where $`\tau r/u`$ and corresponding scaling of the dissipation rate moments as sa :

$$\chi _\tau ^p\tau ^{\mu _p}$$
$`(5)`$
Substituting (1) into (4) one can estimate the GCR flux dissipation rate through characteristics of $`u`$ and $`n`$. We can consider two asymptotic regimes. For sufficiently small GCR energies the flux dissipation rate has been dominated by the concentration dissipation:

$$\chi _\tau ^p[\frac{1}{\tau }_0^\tau (\frac{dn}{dt})^2𝑑t]^p\tau ^{\mu _p}$$
$`(6)`$
while for sufficiently large GCR energies the flux dissipation rate has been dominated by the energy (velocity) dissipation:

$$\chi _\tau ^p[\frac{1}{\tau }_0^\tau (\frac{du}{dt})^2𝑑t]^p\tau ^{\mu _p}$$
$`(7)`$
## III The data

As it was mentioned in Introduction we will use the data collected by ACE/CRIS instrument during 2000 year (the year of maximum solar activity). We will consider carbon nuclei $`C`$ ($`z=6`$), due to carbon is the lightest abundant nuclei (after $`H`$ and $`He`$) in the ACE/CRIS collection. As we shall see this fact allow us consider the both asymptotes mentioned in previous Section.

Figure 1a shows scaling of the GCR flux dissipation rate moments $`\chi _\tau ^p`$ (4)-(5) for the C-nuclei with energies from energy range: 56.3-73.4 Mev/nucleon (the lowest energy range in the ACE/CRIS collection for carbon).

Figure 1b shows the scaling exponents $`\mu _p`$ (circles) extracted from figure 1a (as slopes of the straight lines, cf (5)). The solid curve in figure 1b corresponds to the intermittency exponents $`\mu _p`$ obtained for the inertial-convective region of a passive admixture concentration in a laboratory turbulent air flow sa . To support the striking correspondence between the two data sets multiscaling let us calculate also the Extended Self-Similarity (ESS ben ) exponents, $`\beta _p`$, extracted from equation:

$$\chi _\tau ^p\chi _\tau ^3^{\beta _p}$$
$`(8)`$
The ESS of type (8) usually has clearer scaling form than ordinary scaling (5) and covers a wider range of scales (see ben  for a review of ESS and for examples). Figure 2a shows the ESS of the GCR flux dissipation rate moments $`\chi _\tau ^p`$ (8) using log-log scales. The straight lines (the best fit) are drawn to indicate the ESS (8). Figure 2b shows the ESS exponents $`\beta _p`$ (circles) extracted from figure 2a (as slopes of the straight lines, cf (8)). The solid curve in figure 2b corresponds to the intermittency exponents $`\beta _p`$ obtained for the inertial-convective region of passive admixture concentration in different laboratory turbulent flows sa .

Now let us turn to the data corresponding to energy range with the highest for the carbon nuclei energies observed by ACE/CRIS. Figure 3a shows scaling of the GCR flux dissipation rate moments $`\chi _\tau ^p`$ (4)-(5) for the C-nuclei with energies from energy range: 183.1-198.7 MeV/nucleon.

Figure 3b shows the scaling exponents $`\mu _p`$ (circles) extracted from figure 3a (as slopes of the straight lines, cf (5)). The solid curve in figure 3b corresponds to the intermittency exponents $`\mu _p`$ calculated using the She-Leveque model sl , which is in very good agreement with the data for velocity (kinetic energy dissipation) field intermittency obtained in inertial interval for isotropic fluid turbulence. Figures 4a and 4b shows corresponding Extended Self-Similarity (ESS) properties

$$\chi _\tau ^p\chi _\tau ^4^{\beta _p}$$
$`(9)`$
observed for the data. Again, as in figure 3b, correspondence to the fluid turbulence (energy dissipation rate - the solid curve in figure 4b) is very good.

The interval of time scales under consideration is 1-30 hours. It seems to be useful to compare the observed properties of GCR flux dissipation rate with relevant properties of the local interplanetary magnetic field. Figure 5 shows energy spectrum of 3D magnetic field measured by ACE/MAG magnetometer in the same time intervals as the ACE/CRIS data. The slope -5/3 indicates Kolmogorov-like scaling (cf gold ). Let us recall that the turbulent fluid dissipation rates used for comparison in figures 1-4 were obtained just for inertial (inertial-convection) interval of scales where the Kolmogorov-like scaling should to be expected sa ,sl .

## IV Discussion

The results of previous Section seems to be in agreement with the qualitative estimates made in Section 2: for relatively small energies the flux dissipation rate behaves as one dominated by the GCR particles concentration dissipation rate, whereas for relatively large energies the GCR flux dissipation rate seems to be dominated by velocity (kinetic energy) dissipation.

Since CRIS instrument measures the differential flux ($`J`$) only, but not the velocity of GCR ($`u`$) and their concentration ($`n`$) separately (cf (1)) we have no idea about magnitudes of the fluctuations of $`u`$ and $`n`$ themselves. Therefore existence of the two asymptotes (6) and (7) is still a pure phenomenology (for the intermediate energies the scaling is deformed by competition between the two different mechanisms of the flux dissipation). The same reason (unknown $`u`$ and $`n`$ for the GCR) makes it impossible to perform any theoretical estimates for the scales of the GCR flux dissipation rate multiscaling.

Although correspondence between the motions in MHD and hydro cases has been reported earlier (see, for instance, very recent paper cho2  and references therein) it is much more difficult to understand such detail correspondence between cosmic rays and fluid turbulence. Indeed, as it follows from previous two sections the observed GCR particles dissipation rate exhibits intermittent properties of a passive admixture convected by turbulent motion of a classic (isotropic, incompressible) non-magnetic fluid. While the solar wind (including the interplanetary magnetic field - IMF) is certainly very different from such a fluid. Moreover, interaction of the electrically charged GCR particles with the solar wind and the magnetic field is expected to have strong resonance features (see, for instance, chan -clv ,cho2  and references therein). Therefore, even the fact that many characteristics of solar wind and IMF turbulence are known to be similar to those of turbulence of classic non-magnetic fluids (see hor1 -hor2  and references therein) the presumably resonance-like interactions of GCR particles with solar wind and IMF seems to be non-consistent with the observed here picture, that request reconsideration of our approach to stochastic short-term convection of the low-energy GCR particles (heavy nuclei) in the heliosphere and their short-term interactions with solar wind and IMF. This reconsideration represents a serious challenge for a wide scientific community: both for astrophysicists and for plasma experts.

The author is grateful to ACE/CRIS and to ACE/MAG instrument teams as well as to the ACE Science Center for providing the data. Numerous discussions on the subject with K.R. Sreenivasan as well as the Referees’ comments and suggestions were very useful for this investigation.

Figure Captions

Figure 1a. The GCR flux dissipation rate moments $`\chi _\tau ^p`$ against $`\tau `$ for the C-nuclei with energies from energy range: 56.3-73.4 Mev/nucleon. The time interval $`\tau `$ is measured in hours, whereas the GCR differential flux $`J`$ is measured in $`10^{7.5}particles/m^2ssrMev`$. Log-log scales are chosen in the figure for comparison with scaling equation (5). The straight lines (the best fit) are drawn to indicate the scaling. The upper data sets correspond to larger $`p=2,3,4,5`$

Figure 1b. The scaling exponents $`\mu _p`$ (circles) extracted from figure 1a. The solid curve corresponds to the intermittency exponents $`\mu _p`$ obtained for the inertial-convective region of a passive admixture concentration dissipation rate in a laboratory turbulent air flow sa .

Figure 2a. The ESS of the GCR flux dissipation rate moments $`\chi _\tau ^p`$ against $`\chi _\tau ^3`$ in log-log scales (8) for energy range: 56.3-73.4 Mev/nucleon. The straight lines (the best fit) are drawn to indicate the ESS (8). The upper data sets correspond to larger $`p=2,3,4,5`$.

Figure 2b. The ESS exponents $`\beta _p`$ (circles) extracted from figure 2a. The solid curve corresponds to the intermittency exponents $`\beta _p`$ obtained for the inertial-convective region of passive admixture concentration dissipation rate in different fluid turbulent flows sa .

Figure 3a. The GCR flux dissipation rate moments $`\chi _\tau ^p`$ (4)-(5) against $`\tau `$ for the C-nuclei with energies from energy range: 183.1-198.7 MeV/nucleon.

Figure 3b. The scaling exponents $`\mu _p`$ (circles) extracted from figure 3a. The solid curve in figure 3b corresponds to the intermittency exponents $`\mu _p`$ calculated using the She-Leveque model for isotropic fluid turbulence sl .

Figure 4a. The ESS of the GCR flux dissipation rate moments $`\chi _\tau ^p`$ against $`\chi _\tau ^3`$ in log-log scales (9) for energy range: 183.1-198.7 MeV/nucleon. The straight lines (the best fit) are drawn to indicate the ESS (9).

Figure 4b. The ESS exponents $`\beta _p`$ (circles) extracted from figure 4a. The solid curve corresponds to the intermittency exponents $`\beta _p`$ calculated using the She-Leveque model.

Figure 5. Energy spectrum of 3D magnetic field measured by ACE/MAG magnetometer in the discussed range of the time-scales.


text:
# Time machine (1988–2001)

> “…don’t drive yourself crazy trying to resolve the paradoxes of time travel.”
>
> L. S. de Camp

The time machine (in its science fiction incarnation) was born back in 19-th century , i. e. well before the concept of spacetime appeared. However, it was relativity that provided an adequate language for discussing and studying time machines. In particular, according to special relativity, our world is the Minkowski space and the life history of a pointlike particle is a curve in this space, see Fig. 1.

In these terms one can formulate, in quite a meaningful way, the question about the existence of time machines, or — which is the same — about the possibility of an observer meeting his older self. Namely, one can ask, whether a world line of a particle may have a self-intersection. The answer is straightforward. To make a loop a curve must somewhere leave the null cone as shown in Fig. 1. A particle with such a world line would exceed the speed of light and, insofar as no tachyons have been found yet, the time machine is impossible in special relativity.

In general relativity the situation is much less trivial.

According to this theory, our spacetime must be a smooth Lorentzian manifold, i. e. in sufficiently small regions it must be ‘approximately Minkowskian’, but at large scale (as long as the Einstein equations with a reasonable matter source hold), it may have any geometry and topology. There may be holes, handles, see Fig. 2a; almost whatever one wants. And the null cones need not any more look all in the same direction. In particular, spacetimes are conceivable with closed or self-intersecting timelike curves.

A simple example is the Minkowski space rolled into a cylinder, see Fig. 2b. Locally everything is fine in this spacetime, but due to its non-trivial global structure, an observer can meet his younger self. Such spacetimes , however, do not deserve the name ‘time machine’, because the closed timelike curves *always* exist in them. They are not *created* at some moment. In this sence, a better example is the Deutsch–Politzer (DP) space , which we shall describe as a result of the following surgery. Make two cuts in the Minkowski space

and glue the upper bank of lower cut to the lower bank of the upper cut and vice verse, as shown in Fig. 3a, so that a cylinder appears attached to the plane. It is convenient to draw the resulting spacetime still as the Minkowski plane, see Fig. 3b, and just to keep in mind that any curve reaching the lower segment must be continued from the upper and, vice verse.

The DP space contains timelike loops, which all lie to the future of a good (causal) region, so it fits the name ‘time machine’. This model, of course, in unsuitable for tackling the problem of *building* a time machine (it gives no clue into what could *force* the spacetime to evolve into the time machine instead of just remaining Minkowskian, see below), but is quite adequate in studying what can happen, if for some reason or other a time machine did appear.

One can often read that the existence of a time machine would lead to awful paradoxes contradicting common sense and the notion of free will, or even proving the impossibility of the time machine. Now we can verify this assertion. Let us begin with the notorious grandfather paradox, which is a story of a researcher $``$, who asks his grandson $`𝒢`$ to trip (by time machine) to the past and to kill him ($``$) in infancy. The grandson tries to grant his wish

(see Fig. 4a), but the problem is: if $`𝒢`$ succeeds than the baby does not grow up and, correspondingly, does not have any grandchildren that could kill him. Which proves that the grandson *fails* to kill the baby. The circumstances are not detailed, so there is much room for explanations *just why* he fails.

Say, the gun can jam, or the shooter can miss. However, if we arm grandsons again and again (hundreds of times) and *always* only bad shooters get reliable guns, this, of course, would look miraculous. Why cannot we give a good gun to a marksman? Is not this a restriction on our freedom of will?

In fact, the solution of the grandfather paradox is quite simple. The point is that only those persons are *available* for our experiment whose grandfathers by whatever reasons were *not* killed. So, what prevents us from arming a sufficiently good shooter with a sufficiently good gun is not problems with our freedom of will111The grandson, being born with already known future life, is such an unusual object that we need not worry about *his* freedom of will., but just the shortage of such shooters. And the better our guns are the fewer marksmans are at our disposal.

A more detailed analysis shows  that to avoid such simple solutions one must consider experiments performed *before* the time machine appeared . Indeed, suppose one kindly requests a person to wait until a time machine is built (at some time in the *future*), to enter it, and to kill his (the traveler’s) younger self (see Fig. 4b). The situation in the time machine is exactly as paradoxical as in the grandfather case — if the test person is dead when he meets his younger self, he cannot kill the latter, so why he is dead? And if he is alive and shooting, then how come the victim survived? This time, however, the initial conditions do not depend on what occurs in the time machine and the paradox has no simple resolutions. Moreover, being formulated in terms of pointlike particles (which allows one to check *all* possibilities) it has no resolutions *at all* .

How to fight the time travel paradox? The most obvious way out would be to forbid the time machines by fiat, to *postulate* causality. But this is not so easy. In 1988 Morris, Thorne, and Yurtsever (MTY) published a paper, where they considered a spacetime $`M`$ shown in Fig. 5.

Geometrically $`M`$ can be obtained by removing two close cylinders $`C_1`$ and $`C_2`$ from the Minkowski space and gluing together the boundaries $`_{1,2}\mathrm{Bd}C_{1,2}`$ of the resulting holes (the junction is then smoothed out by curving appropriately a close vicinity of $`_1=_2`$). The cylinders are taken to be parallel to the $`t`$-axis except that $`C_2`$ has a bend at, say, $`1<t<1`$. The gluing must respect the following rule

$$p=p^{}\tau (p)=\tau (p^{}),$$
($``$)
where $`\tau (x)`$ is defined to be the length of the longest (recall that the metric is Lorentzian) timelike curve lying in $`_{1,2}`$ and connecting $`x`$ with the surface $`t=2`$. Due to the bend, $`\tau `$ grows with time slower on $`_2`$ than on $`_1`$ (the ‘twin paradox’). Correspondingly, some of the identified points were causally related in the ‘initial’ (Minkowski) space. So, there are timelike loops in $`M`$. Physically, $`M`$ presents a spacetime in which a wormhole (see Fig. 5b) evolves so that the distance (in the ambient flat space) between its mouths changes at $`1<t<1`$ without significant changes in the form or length of the throat (it is this last condition that necessitates the rule ($``$Time machine (1988–2001))).

The MTY time machine cannot be banished as easily as the DP one, because we have (in assumption that wormholes exist, they are large and stable, etc.) a specific prescription of how to build it. Precisely this was the point of : assuming the time machines *are* prohibited, exactly how this prohibition is enforced? Suppose, one finds a wormhole, pushes one of its mouths, then pulles it back. What *exactly* will go wrong? What will prevent a closed causal curve from appearance? It was such statement of the problem that gave birth to the time machine as an element of the physical, rather than philosophical, or science fiction realm.

A mechanism that could prevent the appearance of a closed timelike curve was sought for more than a decade, but has never been found . In this sense causality still remains unprotected . But does it mean that an advanced civilization can build a time machine and face the ensuing paradoxes? No, it does not. The point is that there is a great difference between *observing* a time machine appearing by itself and *creating* one. And while the former seemingly is possible, the latter is not. The difference is well exemplified by the DP space. We can prepare a flat region in our spacetime and, for all we know, a time machine well may appear there. But just as well it may not, and we shall have a usual Minkowski space. The choice is up to the spacetime and we can do nothing to affect it. Surprisingly, exactly the same is true for *any* time machine, be it based on a wormhole , a cosmic string , or whatever else . A theorem  proved in 2001 says: *Any spacetime obeying local laws and having no causal loops in its past, has a causal maximal extension.* Which means (since the Einstein equations are local) that within general relativity a time machine cannot be built.


text:
# Chapter 1 Search for TeV annihilation radiation from supersymmetric dark matter in nearby galaxies

## Chapter 1 Search for TeV annihilation radiation from supersymmetric dark matter in nearby galaxies

### Abstract

During the 2002-2003 observing season the Whipple 10m imaging atmospheric Čherenkov telescope was used to search for dark matter annihilation radiation in four nearby galaxies: M32, M33, Draco, and Ursa Minor. Scientific motivations for this choice of targets are discussed as well as accumulated exposure. The analysis results are to be reported in the conference presentation.

### 1. Introduction

The lightest supersymmetric particle, the neutralino, with a mass in the range $`50`$ GeV - $`5`$ TeV, is a plausible candidate for non-baryonic cold dark matter (CDM) , which can be detected indirectly via its annihilation products . The annihilation rate is proportional to the square of the neutralino density integrated along the line of sight, suggesting strong enhancement in the direction of dark matter clumps. The increasing dark matter density profile toward the center of Milky Way (MW) galaxy and its proximity made the Galactic Center a natural choice for such searches . The CDM annihilation flux from even nearby extragalactic objects, such as dwarf galaxies located only ten times father, would have to be suppressed by a factor of hundred. However, due to the squared neutralino density dependence, the distance penalty can be overcome by the equal increase of DM density in the cores of MW satellite and local group galaxies. Although little is known about the distribution of DM in the central $`1`$ \- $`10`$ pc of Galactic Nuclei (GN), which is the most relevant region for generating an annihilation signal, the current paradigm suggests that it is a strong function of the merger history of a galaxy  and its long-term evolution that could be affected by the presence and growth of central black holes (BHs), e.g. . Because of this, D. Merritt argues that MW was perhaps unfortunate, and a DM spike in its center is unlikely . On the other hand, if stellar, and stellar velocity distributions can be used as a guide to the distribution of dark matter at the very centers of galaxies, we suggest the following arguments to justify the choice of several extragalactic sources as plausible candidates for observation of DM annihilation.

### 2. M32

M32 is the closest compact elliptical galaxy believed to be formed in a rare cosmological event from a low luminosity spiral galaxy when it plunged into the central region of M31 and most of the outer stellar component was tidally stripped . Stellar kinematic data as well as gas-dynamical studies strongly support the presence of a single supermassive compact object, $`3.6\times 10^6M_o`$ in the center of M32 . The core of M32 has a relatively homogeneous stellar population that can be modeled as being coeval and of intermediate age ($`4`$ Gyr) . Lauer et al. estimate M32 core relaxation time scale as $`23`$ Gyr . These data suggest that the nucleus of M32 was unlikely to have undergone a recent merging event when a massive black hole binary could have depleted the central stellar density by evacuating stars and destroying potential dark matter cusps in the galaxy core . The data also imply that enough time has passed for collisional two-body relaxation of stellar population around a black hole to form a cusp in the Stellar Density Profile (SDP), $`1/r^\alpha `$, with $`\alpha =3/27/4`$ correspondent to Bahcall-Wolf solution . Optical and infrared data indicate SDP compatible with $`\alpha `$ in the range $`1.4`$ to $`1.9`$ and still rising at the resolution limit of $`0.07`$ pc. The exact spectral index and the possibility of its break at a distance $`1`$ pc from the GN center is being debated . The stellar density at the very center of M32’s nucleus likely exceeds $`10^7M_o`$ per pc3, which is the highest known to us in nearby systems .

After initial violent relaxation during large scale structure formation non-interacting dark matter has no means of self-evolution on subgalactic scales because it cannot cool down without gravitational coupling to baryonic matter. This is true at least for spherically symmetrical systems. Thus it is likely that during formation of the galactic core the DM Density Profile (DMDP) should follow SDP with a characteristic relaxation time comparable to stellar relaxation that, in turn, is decreased by the heat transfer associated with the DM component in the nucleus. In addition, the slow growth of the central BH should cause adiabatic compression of DM with the adiabatic invariant $`rM_{BH}=`$const . If almost all BH mass is acquired by slow growth, the amount of compression could be substantial, resulting in the appearance of a spike in the DMDP with $`\alpha >2.25`$ . If, however, galactic mergers have been dominantly responsible for BH growth this may not only destroy potentially created DM spikes, it will also drastically suppress the amount of adiabatic DM compression around the BH. Probably both the effect of adiabatic BH growth and the evolution of DMDP driven by baryonic stellar and gas components in the cores of the galaxies take place at the same time and cannot be completely disentangled . For the case of M32 both scenarios of DM evolution seem to favor a cusp in DMDP with $`\alpha 1.5`$.

### 3. M33

The neutralino annihilation flux from a “cuspy” DMDP ($`\alpha >1.5`$) is formally divergent at small scales. The physically justified minimal radius is set by the equality between the annihilation rate and the rate of supply of neutralinos into annihilation region. The latter is related to the BH growth rate and/or relaxation rate of baryonic & DM components in GN. The general effect of the central BH is to increase the stellar velocity dispersion within the radius of its gravitational influence , and consequently increase GN relaxation time, limiting neutralino annihilation flux dynamically. From this point of view a galaxy with a small BH and very rapid GN evolution scale, such as M33, may be preferable for observations.

M33 is a normal low-luminosity dark matter dominated bulgeless spiral galaxy with dark halo $`5.1\times 10^{11}M_o`$ . The mass of the BH in its center is less than $`1.5\times 10^3M_o`$ , and the stellar population in its nucleus can be modeled by two bursts of star formation $`2`$ and $`0.5`$ Gyr ago . The nucleus of M33 hosts the most luminous steady X-ray source in the Local Group that is also associated with a radio source and reminiscent of the galactic microquasar GRS 1915+105 . Small ($``$%$`10`$) but significant X-ray spectrum and time variability have also been reported . M33 is remarkable for its very small galaxy nucleus relaxation time $`3`$ Myr due to very dense, $`>2\times 10^6M_o`$ per pc3, stellar core ($`1`$ pc) and extremely low velocity dispersion .

### 4. Draco and Ursa Minor dwarf galaxies

Dark matter rich MW satellite Dwarf Galaxies (DGs) are favored for observations due to their proximity, very low baryonic content, and/or possibility of self-interacting DM. Within non-interacting DM scenarios the cores of these galaxies are “frozen” due to particularly low stellar and gas densities and consequently slow evolution rates. Their typical core relaxation times, $`5\times 10^2`$ Gyr, exceed Hubble time. Thus, the DGs could not have evolutionarily built any sub-parsec cusp or spike-like structure in their DMDP for efficient annihilation of neutralinos, and it is also difficult to imagine any DM evolutionary scenario on the scale of $`1`$ pc, a typical separation between stars in the nuclei of DGs. Tyler argues that these perhaps might be the exact conditions which would preserve intact any sub-parsec DM structures from destruction by galaxy merger events or by particularly high abundance of stars in GN . The initial perturbations in DMDP with the scales smaller than the Jeans instability scale could have formed during violent relaxation and survived without trapping a substantial amount of baryons . These invisible and very slowly evolving structures could have been gravitationally trapped in the MW halo or in the centers of dwarf galaxies. If so, the quantitative prediction of the annihilation flux from DM cusps will strongly depend on the cusp evolution mechanism because it will determine the maximal annihilation rate. This scenario becomes particularly interesting if DM is self-interacting indeed . Among many other puzzling observations of dwarf galaxies there are indications of “clumpy” distribution of matter. For example, Ursa Minor dwarf has distinct stellar lumps within $`10^{}`$ circle . The study of stellar proper motion suggests that the lump crossing time is $`5`$ Myr and six individual lumps shouldn’t exist any longer than this . Is it a coincidence to observe six of them at the same time?

### 5. Data

Table 1 shows the range of zenith angles at which all proposed targets were observed and the total exposure so far accumulated by the VERITAS collaboration for each source.

### Acknowledgments

The authors thank E. Roache and J. Melnick for their technical assistance. We gratefully acknowledge support from the VERITAS collaboration and University of Utah. This recearch is supported by National Science Foundation under NSF Grant #0079704 and by grants from the U.S. Department of Energy, Enterprise Ireland and PPARC in the UK.

### 6. References

1. Bahcall J.N. and Wolf R.A. 1976, ApJ 209, 214

2. Bekki K. et al. 2001, ApJ 557, L39

3. Bergstrom L., Ullio P. 1997, Nucl. Phys. B504, 27

4. Bergstrom L., Ullio P., & Buckley J. 1998, Astroparticle Phys. 9, 137

5. Cesarini A. et al. 2003, astro-ph/0305075

6. Corbelli E. & Salucci P. 2000, MNRAS 311, 441

7. Corbin M.R. et al. 2001, AJ 121, 2549

8. Del Burgo C. et al. 2001, MNRAS 321, 227

9. Dubus G. & Rutledge R.E. 2002, MNRAS 336, 901

10. Eskridge P.B. & Schweitzer A.E. 2001, AJ 122, 3106

11. Gebhardt K. et al. 2001, AJ 122, 2469

12. Goldberg H. 1983, Phys. Rev. Lett. 50, 1419

13. Gondolo P. & Silk J. 1999, Phys. Rev. Lett. 83, 1719

14. Hennawi J.F. & Ostriker J.P. 2002, ApJ 572, 41

15. Joseph C.L. et al. 2001, ApJ 550, 668

16. Jungman G., Kamionkowski M., & Griest K. 1996, Phys. Rev. 267, 195

17. Kleyna J.T. et al. 1998, AJ 115, 2359

18. Klypin A. et al. 1999, ApJ 522, 82

19. Kosak K. et al. 2003, This proceedings

20. Landau L.D. & Lifshitz E.M. 1960, in “Mechanics” (Oxford, Pergamon Press)

21. La Parola V. et al. 2003, ApJ 583, 758

22. Lauer T.R. et al. 1998, AJ 116, 2263

23. Merritt D. & Quinlan G.D. 1998, ApJ 498, 625

24. Merritt D., Ferrarese L, Joseph C.L. 2001, Science 293, 1116

25. Merritt D. et al. 2002, Phys. Rev. Lett. 88, 1301

26. Merritt D. 2003, astro-ph/0301365

27. Milosavljevic M. et al. 2002, MNRAS 331, L51

28. Moore B. et al. 1999, MNRAS 310, 1147

29. Stephens A.W. & Frogel J.A. 2002, AJ 124, 2023

30. Tremaine S. et al. 1994, AJ 107, 634

31. Tyler G. 2002, Phys. Rev. D66, 3509

32. Ullio P., Bergstrom L. 1998, Phys. Rev. D57, 1962


text:
# Locally preferred structure in simple atomic liquids

## I Introduction

A liquid is said to be supercooled when it is possible to cool it below its melting temperature $`T_m`$ without crystallizing. The supercooled liquid phase is metastable with respect to the underlying crystal, and it is characterized by a dramatic increase of the viscosity and the relaxation times upon lowering the temperature, an increase that eventually leads to glass formation. These dramatic changes in the dynamic properties are not accompanied by strong signatures in the structural quantities, such as the static structure factor. Yet, it has been suggested that both supercooling and glass formation were deeply connected to the structure of the liquid, more precisely to a competition between extension of a local liquid order, different than that of the crystal, and global constraints associated with tiling of the entire space sadoc99 ; nelson89 ; kivelson95 . This competition has been termed geometric (or topological) frustration sadoc99 ; nelson89 .

Some fifty years ago, Frank put forward the following argument to explain supercooling of liquids frank52 . If one considers atomic liquids in which atoms interact via spherically symmetric potentials like the Lennard-Jones potential, the local arrangement of the atoms that is preferred is not the clusters associated with crystalline order (face-centered cubic and hexagonal close-packed lattices), but a polytetrahedral packing, the icosahedron formed by $`13`$ atoms. The energy of such an icosahedral arrangement interacting via the Lennard-Jones potential is indeed $`8.4\%`$ lower than the close-packed crystalline clusters. Local icosahedral order should then be prevalent in liquids but, because of the $`5`$-fold rotational symmetry of the icosahedron, it cannot tile the entire space and form a crystal: this is a manifestation of geometrical frustration. Crystallization then requires a rearrangement of the local structure of the liquid, which leads to a strong first-order freezing transition and allows supercooling of the liquid. Since then, there has been a large body of experimental and simulation work confirming the prevalence of local icosahedral, or more generally polytetrahedral, order in atomic liquids and metallic glasses nelson89 ; filipponi99 ; reichert00 ; shenk02 ; steinhardt83 ; jonsson88 ; kondo90 ; yonezawa91 ; shumway95 ; jund97 ; dzugutov02 ; doye03 . The tendency to form icosahedral order has been shown to increase as the temperature is lowered nelson89 ; jonsson88 ; kondo90 ; yonezawa91 ; shumway95 ; jund97 ; dzugutov02 ; doye03 .

Frank’s argument has, however, one shortcoming doye01 : the $`13`$-atom cluster considered is isolated, so that most of the energetics is related to the surface, a situation that of course does not occur in bulk liquids. How can one remedy this problem? An a priori easy way would be to study directly the local arrangement of the atoms in a bulk liquid. In such a case, however, it is found that the proportion of icosahedra among all the $`13`$-atom groups is very small, typically a few percent kondo90 . This is indeed to be expected since geometrical frustration is present, which opposes the growth of icosahedral order and distorts the local polytetrahedral arrangements nelson89 ; sachdev92 . How to disentangle then the determination of the locally preferred structure from frustration effects? The method we propose in this work is to consider the influence of the bulk liquid on a given $`13`$-atom cluster at some mean-field level, so that surface effects can be reduced and made more realistic for describing a condensed phase, whereas geometrical frustration is strongly inhibited. The main advantage of this method is that it can be extended to study the locally preferred structure of molecular liquids, for which a priori topological arguments do not easily provide the symmetries of all possible local arrangements of the molecules, nor the nature of geometrical frustration.

More specifically, in this article we consider the ground-state of a cluster of $`13`$ atoms interacting via a Lennard-Jones potential (smaller and larger clusters are also considered); the atoms are placed in a cavity and are subject to an external field that mimics the interaction with the rest of the liquid. The structure of the outside liquid only enters the calculation via the bulk pair distribution function (known from previous simulation studies). By means of an optimization algorithm wales97 ; doye98  we find that the ground-state of the cluster, i.e., the global minimum of the (effective) energy surface formed by the intra-cluster interactions and the external field, is of icosahedral symmetry, therefore generalizing Frank’s result. For sake of comparison we consider in addition other boundary conditions for the cluster. These conditions describe different types of environments: free boundary conditions (isolated cluster), periodic boundary conditions (periodic tiling of space), and icosahedral-like boundary conditions (hypothetical non-frustrated system).

## II Method and choice of the boundary conditions

We consider a system of atoms interacting via a pair-wise additive spherically symmetric potential $`v(r)`$, where $`r`$ is the center-to-center distance. A number $`N`$ of atoms are placed in a spherical cavity $`𝒞`$ of radius $`R_C`$, that we envisage as surrounded by bulk liquid made of the same atoms, and characterized by the temperature $`T`$ and the density $`\rho `$. As explained above, we do not want to fully account for the liquid structure because geometrical frustration would obscure the nature of the local order. We rather resort to a mean-field type of description in which the liquid outside the cavity is considered as a continuum, characterized by a (known) pair distribution function $`g(r;T,\rho )`$ that is not affected by the fact that a cavity has been carved out. The potential energy acting on a given atom at position $`𝐫`$ inside the cavity due to the outside liquid is thus described as

$$W(𝐫;R_C)=\frac{\rho }{2}_{𝐫^{}𝒞}d^3𝐫^{}g(|𝐫𝐫^{}|)v(|𝐫𝐫^{}|),$$
(1)
where the integral is over all positions outside the cavity. Taking the center of the cavity as origin and transforming to spherical coordinates, one finds after standard manipulations

$$W(r;R_C)=\pi \rho \left\{_{R_Cr}^{R_C+r}𝑑xx^2g(x)v(x)\left[1u(x;r,R_C)\right]+2_{R_C+r}^{\mathrm{}}𝑑xx^2g(x)v(x)\right\},$$
(2)
where

$$u(x;r,R_C)=\left(\frac{R_C^2r^2x^2}{2rx}\right).$$
(3)
Note that $`u(x=R_Cr)=u(x=R_C+r)=1`$. The detail of the calculation is given in the Appendix A.

In view of the implementation of the optimization algorithm wales97 ; doye98 , we need an explicit expression for the first derivatives of the external potential, both with respect to $`R_C`$ and $`r`$; we obtain

$$\frac{}{R_C}W(r;R_C)|_r=\pi \rho \left(\frac{R_C}{r}\right)_{R_Cr}^{R_C+r}𝑑xxg(x)v(x),$$
(4)
and

$$\frac{}{r}W(r;R_C)|_{R_C}=\pi \rho _{R_Cr}^{R_C+r}𝑑xx^2g(x)v(x)\left(\frac{R_C^2+r^2x^2}{2xr^2}\right).$$
(5)
The second derivatives, needed for the calculation of the Hessian matrix and the study of the transition states wales97 ; doye98 , can be obtained in a similar way.

The total potential energy for the $`N`$ atoms of the embedded cluster is the sum of the atom-atom interaction potentials inside the cavity and of the external potential,

$$𝒰(\left\{𝐫_j\right\}_{1,\mathrm{},N};R_C)=\underset{i>j=1}{\overset{N}{}}v(|𝐫_i𝐫_j|)+\underset{j=1}{\overset{N}{}}W(r_j;R_C),$$
(6)
where $`W(r;R_C)`$ is given by Eqs. (2) and (3). Finding the ground-state configuration for the $`N`$-atom cluster embedded in a liquid-like environment (at a given $`T`$ and $`\rho `$) amounts to determining the global minimum of $`𝒰`$ with respect to variations of the positions of the $`N`$ atoms. (Note that $`T`$ and $`\rho `$ only enter trough the external potential $`W(r;R_C)`$, both explicitly (see Eq. (2)) and through the state-dependence of the pair distribution function $`g(r)`$.) It is worth stressing that, contrary to the case of an isolated cluster for which the radius of the cavity $`R_C`$ is merely fixed to avoid evaporation of the atoms, $`R_C`$ becomes a relevant variable in a liquid-like environment: to preserve a realistic description of the liquid, $`R_C`$ should adjust to global contractions or expansions of the $`N`$-atom cluster, which should then be taken into account in the minimization procedure. This point will be further discussed below.

To obtain the lowest energy minimum of $`𝒰`$, we have used a slightly modified version of the basin-hopping algorithm introduced by Wales and co-workers wales97 ; doye98 . The algorithm consists of a constant-temperature Monte-Carlo simulation performed with an acceptance criterion based not upon the energy of the proposed new configuration, but upon the energy of the closest minimum of the potential energy surface, obtained by a local minimization starting from that configuration. This algorithm turns out to be a very efficient method for exploring directly the minima of the potential energy surface, and it allows one to locate the ground-state with relatively little effort.

Finally we have also considered other boundary conditions for the $`N`$-atom cluster:

i) Free boundary conditions (isolated cluster).– This is the standard case studied in the literature. It simply corresponds to the above situation in which the external potential is set to zero:

$$W(r;R_C)=0.$$
(7)
ii) Periodic boundary conditions (periodic replication of the local cluster).– Each atom of the cluster now interacts also with the images of the other atoms of the cluster

$$W(𝐫)=\frac{1}{2}\underset{j𝒞}{}v(|𝐫𝐬_j|),$$
(8)
where $`𝐬_j`$ is the position of the image of atom $`j`$ selected through the minimum image criterion, i.e., among all possible images of atom $`j`$, only the closest is selected. We have used a cubic elementary cell.

iii) Icosahedral-like boundary conditions (hypothetical non-frustrated system).– As mentioned in the Introduction, icosahedral order cannot be extended to the entire space. One can however introduce icosahedral-like boundary conditions by embedding the $`13`$-atom cluster in the center of a large $`147`$-atom cluster with icosahedral $`(I_h)`$ symmetry in which the central atom and its first layer have been removed. (Recall that atomic clusters are characterized by “magic numbers” of atoms for which the global minimum has polytetrahedral symmetry northby87 ; harris84 ; in particular, icosahedral symmetry is obtained for $`13`$, $`55`$, and $`147`$ atoms.) The external potential is now written as

$$W(𝐫)=\frac{1}{2}\underset{kS_2^{ico},S_3^{ico}}{}v(|𝐫𝐫_k|),$$
(9)
where $`S_2^{ico}`$ and $`S_3^{ico}`$ are the second and third shells of the $`147`$-atom icosahedron. For completeness, we have also considered the $`55`$-atom Mackay icosahedron, as well as the correction due to embedding the $`55`$-atom or $`147`$-atom icosahedron in bulk liquid with the resulting external potential treated at a mean-field level (see above).

## III Global minimum of the Lennard-Jones clusters

We specialize our investigation to the case of the Lennard-Jones pair potential

$$v_{LJ}(r)=4ϵ\left[\left(\frac{\sigma }{r}\right)^{12}\left(\frac{\sigma }{r}\right)^6\right],$$
(10)
where $`ϵ`$ and $`2^{1/6}\sigma `$ are the well-depth and the separation at the minimum of the potential, respectively. In what follows we set $`\sigma =ϵ=1`$.

In order to evaluate the liquid-like external potential acting on the $`N`$-atom cluster, Eqs. (2) and (3), and its derivatives, Eqs. (4) and (5), one needs a model for the pair distribution function $`g(r)`$. We use the $`7`$-parameter parametrization of Verlet’s Molecular Dynamics simulation data on the Lennard-Jones liquid verlet68  proposed by Matteoli and Mansoori matteoli95 :

$$g(y)=1+y^m\left[g(d)1\lambda \right]+\left[\frac{(y1+\lambda )}{y}\right]\mathrm{exp}[\alpha (y1)]\mathrm{cos}[\beta (y1)],$$
(11)
for $`m1,y1`$, and

$$g(y)=g(d)\mathrm{exp}[\theta (y1)^2],$$
(12)
for $`y<1`$. Here $`y=r/d`$ is the dimensionless intermolecular distance where $`d=2^{1/6}\sigma `$, and $`h`$, $`m`$, $`\lambda `$, $`\alpha `$, $`\beta `$, $`\theta `$, $`g(d)`$ are adjustable parameters. The terms $`y^m`$ and $`\mathrm{exp}[\theta (y1)^2]`$ describe the decay of the first peak, while the term $`\mathrm{exp}[\alpha (y1)]\mathrm{cos}[\beta (y1)]`$ provides the damped oscillations observed at larger distances matteoli95 . We have taken the values $`h=1.065`$, $`m=13.42`$, $`g(d)=2.830`$, $`\lambda =0.9310`$, $`\alpha =1.579`$, $`\beta =6.886`$, $`\theta =135.9`$ that allow to reproduce the pair distribution function for the liquid at $`\rho =0.880`$ and $`T=1.095`$ (in usual reduced Lennard-Jones units) verlet68 ; matteoli95 .

The resulting external potential $`W(r;R_C)`$ is shown in Fig. 1 for several values of the cavity radius $`R_C`$. The shape of the $`r`$-dependence changes with $`R_C`$ so that no rescaling of the curves is possible. For $`R_C<1.2\sigma `$, the potential has a minimum at $`r=0`$ because the central atom sits at the minimum – or very close to it – of the pair interactions due to liquid atoms at the boundary of the cavity; this is no longer true for larger cavity radii, and $`W`$ decreases monotonically with $`r`$, the most favorable position inside the cavity being at its edge where the attractive interaction due to the nearby liquid particles is the strongest. By construction, when $`R_C0`$, $`W(r=0,R_C)`$ becomes equal to the total potential energy of the Lennard-Jones liquid at the considered state point, $`E5.7ϵ`$, whereas when $`R_C\mathrm{}`$, $`W(r=R_C,R_C)`$ is equal to half this energy.

As we discussed it above, the radius of the cavity $`R_C`$ must be adjusted to global contractions or expansions of the cluster. A reasonable way to implement this is to take at each minimization step, i. e., for each configuration of the $`N`$ atoms,

$$R_C=r_{max}+\mu \sigma ,$$
(13)
where $`r_{max}`$ is the distance of the outermost atom from the center of the cavity, and $`\mu `$ is a constant chosen to account for the fact that repulsive interactions between atoms make very unlikely the presence of “bulk liquid” atoms when their centers are too close to those of the cavity atoms; we have taken $`\mu =0.5`$, but we have checked that the results are independent of the actual value, in runs with different values of $`\mu `$ between $`0.1`$ and $`1.0`$.

A typical minimization run for a $`13`$-atom cluster in the presence of a mean-field liquid-like environment is shown in Fig. 2, where we have plotted the evolution of the energy and its intra-cluster and external-field contributions (bottom), together with the evolution of the cavity radius $`R_C`$ (top). During the optimization run that starts from a random configuration of $`13`$ atoms, the structure of the cluster becomes more and more compact, and its energy decreases. The final optimized configuration is found to have icosahedral symmetry. Both the final cluster radius, $`r_{max}1.08\sigma `$, and the final intra-cluster energy, $`U_{intra}44.327ϵ`$ are identical to those found for the ground-state of the isolated $`13`$-atom cluster northby87 . The total energy of the icosahedral cluster in the presence of a liquid-like environment is however much lower (by almost a factor two) because of the external field that compensates for the deficiency of nearest neighbors of the $`12`$ surface atoms. We have repeated the procedure for different starting random configurations and we have always obtained the icosahedron with $`r_{max}1.08\sigma `$ as the global minimum. The same is true for a whole range of liquid density $`\rho `$ and temperature $`T`$ ($`0.65\rho 0.88,0.6T3.6`$).

In addition to carrying out a global optimization procedure, one may also make a calculation in the spirit of Frank’s pioneering work frank52 . We have considered two potential candidates for the ground-state configuration of the $`13`$-atom cluster, namely the icosahedron and the cuboctahedral cluster with $`O_h`$ symmetry that is associated with the face-centered-cubic close-packed lattice and we have compared their energies in a liquid-like environment. The energies of the two, $`I_h`$-symmetric and $`O_h`$-symmetric, clusters are shown in Fig. 3 as a function of the cavity radius $`R_C`$ (or, equivalently, as a function of $`r_{max}`$, a unique distance being enough to fully determine the whole cluster once the symmetry, $`I_h`$ or $`O_h`$, is chosen.). The intra-cluster contribution to the energy has a minimum for $`r_{max}1.08\sigma `$ for the icosahedron, and $`r_{max}1.10\sigma `$ for the $`O_h`$ cluster, whereas the external-field contribution monotonously increases with $`R_C`$ in both cases. One then finds that the icosahedral cluster has a minimum total energy for $`r_{max}1.08\sigma `$ and that this energy is $`4.8\%`$ lower than the lowest energy found for the $`O_h`$ cluster when $`r_{max}1.10\sigma `$. When compared to Frank’s result $`(8.4\%)`$, the relative energy difference between the two types of clusters is thus not drastically modified: a mean-field liquid-like environment only slightly reduces the relative stability of the icosahedral order. Finally, as a mere check of our global optimization procedure, we have verified that the icosahedral ground-state found here is identical to that discussed above.

The influence of the boundary conditions on the ground-state of a $`13`$-atom cluster can be investigated by using again the optimization algorithm (always starting with a random initial configuration) with the appropriate conditions, free, periodic, and icosahedral-like, described in the previous Section. As already well-known, the ground-state of the isolated cluster is an icosahedron and, as anticipated, that of the cluster in the presence of icosahedral-like boundary conditions is also an icosahedron. The ground-state energies are shown in Table 1, and can be compared to that of the icosahedral cluster in a mean-field liquid-like environment. Not surprisingly, this latter is much lower than that of the isolated cluster (see above), but it is higher than that of icosahedra embedded in larger icosahedral structures.

We note on passing that, in the case of an icosahedral-like environment, the change in the structure of the $`13`$-atom cluster from random to icosahedral during the optimization run starts from the outside and propagates inward. This is reminiscent of what has been observed in the simulation of gold nanoclusters nam02 . There, it has been found that, just after freezing, ordered nanosurfaces with five-fold rotational symmetry are formed, while interior atoms remain in a disordered state. On lowering the temperature, the crystallization of the interior atoms proceeds from the surface toward the core region, eventually producing an icosahedral structure nam02 . This is at variance with the classical picture of homogeneous nucleation and rather represents a surface-induced (heterogeneous) crystallization.

Periodic boundary conditions lead to a quite different picture. Since icosahedra cannot tile space by periodic replication, such conditions should favor the symmetries that allow a complete filling of space with true long-range order. It is indeed what we have found: the global minimum is then a cuboctahedral cluster with $`O_h`$ symmetry, that leads to a face-centered-cubic close-packed lattice when periodically replicated (see Fig. 4); the corresponding ground-state energy is given in Table 1, and it is found lower than that of an isolated icosahedron — because of the lack of neighbors already mentioned for this latter case — but higher than icosahedra in either liquid-like or icosahedral-like environment.

Finally, we have considered the effect of varying the number of atoms present in the cluster. The main motivation for this study is to check that without a-priori knowledge of the preferred local structure in the presence of a mean-field liquid-like environment, hence of the number of atoms involved in this structure, the global minimization method will help select the proper preferred configuration. This will be important when considering molecular liquids. We have thus studied the global minima of $`N`$-atom clusters with $`N`$ ranging from $`2`$ to $`23`$. As shown in Fig. 5, the energy per-atom, the only relevant quantity for comparing local structures in a liquid-like environment, is lowest for the $`N=13`$ icosahedral ground-state. For small $`N`$, the mean-field description of the liquid environment is probably too crude to give sensible results, but it is nonetheless significant that for a large range of $`N`$, the icosahedral cluster is properly selected as the locally preferred structure of the Lennard-Jones liquid.

## IV Conclusion

Local icosahedral order has been found both in bulk condensed phases and in clusters formed by spherical particles. The connection between bulk and cluster studies is however obscured by two facts: first, geometrical frustration strongly hinders the spatial extension of local icosahedral order and distorts the local icosahedra in bulk conditions; second, the energetics of isolated clusters is partly determined by surface effects that are of course absent in the bulk. In this work we have tried to bypass those problems. We have proposed to determine the locally preferred structure of a liquid by finding the ground-state configuration of $`N`$-particle clusters embedded in a liquid-like continuum, characterized by the proper density and pair distribution function of the bulk liquid at the chosen thermodynamic state point. This mean-field-like procedure minimizes the surface effects without introducing full blown geometrical frustration.

In terms of potential energy surface, we have therefore introduced an effective energy surface that contains the usual intra-cluster potential energy contribution plus an external field that accounts for the interaction with the outside liquid at a mean-field level. By a global optimization algorithm we have then located the lowest-energy minimum of the effective energy surface. For Lennard-Jones pair interactions, we have found that the locally preferred structure is indeed an icosahedron, the effect of the liquid-like environment being to only slightly reduce the relative stability of the icosahedral structure when compared to Frank’s calculation for an isolated cluster frank52 . We have also shown the importance of the boundary conditions used for the cluster: whereas icosahedral-like boundary conditions stabilize even more the local icosahedral cluster, periodic boundary conditions make the cuboctahedral cluster more stable, a consequence of geometrical frustration that prevents tiling of space by icosahedra, and therefore favors long-range order associated with face-centered cubic or hexagonal close-packed lattices.

The present findings for the Lennard-Jones liquid suggest that the proposed method could be efficient as well for determining the locally preferred structure of molecular liquids, in cases where both translational and rotational degrees of freedom are involved, and a-priori knowledge about the putative local order is scarce.

###### Acknowledgements.

We thank D. J. Wales and F. Sciortino for fruitful comments.

*

## Appendix A Mean-field liquid external potential

Here we give some details about the calculation of the mean-field liquid external potential $`W(r;R_C)`$ discussed in Sect. II.

From Eq. (1) the external potential felt by an atom at $`𝐫`$, where the origin of the coordinates is chosen as the center of the cavity, can be written as

$$W(r,R_C)=\frac{\rho }{2}_{|𝐱+𝐫|>R_C}d^3xg(x)v(x),$$
(14)
where we have changed the integration variable from $`𝐫^{}`$ to $`𝐱=𝐫𝐫^{}`$. We now rotate the reference system such that $`𝐫/r=\widehat{𝐳}`$ and we translate it so that the new origin is at $`P`$, the position of the center of the atom under consideration. The new geometry of the problem is shown in Fig. 6. One can easily convince oneself that

$$W(r;R_C)=\frac{\rho }{2}_0^\pi 𝑑\varphi \mathrm{sin}\varphi _0^{2\pi }𝑑\theta _{s(\varphi )}^{\mathrm{}}𝑑xx^2g(x)v(x).$$
(15)
Given a point inside the cavity, $`s`$ is its distance from the surface in the direction $`𝐱`$. It is simple to check that $`s`$ is a function of $`\varphi `$ only and is solution of the quadratic equation

$$s^2+(2r\mathrm{cos}\varphi )s+(r^2R_C)^2=0;$$
(16)
the correct solution is

$$s(\varphi )=r\mathrm{cos}\varphi +\sqrt{r^2\mathrm{cos}^2\varphi +(R_C^2r^2)}.$$
(17)
Introducing the variable $`u=\mathrm{cos}\varphi `$ we obtain

$$W(r;R_C)=\pi \rho _1^1𝑑u_{s(u)}^{\mathrm{}}𝑑xx^2g(x)v(x),$$
(18)
with

$$s(u)=ru+\sqrt{r^2(u^21)+R_C^2}.$$
(19)
$`s(u)`$ is such that $`s(1)=R_C+r`$, and $`s(1)=R_Cr`$. We now change variable from $`u`$ to $`s`$, and integrate by parts so that

$`W(r;R_C)`$ $`=`$ $`\pi \rho {\displaystyle _{R_C+r}^{R_Cr}}𝑑s{\displaystyle \frac{du}{ds}}{\displaystyle _\rho ^{\mathrm{}}}𝑑xx^2g(x)v(x)`$ (20)
$`=`$ $`\pi \rho \left\{\left[u(s){\displaystyle _\rho ^{\mathrm{}}}𝑑xx^2g(x)v(x)\right]_{R_C+r}^{R_Cr}{\displaystyle _{R_Cr}^{R_C+r}}𝑑ss^2u(s)g(s)v(s)\right\},`$
where

$$u(s)=\frac{R_C^2r^2s^2}{2rs},$$
(21)
and $`u(s=R_Cr)=u(s=R_C+r)=1`$. Substituting Eq. (21) in Eq. (20) and rearranging we obtain Eqs. (2) and (3).