text:
# Strong-Electroweak Unification at About 4 TeV

## Abstract

I show how an $`SU(N)^M`$ quiver gauge theory can accommodate the standard model with three chiral families and unify all of $`SU(3)_C`$, $`SU(2)_L`$ and $`U(1)_Y`$ couplings with high accuracy at one unique scale estimated as $`M4`$ TeV.

Conformal invariance in two dimensions has had great success in comparison to several condensed matter systems. It is an interesting question whether conformal symmetry can have comparable success in a four-dimensional description of high-energy physics.

Even before the standard model (SM) $`SU(2)\times U(1)`$ electroweak theory was firmly established by experimental data, proposals were made  of models which would subsume it into a grand unified theory (GUT) including also the dynamics of QCD. Although the prediction of SU(5) in its minimal form for the proton lifetime has long ago been excluded, ad hoc variants thereof  remain viable. Low-energy supersymmetry improves the accuracy of unification of the three 321 couplings and such theories encompass a “desert” between the weak scale $`250`$ GeV and the much-higher GUT scale $`2\times 10^{16}`$ GeV, although minimal supersymmetric $`SU(5)`$ is by now ruled out.

Recent developments in string theory are suggestive of a different strategy for unification of electroweak theory with QCD. Both the desert and low-energy supersymmetry are abandoned. Instead, the standard $`SU(3)_C\times SU(2)_L\times U(1)_Y`$ gauge group is embedded in a semi-simple gauge group such as $`SU(3)^N`$ as suggested by gauge theories arising from compactification of the IIB superstring on an orbifold $`AdS_5\times S^5/\mathrm{\Gamma }`$ where $`\mathrm{\Gamma }`$ is the abelian finite group $`Z_N`$. In such nonsupersymmetric quiver gauge theories the unification of couplings happens not by logarithmic evolution over an enormous desert covering, say, a dozen orders of magnitude in energy scale. Instead the unification occurs abruptly at $`\mu =M`$ through the diagonal embeddings of 321 in $`SU(3)^N`$. The key prediction of such unification shifts from proton decay to additional particle content, in the present model at $`4`$ TeV.

Let me consider first the electroweak group which in the standard model is still un-unified as $`SU(2)\times U(1)`$. In the 331-model where this is extended to $`SU(3)\times U(1)`$ there appears a Landau pole at $`M4`$ TeV because that is the scale at which $`\mathrm{sin}^2\theta (\mu )`$ slides to the value $`\mathrm{sin}^2(M)=1/4`$. It is also the scale at which the custodial gauged $`SU(3)`$ is broken in the framework of .

Such theories involve only electroweak unification so to include QCD I examine the running of all three of the SM couplings with $`\mu `$ as explicated in e.g. . Taking the values at the Z-pole $`\alpha _Y(M_Z)=0.0101,\alpha _2(M_Z)=0.0338,\alpha _3(M_Z)=0.118\pm 0.003`$ (the errors in $`\alpha _Y(M_Z)`$ and $`\alpha _2(M_Z)`$ are less than 1%) they are taken to run between $`M_Z`$ and $`M`$ according to the SM equations

$`\alpha _Y^1(M)`$ $`=`$ $`(0.01014)^1(41/12\pi )\mathrm{ln}(M/M_Z)`$ (1)
$`=`$ $`98.6191.0876y`$ (2)
$`\alpha _2^1(M)`$ $`=`$ $`(0.0338)^1+(19/12\pi )\mathrm{ln}(M/M_Z)`$ (3)
$`=`$ $`29.586+0.504y`$ (4)
$`\alpha _3^1(M)`$ $`=`$ $`(0.118)^1+(7/2\pi )\mathrm{ln}(M/M_Z)`$ (5)
$`=`$ $`8.474+1.114y`$ (6)
where $`y=\mathrm{log}(M/M_Z)`$.

The scale at which $`\mathrm{sin}^2\theta (M)=\alpha _Y(M)/(\alpha _2(M)+\alpha _Y(M))`$ satisfies $`\mathrm{sin}^2\theta (M)=1/4`$ is found from Eqs.(2,4) to be $`M4`$ TeV as stated in the introduction above.

I now focus on the ratio $`R(M)\alpha _3(M)/\alpha _2(M)`$ using Eqs.(4,6). I find that $`R(M_Z)3.5`$ while $`R(M_3)=3`$, $`R(M_{5/2})=5/2`$ and $`R(M_2)=2`$ correspond to $`M_3,M_{5/2},M_2400\mathrm{G}\mathrm{e}\mathrm{V},4\mathrm{T}\mathrm{e}\mathrm{V},\mathrm{and}140\mathrm{T}\mathrm{e}\mathrm{V}`$ respectively. The proximity of $`M_{5/2}`$ and $`M`$, accurate to a few percent, suggests strong-electroweak unification at $`4`$ TeV.

There remains the question of embedding such unification in an $`SU(3)^N`$ of the type described in . Since the required embedding of $`SU(2)_L\times U(1)_Y`$ into an $`SU(3)`$ necessitates $`3\alpha _Y=\alpha _H`$ the ratios of couplings at $`4`$ TeV is: $`\alpha _{3C}:\alpha _{3W}:\alpha _{3H}::5:2:2`$ and it is natural to examine $`N=12`$ with diagonal embeddings of Color (C), Weak (W) and Hypercharge (H) in $`SU(3)^2,SU(3)^5,SU(3)^5`$ respectively.

To accomplish this I specify the embedding of $`\mathrm{\Gamma }=Z_{12}`$ in the global $`SU(4)`$ R-parity of the $`𝒩=4`$ supersymmetry of the underlying theory. Defining $`\alpha =\mathrm{exp}(2\pi i/12)`$ this specification can be made by $`\mathrm{𝟒}(\alpha ^{A_1},\alpha ^{A_2},\alpha ^{A_3},\alpha ^{A_4})`$ with $`\mathrm{\Sigma }A_\mu =0(\mathrm{mod12})`$ and all $`A_\mu 0`$ so that all four supersymmetries are broken from $`𝒩=4`$ to $`𝒩=0`$.

Having specified $`A_\mu `$ I calculate the content of complex scalars by investigating in $`SU(4)`$ the $`\mathrm{𝟔}(\alpha ^{a_1},\alpha ^{a_2},\alpha ^{a_3},\alpha ^{a_3},\alpha ^{a_2},\alpha ^{a_1})`$ with $`a_1=A_1+A_2,a_2=A_2+A_3,a_3=A_3+A_1`$ where all quantities are defined (mod 12).

Finally I identify the nodes (as C, W or H) on the dodecahedral quiver such that the complex scalars

$$\mathrm{\Sigma }_{i=1}^{i=3}\mathrm{\Sigma }_{\alpha =1}^{\alpha =12}(N_\alpha ,\overline{N}_{\alpha \pm a_i})$$
(7)
are adequate to allow the required symmetry breaking to the $`SU(3)^3`$ diagonal subgroup, and the chiral fermions

$$\mathrm{\Sigma }_{\mu =1}^{\mu =4}\mathrm{\Sigma }_{\alpha =1}^{\alpha =12}(N_\alpha ,\overline{N}_{\alpha +A_\mu })$$
(8)
can accommodate the three generations of quarks and leptons.

It is not trivial to accomplish all of these requirements so let me demonstrate by an explicit example.

For the embedding I take $`A_\mu =(1,2,3,6)`$ and for the quiver nodes take the ordering:

$$CWHCW^4H^4$$
(9)
with the two ends of (9) identified.

The scalars follow from $`a_i=(3,4,5)`$ and the scalars in Eq.(7)

$$\mathrm{\Sigma }_{i=1}^{i=3}\mathrm{\Sigma }_{\alpha =1}^{\alpha =12}(3_\alpha ,\overline{3}_{\alpha \pm a_i})$$
(10)
are sufficient to break to all diagonal subgroups as

$$SU(3)_C\times SU(3)_W\times SU(3)_H$$
(11)
The fermions follow from $`A_\mu `$ in Eq.(8) as

$$\mathrm{\Sigma }_{\mu =1}^{\mu =4}\mathrm{\Sigma }_{\alpha =1}^{\alpha =12}(3_\alpha ,\overline{3}_{\alpha +A_\mu })$$
(12)
and the particular dodecahedral quiver in (9) gives rise to exactly three chiral generations which transform under (11) as

$$3[(3,\overline{3},1)+(\overline{3},1,3)+(1,3,\overline{3})]$$
(13)
I note that anomaly freedom of the underlying superstring dictates that only the combination of states in Eq.(13) can survive. Thus, it is sufficient to examine one of the terms, say $`(3,\overline{3},1)`$. By drawing the quiver diagram indicated by Eq.(9) with the twelve nodes on a “clock-face” and using $`A_\mu =(1,2,3,6)`$ in Eq.(8) I find five $`(3,\overline{3},1)`$’s and two $`(\overline{3},3,1)`$’s implying three chiral families as stated in Eq.(13).

After further symmetry breaking at scale $`M`$ to $`SU(3)_C\times SU(2)_L\times U(1)_Y`$ the surviving chiral fermions are the quarks and leptons of the SM. The appearance of three families depends on both the identification of modes in (9) and on the embedding of $`\mathrm{\Gamma }SU(4)`$. The embedding must simultaneously give adequate scalars whose VEVs can break the symmetry spontaneously to (11). All of this is achieved successfully by the choices made. The three gauge couplings evolve according to Eqs.(2,4,6) for $`M_Z\mu M`$. For $`\mu M`$ the (equal) gauge couplings of $`SU(3)^{12}`$ do not run if, as conjectured in  there is a conformal fixed point at $`\mu =M`$.

The basis of the conjecture in  is the proposed duality of Maldacena which shows that in the $`N\mathrm{}`$ limit $`𝒩=4`$ supersymmetric $`SU(N)`$gauge theory, as well as orbifolded versions with $`𝒩=2,1`$ and $`0`$ become conformally invariant. It was known long ago  that the $`𝒩=4`$ theory is conformally invariant for all finite $`N2`$. This led to the conjecture in  that the $`𝒩=0`$ theories might be conformally invariant, at least in some case(s), for finite $`N`$. It should be emphasized that this conjecture cannot be checked purely within a perturbative framework. I assume that the local $`U(1)`$’s which arise in this scenario and which would lead to $`U(N)`$ gauge groups are non-dynamical, as suggested by Witten, leaving $`SU(N)`$’s.

This is a non-gravitational theory with conformal invariance when $`\mu >M`$ and where the Planck mass it taken to be infinitely large. The ubiquitous question is: What about gravity which breaks conformal symmetry in the ultraviolet (UV)? This is a question about the holographic principle for flat spacetime.

From the phenomenological viewpoint the equal couplings of $`SU(3)^{12}`$ can, instead of remaining constant at energies $`\mu >M`$, decrease smoothly by asymptotic freedom to a conformal fixed point as $`\mu \mathrm{}`$. This possibility is less restrictive and may fit in better with the AdS/CFT correspondence. The desert resides in the unexplored domain of the orders of magnitude in energy scale between 4 TeV and the gravitational scale, $`M_{Planck}`$.

As for experimental tests of such a TeV GUT, the situation at energies below 4 TeV is predicted to be the standard model with a Higgs boson still to be discovered at a mass predicted by radiative corrections  to be below 267 GeV at 99% confidence level.

There are many particles predicted at $`4`$ TeV beyond those of the minimal standard model. They include as spin-0 scalars the states of Eq.(10). and as spin-1/2 fermions the states of Eq.(12), Also predicted are gauge bosons to fill out the gauge groups of (11), and in the same energy region the gauge bosons to fill out all of $`SU(3)^{12}`$. All these extra particles are necessitated by the conformality constraints of  to lie close to the conformal fixed point.

One important issue is whether this proliferation of states at $`4`$ TeV is compatible with precision electroweak data in hand. This has been studied in the related model of  in a recent article. Those results are not easily translated to the present model but it is possible that such an analysis including limits on flavor-changing neutral currents could rule out the entire framework.

As alternative to $`SU(3)^{12}`$ another approach to TeV unification has as its group at $`4`$ TeV $`SU(6)^3`$ where one $`SU(6)`$ breaks diagonally to color while the other two $`SU(6)`$’s each break to $`SU(3)_{k=5}`$ where level $`k=5`$ characterizes irregular embedding. The triangular quiver $`CWH`$ with ends identified and $`A_\mu =(\alpha ,\alpha ,\alpha ,1)`$, $`\alpha =\mathrm{exp}(2\pi i/3)`$, preserves $`𝒩=1`$ supersymmetry. I have chosen to describe the $`𝒩=0`$ $`SU(3)^{12}`$ model in the text mainly because the symmetry breaking to the standard model is more transparent.

The TeV unification fits $`\mathrm{sin}^2\theta `$ and $`\alpha _3`$, predicts three families, and partially resolves the GUT hierarchy. If such unification holds in Nature there is a very rich level of physics one order of magnitude above presently accessible energy.

Is a hierarchy problem resolved in the present theory? In the non-gravitational limit $`M_{Planck}\mathrm{}`$ I have, above the weak scale, the new unification scale $`4`$ TeV. Thus, although not totally resolved, the GUT hierarchy is ameliorated. The gravitational hierarchy problem is not addressed.

My final remark is on the non-appearance of relevant deformations which break conformal invariance above 4 TeV. This is an assumption I make by analogy to several other systems in Nature with a large scaling region, e.g. superfluid helium where there is a comparable non-appearance of relevant operators over many orders of magnitude in scale size.

This work was supported in part by the Office of High Energy, US Department of Energy under Grant No. DE-FG02-97ER41036.


format: text

text:
# 1 Introduction

## 1 Introduction

The self-organizing map (SOM) algorithm (Kohonen 1982) served both as model for topology-preserving primary sensory processing in the cortex (Obermayer et al. 1992), and for technical applications (Ritter et al. 1992). Self-organizing feature maps map an input space, such as the retina or skin receptor fields, into a neural layer by feedforward structures with lateral inhibition. Defining properties are topology preservation, error tolerance, plasticity, and self-organized formation by a local process. Compared to other clustering algorithms and vector quantizers its apparent advantage for data visualization and exploration is its approximative topology preservation. In contrast to the Elastic Net (Durbin & Willshaw 1987) and the Linsker (1989) Algorithm, which are performing gradient descent in a certain energy landscape, the Kohonen algorithm lacks an energy function in the general case of a continuous input distribution. Although the learning process can be described in terms of a Fokker-Planck equation (Ritter & Schulten 1988), the expectation value of the learning step is a nonconservative force (Obermayer et al. 1992) driving the process so that it has no associated energy function. Despite a lot of research, the relationships between the Kohonen model and its variants to general principles remain an open field (Kohonen 1991).
To appear in Neural Computation

### 1.1 Kohonen’s Self Organizing Feature Map

Kohonen’s Self Organizing Map is defined as follows: Every stimulus $`𝐯`$ of an input space $`V`$ is mapped to a “center of excitation”, or winner

$`𝐬=\text{argmin}_{𝐫R}|𝐰_𝐫𝐯|,`$ (1)
where $`|.|`$ denotes the Euclidian distance in input space. In the Kohonen model, the learning rule for each synaptic weight vector $`𝐰_𝐫`$ is given by

$`\delta 𝐰_𝐫=\eta g_{\mathrm{𝐫𝐬}}(𝐯𝐰_𝐫),`$ (2)
where $`g_{\mathrm{𝐫𝐬}}`$ defines the neighborhood relation in $`R`$, and will throughout this paper be a Gaussian function of the Euclidian distance $`|𝐫𝐬|`$ in the neural layer. Topology preservation is enforced by the common update of all weight vectors whose neuron $`𝐫`$ is adjacent to the center of excitation $`𝐬`$; the adjacency function $`g_{\mathrm{𝐫𝐬}}`$ prescribes the topology in the neural layer. The speed of learning $`\eta `$ usually is decreased during the process.

### 1.2 The Winner Relaxing Kohonen Algorithm

We now consider an energy function $`V`$ first proposed in (Ritter et al. 1992). If we have a discrete input space, the potential function for the expectation value of the learning step is given by

$`V(\{𝐰\})={\displaystyle \frac{1}{2}}{\displaystyle \underset{\mathrm{𝐫𝐬}}{}}g_{\mathrm{𝐫𝐬}}^\gamma {\displaystyle \underset{\mu |𝐯^\mu F_𝐬(\{𝐰\})}{}}p(𝐯^\mu )|𝐯^\mu 𝐰_𝐫|^2,`$ (3)
where $`F_𝐬(\{𝐰\})`$ is the cell of the Voronoi tesselation (or Dirichlet tesselation) of input space defined by (1). For discrete input space, where $`p(𝐯)`$ is a sum over delta peaks $`\delta (𝐯𝐯^\mu )`$, the first derivative w.r.t. $`𝐰_𝐫`$ is not continuous at all weight vectors where the borders of the voronoi tesselation are shifting over one of the input vectors (Fig. 1).

However, (3) requires the assumption that none of the borders of the Voronoi tesselation is shifting over a pattern vector $`𝐯^\mu `$, which may be fulfilled in the final convergence phase for discrete input spaces, but becomes problematic if there are more receptor positions than neurons. If $`p(𝐯)`$ is continuous, the sum over $`\mu `$ becomes an integral, and with every stimulus vector update the surrounding Voronoi borders slide over stimuli (which means they become represented by annother weight vector), so that there is no global energy function for the general case.

We remark that replacing the crisp (or hard) winner selection (1) by a soft-winner $`𝐬=\text{argmin}_𝐫_𝐫^{}g_{\mathrm{𝐫𝐫}^{}}|𝐰_𝐫^{}𝐯|^2`$ minimizes (3) even in the continuous case (Graepel et al 1997, Heskes 1999). This is a formally elegant approach if one wants to ensure the existence of an energy function and accepts to modify the winner selection.

However, to motivate the Winner Relaxing learning, we return to the hard winner selection scheme (1) and take up the learning rule given by Kohonen (1991). Our use of this ansatz however is justified here only a posteriori by its use for adjusting the magnification.

From the shift of the borders of the Voronoi tesselation $`F_𝐬(\{𝐰\})`$ (see Fig. 1) in evaluation of the gradient with respect to a weight vector $`𝐰_𝐫`$, Kohonen (1991) derived for the (approximated) gradient descent in $`V`$ the additive term $`\frac{1}{2}\eta \delta _{\mathrm{𝐫𝐬}}_{𝐫^{^{}}𝐬}g_{𝐫^{^{}}𝐬}(𝐯𝐰_𝐫^{^{}})`$ extending (2) for the winning neuron. As it implied an additional elastic relaxation, it was straightforward to call it ‘Winner Relaxing’ (WR) Kohonen algorithm (Claussen 1992). In the remainder we study the (generalized) Winner Relaxing Kohonen algorithm, or Winner Relaxing Self-Organizing Map (WRSOM), introduced firstly in (Claussen 1992), in the form

$`\delta 𝐰_𝐫`$ $`=`$ $`\eta \{(𝐯𝐰_𝐫)g_{\mathrm{𝐫𝐬}}^\gamma \lambda \delta _{\mathrm{𝐫𝐬}}{\displaystyle \underset{𝐫^{^{}}𝐬}{}}g_{𝐫^{^{}}𝐬}^\gamma (𝐯𝐰_𝐫^{^{}})\},`$ (4)
where $`𝐬`$ is the center of excitation for incoming stimulus $`𝐯`$, and $`g_{\mathrm{𝐫𝐬}}^\gamma `$ is a Gaussian function of distance in the neural layer with characteristic length $`\gamma `$. Here $`\lambda `$ is a free parameter of the algorithm. The original algorithm (associated with the potential function) proposed by Kohonen in 1991 is obtained for $`\lambda =+1/2`$, whereas the classical Self Organizing Map Algorithm is obtained for $`\lambda =0`$. The influence of $`\lambda `$ on the magnification behaviour is the central issue of this paper.

### 1.3 The Magnification Factor

The magnification factor is defined as the density of neurons $`𝐫`$ (i.e. the density of synaptic weight vectors $`𝐰_𝐫`$) per unit volume of input space, and therefore is given by the inverse Jacobi determinant of the mapping from input space to neuron layer: $`M=|J|^1=|det(d𝐰/d𝐫)|^1`$. We assume the input space to be continuous and of same dimension as the neural layer, and the map to be noninverting $`(J>0)`$.

The magnification factor quantifies the networks’ response to a given probability density of stimuli $`P(𝐯)`$. To evaluate $`M`$ in higher dimensions, one in general has to compute the equilibrium state of the whole network and needs therefore the complete global knowledge on $`P(𝐯)`$, except for separable cases. For one-dimensional mappings the magnification factor can follow an universal magnification law, that is, $`M(\overline{𝐰}(𝐫))`$ is a function of the local probability density $`P`$ only, independent of both the location $`𝐫`$ in the neural layer and the location $`\overline{𝐰}(𝐫)`$ in input space. Hereby it is nontrivial whether there exists a power law or not; the Elastic Net obeys an universal magnification law that remarkably is not a power law (Claussen & Schuster 2002) due to a nonvanishing elastic tension in regions of small input density. For the classical Kohonen algorithm the magnification law is given by a power law $`M(\overline{𝐰}(𝐫))P(\overline{𝐰}(𝐫))^\rho `$ with exponent $`\rho =\frac{2}{3}`$ (Ritter & Schulten 1986). See Table 1 for an overview. For a discrete neural layer and different neighborhood kernels corrections apply (Ritter 1991, Ritter et al. 1992, Dersch& Tavan 1995).

As the brain is assumed to be optimized by evolution for information processing, one could conjecture that maximal mutual information can define an extremal principle governing the setup of neural structures. For feedforward neural structures with lateral inhibition, an algorithm of maximal mutual information has been defined by Linsker (1989) using the gradient descend in mutual information. It requires computationally costly integrations, and has a highly nonlocal learning rule; therefore it is neither favourable as a model for biological maps, nor feasible for technical applications. Due to realization constraints, both technical applications and cortical networks (Plumbley 1999) are not necessarily capable of reaching this optimum. Even if one had experimental data of the magnification behaviour, the question from what self-organizing dynamics neural structures emerge, remains. Overall it is desirable to find learning rules that minimize mutual information in a simpler way.

An optimal map from the view of information theory would reproduce the input probability exactly ($`MP(𝐯)^\rho `$ with $`\rho =1`$), being equivalent to the condition that all neurons in the layer are firing with same probability. This defines an equiprobabilistic mapping (van Hulle 2000). An exponent $`\rho =0`$, on the other hand, corresponds to a uniform distribution of weight vectors, or no adaptation at all. So the magnification exponent is a direct indicator, how far a Self Organizing Map algorithm is away from the optimum predicted by information theory.

## 2 Magnification Exponent of the Winner-Relaxing Kohonen Algorithm

We now derive the magnification law of the Winner-Relaxing Kohonen algorithm (4) for the case of a 1D$``$1D map. Note that for higher dimensions analytical results can only be obtained for special degenerate cases of the input probability density and therefore lack generality.

The necessary condition for the final state of the algorithm is that the expectation value of the learning step vanishes for all neurons $`r`$:

$`_{rR}\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}={\displaystyle dvp(v)\delta \overline{w}_r(v)}.`$ (5)
Since this expectation value is equal to the learning step of the pattern parallel rule, (5) is the stationary state condition for both serial and parallel updating, and also for batch updating. Thus we can proceed for these variants simultaneously (As synaptic plasticity is widely assumed to be based on integrative effects, one could claim that a parallel model is sufficient). The update rule (4) can be extended by an additional diagonal term controlled by $`\mu `$111 Whereas the extra term controlled by the parameter $`\mu `$ has been introduced in (Claussen 1992) for pure generality, and will be kept within the derivation, it does not contribute to the magnification. In general, the setting $`\mu =0`$, is recommended (and probably most stable), and the Winner Relaxing Kohonen algorithm thus has only one relevant control parameter $`\lambda `$.:

$`\delta w_r`$ $`=`$ $`\eta \{(vw_r)g_{rs}^\gamma +\mu (vw_r)\delta _{rs}`$ (6)
$`\lambda \delta _{rs}{\displaystyle \underset{r^{^{}}s}{}}g_{r^{^{}}s}^\gamma (vw_r^{^{}})\}.`$
By insertion of the update rule (6) one obtains

$`0`$ $`=`$ $`{\displaystyle dsP(\overline{w}(s))J(s)g_{rs}^\gamma (\overline{w}(s)\overline{w}(r))}`$ (7)
$`+(\mu +\lambda )\underset{0}{\underset{}{{\displaystyle dsP(\overline{w}(s))J(s)\delta _{rs}(\overline{w}(s)\overline{w}(r))}}}`$
$`\lambda {\displaystyle dsdr^{^{}}P(\overline{w}(s))J(s)\delta _{rs}g_{r^{^{}}s}^\gamma (\overline{w}(s)\overline{w}(r^{^{}}))}`$
$`=`$ $`{\displaystyle dsP(\overline{w}(s))J(s)g_{rs}^\gamma (\overline{w}(s)\overline{w}(r))}`$
$`+\lambda P(\overline{w}(r))J(r){\displaystyle dr^{^{}}g_{r^{^{}}r}^\gamma (\overline{w}(r^{^{}})\overline{w}(r))}.`$
The derivation can be performed analoguous to (Ritter, Martinetz & Schulten 1992). In the continuum limit there is always an exactly matching winning weight vector $`\overline{w}_s=v`$. Further the integration variable is substituted, $`\mathrm{d}v=\mathrm{d}\overline{w}_s=J(s)\mathrm{d}s`$, and we define the abbreviation $`\overline{P}:=P(\overline{w}(r))`$. In the first integrand $`\overline{P}J`$ has to be expanded in powers of $`q:=sr`$. Within the second integral $`\overline{P}J`$ is evaluated only at $`r`$. Thus the integration yields in leading order in $`q`$:

$`0`$ $`=`$ $`\gamma ^2\left({\displaystyle \frac{\mathrm{d}\overline{w}}{\mathrm{d}r}}{\displaystyle \frac{\mathrm{d}(\overline{P}J)}{\mathrm{d}r}}+{\displaystyle \frac{1}{2}}\overline{P}J{\displaystyle \frac{\mathrm{d}^2\overline{w}}{\mathrm{d}r^2}}\right)`$
$`+\lambda \overline{P}J{\displaystyle dqg_{0q}^\gamma (\underset{\text{contribution 0}}{\underset{}{q\frac{\mathrm{d}\overline{w}}{\mathrm{d}r}}}+\frac{q^2}{2}\frac{\mathrm{d}^2\overline{w}}{\mathrm{d}r^2})}`$
$`0`$ $`=`$ $`\gamma ^2\left(J{\displaystyle \frac{\mathrm{d}(\overline{P}J)}{\mathrm{d}r}}+{\displaystyle \frac{\overline{P}J}{2}}{\displaystyle \frac{\mathrm{d}J}{\mathrm{d}r}}+\lambda {\displaystyle \frac{\overline{P}J}{2}}{\displaystyle \frac{\mathrm{d}J}{\mathrm{d}r}}\right).`$ (8)
Further we have to require $`\gamma 0,`$ $`P0,`$ $`d\overline{P}/dr0`$. Then the ansatz of an universal local magnification law $`J(r)=J(\overline{P}(r))`$, i.e. $`J`$ depends only on the local value of $`P`$, that may be expected for the one-dimensional case only, requires $`J`$ to fulfill the differential equation

$`0`$ $`=`$ $`{\displaystyle \frac{J}{\overline{P}}}+(1+{\displaystyle \frac{1}{2}}+{\displaystyle \frac{\lambda }{2}}){\displaystyle \frac{\mathrm{d}J}{\mathrm{d}\overline{P}}}`$ (9)
or

$`{\displaystyle \frac{\mathrm{d}J}{\mathrm{d}\overline{P}}}`$ $`=`$ $`{\displaystyle \frac{2}{3+\lambda }}{\displaystyle \frac{J}{\overline{P}}}.`$ (10)
It has a power law solution, (provided that $`\lambda 3`$), which verifies the ansatz made above, $`J`$ being a function of the local density only,

$`M={\displaystyle \frac{1}{J}}P(v)^{{\scriptscriptstyle \frac{2}{3+\lambda }}}.`$ (11)
Thus the magnification exponent is given by $`\frac{2}{3+\lambda }`$ and can be tuned from $`1/2`$ to 1 (see Fig. 2) within the range of stability.

For the $`\lambda =1/2`$ choice of the Winner-Relaxing Kohonen Algorithm the magnification factor follows an exact power law with magnification exponent $`\rho =4/7`$, which is smaller than $`\rho =2/3`$ for the classical Self Organizing Feature Map (Ritter & Schulten 1986), but is still much larger than $`\rho =1/3`$ for Vector Quantization and Neural Gas. In any case, the maps resulting from the choices $`\lambda =1/2`$ and $`\lambda =0`$ are not optimal in terms of information theory.

## 3 Enhancing the Magnification

From this result one would try to invert the Relaxing Effect by choice of negative values for $`\lambda `$, which means to “enforce” the winner. In fact, the choice of $`\lambda =1`$ leads to the magnification exponent 1.

The magnification law (11) is verified numerically as is shown in Table 2. Apart from the fact that the exponent can be varied by a priori parameter choice between $`1/2`$ and $`1`$, the simulations show that our Winner Relaxing Algorithm is able to establish information-theoretically optimal self-organizing maps in the “winner enforcing” case ($`\lambda <0`$).

## 4 Ordering time and Stability Region

At least for the 2D$``$2D case, the Winner-Relaxing Kohonen Algorithm was reported as ‘somewhat faster’ (Kohonen 1991) in the initial ordering process. In a 1D$``$1D sample setup (Claussen 2003), a marginally quicker ordering was observed for negative $`\lambda `$, at least at a relatively high learning rate $`\eta =1`$. As a lot of parameters and the input distribution itself influence the ordering time and decay of fluctuations, different results may be obtained; e.g., a small fraction of input distributions containing topological kinks take much longer to become ordererd, thus minimal, maximal, averaged, and inverse averaged ordering time will deviate.

If one instead investigates the time dependence of the fluctuations, for positive values of $`\lambda `$ a considerably quicker decay is observed (Fig. 4), being consistent with the observation by Kohonen (1991) mentioned above.

These simulations indicate that for obtaining optimal magnification, the price of a longer learning phase may have to be paid. However, this drawback can be circumvented by combining the advantages of both $`\lambda `$ ranges; i.e. using $`\lambda >1`$ in the initial phase to speed up ordering, and switching to $`\lambda =1`$ after a considerable decay of fluctuations (Fig. 4). No complicated time-dependence of this parameter switch has been used, and neither learning rate nor neighborhood have been changed during the simulation.

The last important issue to be addressed is the dependence of stability on the parameter $`\lambda `$, especially at the border $`1`$. Fortunately, the algorithm appears to be stable (in the 1D$``$1D case) in the whole range $`1\lambda +1`$, as shown in Fig. 5. On both borders the Winner Relaxing learning remains stable. Thus, the full range of magnification exponents between $`1/2`$ and $`1`$ can be acheived.

In higher dimensions no universal magnification law is expected, but one can evaluate the output entropy for a given input distribution and network. As shown in Fig. 6, the enhancement of output entropy by Winner Relaxing learning is effective also in the twodimensional case, where however parameters have to be chosen more carefully.

## 5 Discussion

After our first study (Claussen 1992), Herrmann et al. (1995) introduced annother modification of the learning process, which was also applied to the Neural Gas algorithm (Villmann & Herrmann 1998). Their central idea is to use a learning rate $`\eta `$ being locally dependent on the input probability density and also an exponent 1 can be obtained. As the input probability density should not be available to a neural map that self-organizes from stimuli drawn from that distribution, it is estimated from the actual local reconstruction mismatch (being an estimate for the size of the Voronoi cell) and from the time elapsed since the last time being the winner. Both operations require additional memory and computation, and, due to the estimating character, the learning rate has to be bounded in practical use. This localized learning was overall easier applicable and overcame the stability problems of the early approach of conscience learning (deSieno 1988).

Another systematic method, the extended Maximum Entropy Learning Rule, has been introduced by van Hulle (1997). It approximates a map of maximal output entropy for arbitrary dimension, alhough in higher dimensions the handling of the quantization regions becomes less practial (van Hulle 1998). A quite different approach being also capable of generating equiprobabilistic maps is via kernel optimization (van Hulle 1998, 2000, 2002), i.e. neighborhood kernel radii themselves become learning parameters, in addition to the weight vectors defining the kernel centers. Other approaches, also influencing magnification, consider the selection of the winner to be probabilistic, leading to elegant statistical approaches to potential functions, as given by Graepel et al. (1997) and Heskes (1999).

As shown recently (Claussen & Villmann 2004), the Winner Relaxing concept can also be transferred successfully to the Neural Gas, confirming the utility of this class of learning rules.

## 6 Conclusions

The Linsker, Elastic Net and Winner-Relaxing Kohonen algorithms can be derived from an extremal principle, given by information theory, physical motivations, and reconstruction error, respectively. In this paper we have chosen the magnification law to indicate how close the algorithm reaches the adaptation properties of a map of maximal mutual information. The magnification law is one quantitative property that both is accessible by neurobiological experiments and manifests as a quantitative control parameter of a neural map used as vector quantizer in applications. A map of maximal mutual information uses all neurons with same probability, i.e. their firing rate will be equal.

In this work we have investigated the Winner Relaxing approach to establish a new family of vector quantizers. The shift from Kohonen ($`\rho =2/3`$) to Winner Relaxing Kohonen algorithm ($`\rho =4/7`$) seems to be marginal, if the emphasis is laid on the existence of a potential function. If a large magnification exponent is desired, the Winner Relaxing Kohonen Algorithm (with $`\lambda =1`$) combines simple computation with a magnification corresponding to maximal mutual information.

Acknowledgements: The author wants to thank H.G. Schuster for raising attention to the topic and for stimulating discussions.

## References

Claussen, J.C. (born Gruel) (1992). Selbstorganisation Neuronaler Karten. Diploma thesis, Kiel, Germany.

Claussen, J.C. & Schuster, H.G. (2002). Asymptotic level density of the Elastic Net Self-Organizing Map, Proc. ICANN 2002.

Claussen, J.C. (2003). Winner Relaxing and Winner-Enhancing Kohonen Maps: Maximal Mutual information from Enhancing the Winner. Complexity, 8(4), 15-22.

Claussen, J.C. & Villmann, T. (2004). Magnification Control in Winner Relaxing Neural Gas. Neurocomputing, in print.

Dersch, D.R., & Tavan, P. (1995). Asymptotic level density in topological feature maps. IEEE Transactions on Neural Networks, 6, 230-236.

DeSieno, D. (1988). Adding a conscience to competitive learning. In: Proc. ICNN’88, International Conference on Neural Networks, pp. 117-124, IEEE Service Center, Piscataway, NJ.

Durbin, R. & Willshaw, D. (1987). An analogue approach to the Travelling Salesman Problem using an Elastic Net Method. Nature, 326, 689-691.

Herrmann, M., Bauer, H.-U., Der, R. (1995). Optimal Magnification Factors in Self-Organizing Maps. Proc. ICANN 1995, pp. 75-80.

Graepel, T., Burger, M., & Obermayer, K. (1997). Phase transitions in stochastic self-organizing maps. Phys. Rev. E, 56, 3876-3890.

Heskes, T. (1999). Energy functions for self-organizing maps. In: Oja and Karski, ed.: Kohonen Maps, Elsevier.

van Hulle, M. M. (1997). Nonparametric density estimation and regression achieved with topographic maps maximizing the information-theoretic entropy of their outputs. Biological Cybernetics, 77, 49-61.

van Hulle, M. M. (1998). Kernel-based equiprobabilistic topographic map formation. Neural Computation. 10, 1847-1871.

van Hulle, M. M. (2000). Faithful Representations and Topographic Maps. Wiley.

van Hulle, M. M. (2002). Kernel-based topographic map formation by local density modeling. Neural Computation, 14, 1561-1573.

Kohonen, T. (1982). Self-Organized Formation of Toplogically Correct Feature Maps. Biological Cybernetics, 43, 59-69.

Kohonen, T. (1991). Self-Organizing Maps: Optimization Approaches. In: Artificial Neural Networks, ed. T. Kohonen et al. (North-Holland, Amsterdam).

Linsker, R. (1989). How To generate Ordered maps by Maximizing the Mutual Information between Input and Output Signals. Neural Computation, 1, 402-411.

Obermayer, K., Blasdel, G.G., & Schulten, K. (1992). Statistical-mechanical analysis of self-organization and pattern formation during the development of visual maps. Phys. Rev. A, 45, 7568-7589.

Plumbley, M.D. (1999). Do cortical maps adapt to optimize information density?, Network, 10, 41-58.

Ritter, H. & Schulten, K. (1986). On the Stationary State of Kohonen’s Self-Organizing Sensory Mapping. Biol. Cybernetics, 54, 99-106.

Ritter, H. & Schulten, K. (1988). Convergence Properties of Kohonen’s Topology Conserving Maps: Fluctuations, Stability and Dimension Selection. Biological Cybernetics, 60, 59-71.

Ritter, H. (1991). Asymptotic Level Density for a Class of Vector Quantization Processes. IEEE Trans. Neural Networks, 2, 173-175.

Ritter, H., Martinetz, T., & Schulten, K. (1992). Neural Computation and Self-Organizing Maps: An Introduction. Addison-Wesley, NY.

Simic, P.D. (1990). Statistical mechanics as the underlying theory of ‘elastic’ and ‘neural’ optimizations. Network, 1, 89-103.

Villmann, T., & Herrmann, M. (1998). Magnification Control in Neural Maps, Proc. ESANN 1998, pp. 191-196.

Manuscript submitted August 15, 2002; revised June 4 & September 24, 2004; accepted November 1, 2004.


text:
# The Hamiltonian Mean Field Model: from Dynamics to Statistical Mechanics and back

## 1 Introduction

Long-range interactions appear in the domains of gravity paddyhmf ; chavanishmf  and of plasma physics elskenshmf  and make the statistical treatment extremely complex. Additional features are present in such systems at short distances: the gravitational potential is singular at the origin and screening phenomena mask the Coulomb singularity in a plasma. This justifies the introduction of simplified toy models that retain only the long-range properties of the force, allowing a detailed description of the statistical and dynamical behaviors associated to this feature. In this context a special role is played by mean-field models, i.e. models where all particles interact with the same strength. This constitutes a dramatic reduction of complexity, since in such models the spatial coordinates have no role, since each particle is equivalent. However, there are several preliminary indications that behaviors found in mean-field models extend to cases where the two-body potential decays at large distances with a power smaller than space dimension tamahmf ; BMRhmf ; campa .

The Blume-Emery-Griffiths (BEG) mean-field model is discussed in this book BMRhmf  and represents an excellent benchmark to discuss relations between canonical and microcanonical ensembles. Indeed, this model is exactly solvable in both ensembles and is, at the same time, sufficiently rich to display such interesting features as negative specific heat and temperature jumps in the microcanonical ensemble. Since these effects cannot be present in the canonical ensemble, this rigorously proves ensemble inequivalence. However, the BEG model has no dynamics and only the thermodynamical behavior can be investigated. Moreover, it is a spin model where variables take discrete values. It would therefore be amenable to introduce a model that displays all these interesting thermodynamical effects, but for which one would also dispose of an Hamiltonian dynamics with continuous variables, whose equilibrium states could be studied both in the canonical and in the microcanonical ensemble. Having access to dynamics, one could moreover study non equilibrium features and aspects of the microscopic behavior like sensitivity to initial conditions, expressed by the Lyapunov spectrum Ruelleeckman . Such a model has been introduced in Ref. Antoni  and has been called the Hamiltonian Mean Field (HMF) model. In the simpler version, it represents a system of particles moving on a circle, all coupled by an equal strength attractive or repulsive cosine interaction. An extension of it to the case in which particles move on a 2D torus has been introduced in Ref. at  and it has been quite recently realized that all such models are particular cases of a more general Hamiltonian art .

The HMF model, that we introduce in Section 2, is exactly solvable in the canonical ensemble by a Hubbard-Stratonovich transformation. The solution in the microcanonical ensemble can be obtained only under certain hypotheses that we will discuss in Section 3, but detailed information on the behavior in the microcanonical ensemble can be obtained by direct molecular dynamics (MD) simulations. The model has first and second order phase transitions and tricritical points. Its rich phase diagram allows to test the presence of ensemble inequivalence near canonical first order phase transitions and, indeed, we find negative specific heat and temperature jumps in the microcanonical ensemble. Having access to dynamics, one can study metastability of out-of-equilibrium states. This is done in Section 4, where we analyze the emergence of a coherent structure in the repulsive HMF at low energy. Similar features are also discussed in another chapter of this book tsallisrap  for the attractive case near the second order phase transition. Section 5 is devoted to the study of the spectrum of Lyapunov exponents. The maximal exponent has a peak near the phase transition at ; lat1hmf ; celiahmf  and vanishes when increasing the number of particles with a universal scaling law in the whole high energy disordered phase. In a low energy range the Lyapunov spectrum has a thermodynamic limit distribution similar to the one observed for systems with short-range interaction livipolruffo .

## 2 The HMF models

A generic two-body potential in a two dimensional square box of side $`2\pi `$ with periodic boundary conditions, a 2D torus, can be Fourier expanded as

$$V(x,y)=\underset{𝐤=(k_x,k_y)}{}\mathrm{exp}\left(i𝐤𝐫\right)V(k_x,k_y).$$
(1)
A sufficiently rich family of potential functions is obtained if, we restrict to the first two momentum shells $`|k|=1`$ and $`|k|=\sqrt{2}`$, we require that the potential is only invariant under discrete rotations by all multiples of $`\pi /4`$, and we assume that the Fourier coefficients on each shell are the same. This amounts to perform a truncation in the Fourier expansion of the potential (1), as done in studies of spherically symmetric gravitational systems in another chapter of this book grossdh . We get

$$V(x,y)=a+b\mathrm{cos}x+b\mathrm{cos}y+c\mathrm{cos}x\mathrm{cos}y.$$
(2)
As the constant $`a`$ is arbitrary and scaling $`b`$ is equivalent to scale the energy, $`c`$ remains the only free parameter. We consider $`N`$ particles interacting through the two-body potential $`V(x,y)`$ and we adopt the Kac prescription Kachmf  to scale the equal strength coupling among the particles by their number $`N`$. This scaling allows to perform safely the thermodynamic limit, since both the kinetic and the potential energy increase proportional to $`N`$111Kac prescription is, however, unphysical and it would be important to find a viable alternative.. By appropriately redefining the constants $`a=2\epsilon +A,b=\epsilon ,c=A`$ in formula (2) and using Kac prescription one gets the following potential energy

$`V_A`$ $`=`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}\epsilon \left(1\mathrm{cos}(x_ix_j)\right)+\epsilon \left(1\mathrm{cos}(y_iy_j)\right)`$ (4)
$`+A\left(1\mathrm{cos}(x_ix_j)\mathrm{cos}(y_iy_j)\right),`$
with $`(x_i,y_i)]\pi ,\pi ]\times ]\pi ,\pi ]`$ representing the coordinates of $`i`$-th particle and $`(p_{x,i},p_{y,i})`$ its conjugated momentum. The Hamiltonian of the HMF model is now the sum of this potential energy with the kinetic energy

$$K=\underset{i=1}{\overset{N}{}}\left(\frac{p_{x,i}^2+p_{y,i}^2}{2}\right).$$
(5)
We get

$$H_{HMF}=K+V_A.$$
(6)
In the following, we will consider model (6) for $`A=0`$ and both $`\epsilon `$ positive (attractive case) and negative (repulsive case). The $`A0`$ case will always have $`A>0`$ and $`\epsilon =1`$.

## 3 Equilibrium Thermodynamics

In this section, we discuss the equilibrium thermodynamical results for model (4) in the canonical and microcanonical ensembles. Canonical results will be obtained analytically while, for the microcanonical ones, we will mostly rely on molecular dynamics (MD) simulations.

### 3.1 Canonical Ensemble for $`A=0`$

For pedagogical reasons, we will initially limit our analysis to the case $`A=0`$, for which the model reduces to two identical uncoupled systems: one describing the evolution of the $`\{x_i,p_{x,i}\}`$ variables and the other $`\{y_i,p_{y,i}\}`$. Therefore let us rewrite the Hamiltonian associated to one of these two sets of variables, named $`\theta _i`$ in the following. We obtain

$$H_0=\underset{i=1}{\overset{N}{}}\frac{p_i^2}{2}+\frac{\epsilon }{2N}\underset{i,j=1}{\overset{N}{}}[1\mathrm{cos}(\theta _i\theta _j)]=K_0+V_0$$
(7)
where $`\theta _i[\pi ;\pi [`$ and $`p_i`$ are the corresponding momenta. This model can be seen as representing particles moving on the unit circle, or as classical $`XY`$-rotors with infinite range couplings. For $`\epsilon >0`$, particles attract each other and rotors tend to align (ferromagnetic case), while for $`\epsilon <0`$, particles repel each other and spins tend to anti-align (antiferromagnetic case). At short distances, we can either think that particles cross each other or that they collide elastically since they have the same mass.

The physical meaning of this model is even clearer if one introduces the mean field vector

$$𝐌=M\mathrm{e}^{i\varphi }=\frac{1}{N}\underset{i=1}{\overset{N}{}}𝐦_i$$
(8)
where $`𝐦_i=(\mathrm{cos}\theta _i,\mathrm{sin}\theta _i)`$. $`M`$ and $`\varphi `$ represent the modulus and the phase of the order parameter, which specifies the degree of clustering in the particle interpretation, while it is the magnetization for the $`XY`$ rotors. Employing this quantity, the potential energy can be rewritten as a sum of single particle potentials $`v_i`$

$$V_0=\frac{1}{2}\underset{i=1}{\overset{N}{}}v_i\mathrm{with}v_i=1M\mathrm{cos}(\theta _i\varphi ).$$
(9)
It should be noticed that the motion of each particle is coupled to all the others, since the mean-field variables $`M`$ and $`\varphi `$ are determined at each time $`t`$ by the instantaneous positions of all particles.

The equilibrium results in the canonical ensemble can be obtained from the evaluation of the partition function

$$Z=d^Np_id^N\theta _i\mathrm{exp}(\beta H)$$
(10)
where $`\beta =1/(k_BT)`$, with $`k_B`$ the Boltzmann constant and $`T`$ the temperature. The integration domain is extended to the whole phase space. Integrating over momenta, one gets:

$$Z=\left(\frac{2\pi }{\beta }\right)^{N/2}_\pi ^\pi d^N\theta _i\mathrm{exp}\left[\frac{\beta \epsilon N}{2}(1𝐌^2)\right].$$
(11)
In order to evaluate this integral, we use the two dimensional Gaussian identity

$$\mathrm{exp}\left[\frac{\mu }{2}𝐱^\mathrm{𝟐}\right]=\frac{1}{\pi }_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}𝑑𝐲\mathrm{exp}[𝐲^2+\sqrt{2\mu }𝐱𝐲]$$
(12)
where $`𝐱`$ and $`𝐲`$ are two-dimensional vectors and $`\mu `$ is positive. We can therefore rewrite Eq. (11) as

$$Z=\left(\frac{2\pi }{\beta }\right)^{N/2}\mathrm{exp}\left[\frac{\beta \epsilon N}{2}\right]J$$
(13)
with

$$J=\frac{1}{\pi }_\pi ^\pi d^N\theta _i_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}𝑑𝐲\mathrm{exp}[𝐲^2+\sqrt{2\mu }𝐌𝐲]$$
(14)
and $`\mu =\beta \epsilon N`$. We use now definition (8) and exchange the order of the integrals in (14), factorizing the integration over the coordinates of the particles. Introducing the rescaled variable $`𝐲𝐲\sqrt{N/2\beta \epsilon }`$, one ends up with the following expression for $`J`$

$$J=\frac{N}{2\pi \beta \epsilon }_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}𝑑𝐲\mathrm{exp}\left[N\left(\frac{y^2}{2\beta \epsilon }\mathrm{ln}\left(2\pi I_0(y)\right)\right)\right]$$
(15)
where $`I_n`$ is the modified Bessel function of order $`n`$ and $`y`$ is the modulus of $`𝐲`$. Finally, integral (15) can be evaluated by employing the saddle point technique in the mean-field limit (i.e. for $`N\mathrm{}`$). In this limit, the Helmholtz free energy per particle $`f`$ reads as :

$$\beta f=\underset{N\mathrm{}}{lim}\frac{\mathrm{ln}Z}{N}=\frac{1}{2}\mathrm{ln}\left(\frac{2\pi }{\beta }\right)+\frac{\epsilon \beta }{2}+\underset{y}{\mathrm{max}}\left(\frac{y^2}{2\beta \epsilon }\mathrm{ln}(2\pi I_0(y))\right).$$
(16)
The maximum condition leads to the consistency equation

$$\frac{y}{\beta \epsilon }=\frac{I_1(y)}{I_0(y)}.$$
(17)
For $`\epsilon <0`$, there is a unique solution $`\overline{y}=0`$, which means that the order parameter remains zero and there is no phase transition (see Figs. 1c,d). The particles are all the time homogeneously distributed on the circle and the rotors have zero magnetization. On the contrary, in the ferromagnetic case ($`\epsilon >0`$), the solution $`\overline{y}=0`$ is unstable for $`\beta \beta _c=2`$. At $`\beta =\beta _c`$, two stable symmetric solutions appear through a pitchfork bifurcation and a discontinuity in the second derivative of the free energy is present, indicating a second order phase transition.

These results are confirmed by an analysis of the order parameter222This is obtained by adding to the Hamiltonian an external field and taking the derivative of the free energy with respect to this field, evaluated at zero field.

$$M=\frac{I_1(\overline{y})}{I_0(\overline{y})}.$$
(18)
For $`\epsilon >0`$, the magnetization $`M`$ vanishes continuously at $`\beta _c`$ (see Fig. 1a,b), while it is always identical to zero in the antiferromagnetic case (see Fig. 1c,d). Since $`M`$ measures the degree of clustering of the particles, we have for $`\epsilon >0`$ a transition from a clustered phase when $`\beta >\beta _c`$ to a homogeneous phase when $`\beta <\beta _c`$. We can obtain also the energy per particle

$$U=\frac{(\beta f)}{\beta }=\frac{1}{2\beta }+\frac{\epsilon }{2}\left(1M^2\right)$$
(19)
which is reported for $`|\epsilon |=1`$ in Fig. 1. Panels (c) and (d) of Fig. 1 are limited to the range $`U>0`$ because in the antiferromagnetic model a non-homogeneous state, a bicluster, can be generated for smaller energies. The emergence of this state modifies all thermodynamical and dynamical features as will be discussed in section 4.2.

The dynamics of each particle obeys the following pendulum equation of motion

$$\ddot{\theta }_i=M\mathrm{sin}(\theta _i\varphi ),$$
(20)
where $`M`$ and $`\varphi `$ have a non trivial time dependence, related to the motion of all the other particles in the system. Equation (20) has been very successfully used to describe several features of particle motion, like for instance trapping and untrapping mechanisms Antoni . There are also numerical indications Antoni  and preliminary theoretical speculations priv  that taking the mean-field limit before the infinite time limit, the time-dependence would disappear and the modulus and the phase of the magnetization become constant. This implies, as we will discuss in Section 5, that chaotic motion would disappear. The inversion of these two limits is also discussed in the contribution by Tsallis et al tsallisrap .

### 3.2 Canonical ensemble for $`A0`$

As soon as $`A>0`$, the evolution along the two spatial directions is no more decoupled and the system cannot be described in terms of a single order parameter. Throughout all this section $`ϵ=1`$. A complete description of the phase diagram of the system requires now the introduction of two distinct order parameters:

$$𝐌_z=(\frac{_i\mathrm{cos}z_i}{N},\frac{_i\mathrm{sin}z_i}{N})=M_z\mathrm{exp}(i\varphi _z)$$
(21)
where $`z_i=x_i`$ or $`y_i`$ and;

$$𝐏_{x\pm y}=(\frac{_i\mathrm{cos}(x_i\pm y_i)}{N},\frac{_i\mathrm{sin}(x_i\pm y_i)}{N})=P_z\mathrm{exp}(i\psi _z).$$
(22)
It can be shown that on average $`M_xM_yM`$ and $`P_{x+y}P_{xy}P`$: therefore, we are left with only two order parameters.

Following the approach of section 3.1, the canonical equilibrium properties can be derived analytically in the mean-field limit at . We obtain

$$\beta f=\frac{M^2+P^2}{\beta }\mathrm{ln}\left[\frac{G(M,P;A)}{\beta }\right]$$
(23)
with

$$G=_0^{2\pi }𝑑sI_0\left(M+\sqrt{2A}P\mathrm{cos}s\right)\mathrm{exp}(M\mathrm{cos}s)$$
(24)
where $`s`$ is an integration variable.

The energy per particle reads as

$$U=\frac{1}{\beta }+\frac{2+A2M^2AP^2}{2}=K+V_A.$$
(25)
Depending on the value of the coupling constant $`A`$, the single particle potential defined through $`V_A=\frac{1}{2}_iv_i`$ changes its shape, inducing the different clustering phenomena described below. For small values of $`A`$, $`v_i`$ exhibits a single minimum per cell $`(x,y)([\pi ,\pi ],[\pi ,\pi ])`$ (see Fig. 2a for $`A=1`$), while for larger values of $`A`$ four minima can coexist in a single cell (see Fig. 2b for $`A=4`$). In the former case only one clustered equilibrium state can exist at low temperatures, while in the latter case, when all the four minima have the same depth a phase with two clusters can emerge, as described in the following.

Since the system is ruled by two different order parameters $`M`$ and $`P`$, the phase diagram is more complicated than in the $`A=0`$ case and we observe two distinct clustered phases. In the very low temperature regime, the system is in the clustered phase $`CP_1`$: the particles have all the same location in a single point-like cluster and $`MP1`$. In the very large temperature range, the system is in a homogeneous phase ($`HP`$) with particles uniformly distributed, $`MP=O(1/\sqrt{N})`$. For $`A>A_23.5`$, an intermediate two-clusters phase $`CP_2`$ appears. In this phase, due to the symmetric location of the two clusters in a cell, $`MO(1/\sqrt{N})`$ while $`PO(1)`$ art . We can gain good insights on the transitions by considering the line $`T_M`$ (resp. $`T_P`$) where $`M`$ (resp. $`P`$) vanishes and the phase $`CP_1`$ (resp. $`CP_2`$) looses its stability (see Fig. 3 for more details).

We can therefore identify the following four different scenarios depending on the value of $`A`$.

(I) When $`0AA_1=2/5`$, one observes a continuous transition from the phase $`CP_1`$ to $`HP`$. The critical line is located at $`T_M=1/2`$ ($`U_M=3/2+A`$) and a canonical tricritical point, located at $`A_1=2/5`$, separates the 1st order from the second order phase transition regions333 The tricritical point has been identified by finding first the value of $`\overline{P}(M,A,\beta )`$ that minimizes $`f(M,P;A,\beta )`$ and then substituting the solution $`\overline{P}`$ in $`f`$. This reduces $`f`$ to a function of $`M`$ only. Now, the standard procedure described also for the BEG model BMRhmf  can be applied. This consists in finding the value of $`\beta `$ and $`A`$ where both the $`M^2`$ and the $`M^4`$ coefficient of the development of the free energy in powers of $`M`$ vanishes

(II) When $`A_1<A<A_23.5`$, the transition between $`CP_1`$ and $`HP`$ is first order with a finite energy jump (latent heat). Inside this range of values of $`A`$ one also finds microcanonical discontinuous transitions (temperature jumps) as we will see in Section 3.2.

(III) When $`A_2<A<A_35.7`$, the third phase begins to play a role and two successive transitions are observed: first $`CP_1`$ disappears at $`T_M`$ via a first order transition that gives rise to $`CP_2`$; then this two-clusters phase gives rise to the $`HP`$ phase via a continuous transition. The critical line associated to this transition is $`T_P=A/4`$ ($`U_P=3A/4+1`$).

(IV) When $`A>A_3`$, the transition connecting the two clustered phases, $`CP_1`$ and $`CP_2`$, becomes second order.

### 3.3 Microcanonical Ensemble

As we have anticipated, our microcanonical results have been mostly obtained via MD simulations, since we cannot estimate easily the microcanonical entropy $`S`$ analytically for the HMF model444It would be possible using large deviation technique barrethesis .. However, we show below how far we can get, starting from the knowledge of the canonical free energy, using Legendre transform or inverse Laplace transform techniques. The following derivation is indeed valid in general, it does not refer to any specific microscopic model.

The relation that links the partition function $`Z(\beta ,N)`$ to the microcanonical phase-space density at energy $`E=UN`$

$$\omega (E,N)=d^Np_id^N\theta _i\delta (EH)$$
(26)
is given by

$$Z(\beta ,N)=_0^{\mathrm{}}𝑑E\omega (E,N)\mathrm{e}^{\beta E},$$
(27)
where the lower limit of the integral ($`E=0`$) corresponds to the energy of the ground state of the model. Expression (27) can be readily rewritten as

$$Z(\beta ,N)=N_0^{\mathrm{}}dU\mathrm{exp}[N(\beta U+\frac{1}{N}\mathrm{ln}(\omega (E,N))],$$
(28)
which is evaluated by employing the saddle-point technique in the mean-field limit. Employing the definition of entropy per particle in the thermodynamic limit

$$S(U)=\underset{N\mathrm{}}{lim}\left[\frac{1}{N}\mathrm{ln}\omega (U,N)\right],$$
(29)
one can obtain the Legendre transform that relates the free energy to the entropy:

$$\beta F(\beta )=\underset{U}{\mathrm{max}}[\beta U+S(U)]\mathrm{with}\beta =\frac{S}{U}.$$
(30)
Since a direct analytical evaluation of the entropy of the HMF model in the microcanonical ensemble is not possible, we are rather interested in obtaining the entropy from the free energy. This can be done only if the entropy $`S`$ is a concave function of the energy. Then, one can invert (30), getting

$$S(U)=\underset{\beta >0}{\mathrm{min}}[\beta (UF(\beta ))]\mathrm{with}U=\frac{(\beta F)}{\beta }.$$
(31)
However, the assumption that $`S`$ is concave is not true for systems with long range interactions near a canonical first order transition, where a ”convex intruder” of $`S`$ appears BMRhmf ; grossdh , which gives rise to a negative specific heat regime in the microcanonical ensemble555A convex intruder is present also for short range interactions in finite systems, but the entropy regains its concave character in the thermodynamic limit grossdh .. For such cases, we have to rely on MD simulations or microcanonical Monte-Carlo simulations grossbook .

An alternative approach to the calculation of the entropy consists in expressing the Dirac $`\delta `$ function in Eq. (26) by a Laplace transform. One obtains

$$\omega (E,N)=\frac{1}{2i\pi }_i\mathrm{}^{+i\mathrm{}}𝑑\beta e^{\beta E}Z(\beta )$$
(32)
where one notices that $`\beta `$ is imaginary. As the partition function can be estimated for our model, we would just have to analytically continue to complex values of $`\beta `$. By performing a rotation to the real axis of the integration contour, one could then evaluate $`\omega `$ by saddle point techniques. However, this rotation requires the assumption that no singularity is present out of the real axis: this is not in general guaranteed. It has been checked numerically for the $`A=0`$ model barrethesis , allowing to obtain then an explicit expression for $`\omega `$ and confirming ensemble equivalence for this case.

However, when $`A>A_1`$, a first order canonical transition occurs, which implies the presence of a convex intruder and makes the evaluation of $`S(U)`$ through Eq. 31 impossible . For these cases, a canonical description is unable to capture all the features associated with the phase transition. For the microcanonical entropy, we have to rely here on MD simulations.

The MD simulations have been performed adopting extremely accurate symplectic integration schemes mac , with relative energy conservations during the runs of order $`10^6`$. It is important to mention that the CPU time required by our integration schemes, due to the mean-field nature of the model, increases linearly with the number of particles.

Whenever we observe canonically continuous transitions (i.e. for $`AA_1`$), the MD results coincide with those obtained analytically in the canonical ensemble, as shown in Fig. 1 for $`A=0`$. The curve $`T(U)`$ is thus well reproduced from MD data, apart from finite $`N`$ effects. It has been however observed that starting from ”water-bags” initial conditions metastable states can occur in the proximity of the transition tsallisrap .

When $`A>A_1`$, discrepancies between the results obtained in the two ensembles are observable in Fig. 4. For example, MD results in the case $`A=1`$, reported in Fig. 4(a), differ clearly from canonical ones around the transition, exhibiting a regime characterized by a negative specific heat. This feature is common to many models with long-range thir ; compa  or power-law decaying interactions tamahmf  as well as to finite systems with short-range forces grossbook . However, only recently a characterization of all possible microcanonical transitions associated to canonically first order ones has been initiated BMRhmf ; art .

For $`A`$ slightly above $`A_1`$, the transition is microcanonically continuous, i.e. there is no discontinuity in the $`TU`$ relation (this regime presumably extends up to $`A1.2`$). Before the transition, one observes a negative specific heat regime (see Fig. 4a). In addition, as already observed for the Blume-Emery-Griffiths model BMRhmf , microcanonically discontinuous transitions can be observed in the ”convex intruder” region. This means that, at the transition energy, temperature jumps exist in the thermodynamic limit. A complete physical understanding of this phenomenon, which has also been found in gravitational systems chavanishmf , has not been reached. For $`A>1.2`$, i.e. above the ”microcanonical tricritical point”, our model displays temperature jumps. This situations is shown in Fig. 5 for $`A=2`$.

For $`A>A_2`$, we have again a continuous transition connecting the two clustered phases $`CP_1`$ and $`CP_2`$. This is the first angular point in the $`TU`$ relation at $`U3.65`$ in Fig. 4(b). The second angular point at $`U4`$ is the continuous transition, connecting $`CP_2`$ to $`HP`$. The transition at lower energy associated to the vanishing of $`M`$ is continuous in the microcanonical ensemble with a negative specific heat, while discontinuous in the canonical ensemble; the dash-dotted line indicates the transition temperature in the canonical ensemble, derived using the Maxwell construction). The second transition associated to the vanishing of $`P`$ is continuous in both ensembles.

## 4 Dynamical Properties I: Out-of-equilibrium states

### 4.1 Metastable states

Around the critical energy, relaxation to equilibrium depends in a very sensitive way on the initial conditions adopted. When one starts with out-of-equilibrium initial conditions in the ferromagnetic case, one finds quasi-stationary (i.e. long lived) nonequilibrium states. An example is represented by the so-called “water bag” initial condition: all the particles are clustered in a single point and the momenta are distributed according to a flat distribution of finite width centered around zero. These states have a lifetime which increases with the number of particles $`N`$, and are therefore stationary in the continuum limit. In correspondence of these metastable states, anomalous diffusion and Lévy walks levyhmf , long living correlations in $`\mu `$space lat0  and zero Lyapunov exponents lat0a  have been found. In addition, these states are far from the equilibrium caloric curve around the critical energy, showing a region of negative specific heat and a continuation of the high temperature phase (linear $`T`$ vs $`U`$ relation) into the low temperature one. It is very intriguing that these out-of-equilibrium quasi-stationary states indicate a caloric curve very similar to the one found in the region where one gets a canonical first order phase transitions, but a continuous microcanonical one, as discussed in section 3.3. In the latter case, however, the corresponding states are stationary also at finite $`N`$. The coexistence of different states in the continuum limit near the critical region is a purely microcanonical effect, and arises after the inversion of the $`t\mathrm{}`$ limit with the $`N\mathrm{}`$ one lat0 ; lat0a ; tsallisrap .

Similarly, the antiferromagnetic HMF ($`A=0`$, $`\epsilon =1`$), where the particles interact through repulsive forces presents unexpected dynamical properties in the out-of-equilibrium thermodynamics. On the first sight, the thermodynamics of this model seems to be less interesting since no phase transition occurs as discussed above. However, thermodynamical predictions are again in some cases in complete disagreement with dynamical results leading in particular to a striking localization of energy. This aspect, as we will show below, is of course, again, closely related to the long-range character of the interaction and to the fact that such a dynamics is chaotic and self-consistent. We mean by this that all particles give a contribution to the field acting on each of them. One calls this phenomenon, self-consistent chaos Diego . In addition to the toy model that we consider here, we do think that similar emergence of structures, but even more importantly, similar dynamical stabilization of out-equilibrium states could be encountered in other long-range systems, as we briefly describe at the end of this section.

### 4.2 The dynamical emergence of the bicluster in the antiferromagnetic case

In the antiferromagnetic case, Eq. (7) with $`\epsilon =1`$, the intriguing properties appear in the region of very small energies. To be more specific, if an initial state with particles evenly distributed on the circle (i.e. close to the ground state predicted by microcanonical or canonical thermodynamics) and with vanishingly small momenta is prepared, this initial condition can lead to the formation of unstable states. This process, discovered by chance, is now characterized in full detail Antoni ; barreepjB ; drh ; firpoleyvraz .

As shown by Fig. 6, the density of particles is initially homogeneous. However a localization of particles do appear at a given time, in two different points, symmetrically located with respect to the center of the circle. This localized state, that we call bicluster, is however unstable (as shown again by Fig. 6), since both clusters are giving rise to two smaller localized groups of particles: this is the reason for the appearance of the first chevron. However, also this state is unstable, so that the first chevron disappears to give rise again to a localization of energy in two points. This state enhances the formation of a chevron with a smaller width and this phenomenon repeats until the width of the chevron is so small that one does not distinguish anymore its destabilization. Asymptotically one gets a density distribution displaying two sharp peaks located at distance $`\pi `$ on the circle, a dynamically stable bicluster.

As we have shown in Ref. barreepjB , the emergence of the bicluster is the signature of shock waves present in the associated hydrodynamical equations. Indeed, we found a strikingly good description of the dynamics of the particles by a non linear analysis of the associated Vlasov equation, which is mathematically justified BraunHepp  in the infinite $`N`$ limit. The physical explanation of this problem can be summarized as follows. Once the Hamiltonian has been mapped to the Vlasov equation, it is possible to introduce a density $`\rho (\theta ,t)`$ and a velocity field $`v(\theta ,t)`$. Neglecting the dispersion in momentum and relying on usual non linear dynamics hierarchy of time-scales, one ends up with dynamical equations at different orders in a multiscale analysis. The first order corresponds to the linear dynamics and defines the plasma frequency of order one. However, a second timescale appears that is related to the previous one by the relationship $`\tau =\sqrt{U}t`$, where $`U`$ is the energy per particle. When one considers initial conditions with a very small energy density, the two time scales are very different and clearly distinguishable by considering particle trajectories: a typical trajectory corresponds indeed to a very fast motion with a very small amplitude, superimposed to a slow motion with a large amplitude.

This suggests to average over the fast oscillations and leads to the spatially forced Burgers equation

$`{\displaystyle \frac{u}{\tau }}+u{\displaystyle \frac{u}{\theta }}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{sin}2\theta ,`$ (33)
once the average velocity $`u(\theta ,\tau )=v(\theta ,t,\tau )_t`$ is introduced. Due to the absence of dissipative or diffusive terms, Equation (33) supports shock waves and these can be related to the emergence of the bicluster. By applying the methods of characteristics to solve Equation (33), one obtains

$`{\displaystyle \frac{d^2\theta }{d\tau ^2}}+{\displaystyle \frac{1}{2}}\mathrm{sin}2\theta `$ $`=`$ $`0.`$ (34)
which is a pendulum-like equation.

Fig. 7 shows the trajectories, derived from Eq. (34), for particles that are initially evenly distributed on the circle. One clearly sees that two shock waves appear and lead to an increase of the number of particles around two particular sites, which depend on the initial conditions: this two sites correspond to the nucleation sites of the bicluster. Because of the absence of a diffusive term, the shock wave starts a spiral motion that explains the destabilization of the first bicluster and also the existence of the two arms per cluster, i.e. the chevron.

This dynamical analysis allows an even more precise description, since the methods of characteristics show that the trajectories correspond to the motion of particles in the double well periodic potential $`V(\theta ,\tau )=(\mathrm{cos}2\theta )/4`$. This potential is of mean field origin, since it represents the effect of all interacting particles. Therefore, the particles will have an oscillatory motion in one of the two wells. One understands thus that particles starting close to the minimum will collapse at the same time, whereas a particle starting farther will have a larger oscillation period. This is what is shown in Fig. 8 where trajectories are presented for different starting positions. One sees that the period of recurrence of the chevrons corresponds to half the period of oscillations of the particles close to the minimum; this fact is a direct consequence of the isochronism of the approximate harmonic potential close to the minimum. The description could even go one step further by computing the caustics, corresponding to the envelops of the characteristics. They are shown in Fig. 8 and testify the striking agreement between this description and the real trajectories: the chevrons of Fig. 6 correspond to the caustics.

### 4.3 Thermodynamical predictions versus dynamical stabilization

If the above description is shown to be particularly accurate, it does not explain why this state is thermodynamically preferred over others. Indeed, as shown before, thermodynamics predicts that the only equilibrium state is homogeneous. This result has been discussed in Section 3.1, where we have proved that magnetization is zero at all energies for the antiferromagnetic model and this has been also confirmed by Monte Carlo simulations drh . However, since MD simulations are performed at constant energy, it is important to derive analytically the most probable state in the microcanonical ensemble. Since this model does not present ensemble inequivalence, we can obtain the microcanonical results by employing the inverse Laplace transform (32) of the canonical partition function. The microcanonical solution confirms that the maximal entropy state is homogeneous on the circle. It is therefore essential to see why the bicluster state, predicted to be thermodynamically unstable is instead dynamically stable.

The underlying reason rely on the existence of the two very different timescales and the idea is again to average over the very fast one. Instead of using the classical asymptotic expansion on the equation of motions, it is much more appropriate to develop an adiabatic approximation which leads to an effective Hamiltonian that describes very well the long time dynamics. Doing statistical mechanics of this averaged problem, one predicts the presence of the bicluter.

The theory that we have developed relies on an application of adiabatic theory, which in the case of the HMF model is rather elaborate and needs lengthy calculations barreepjB  that we will not present here. An alternative, but less powerful, method to derive similar results has been used in Ref. firpoleyvraz . On the contrary, we would like to present a qualitative explanations of this phenomenon, using a nice (but even too simple !) analogy.

This stabilization of unstable states can be described using the analogy with the inverted pendulum, where the vertical unstable equilibrium position can be made stable by the application of a small oscillating force. One considers a rigid rod free to rotate in a vertical plane and whose point-of-support is vibrated vertically as shown by Fig. 9a. If the support oscillates vertically above a certain frequency, one discovers the remarkable property that the vertical position with the center of mass above its support point is stable (Fig. 9b). This problem, discussed initially by Kapitza Kapitaza , has strong similarities with the present problem and allows a very simplified presentation of the averaging technique we have used.

The equation of motion of the vibrating pendulum is

$$\frac{d^2\theta }{dt^2}+\left(\frac{\omega _0^2}{\omega ^2}a\mathrm{cos}t\right)\mathrm{sin}\theta =0,$$
(35)
where $`\omega _0`$ is the proper linear frequency of the pendulum, $`\omega `$ the driving frequency of the support and $`a`$ the amplitude of excitation. Introducing a small parameter $`\epsilon =\omega _0/\omega `$, one sees that (35) derives from the Lagrangian

$$=\frac{1}{2}\left(\frac{d\theta }{dt}\right)^2+\left(\epsilon ^2a\mathrm{cos}t\right)\mathrm{cos}\theta .$$
(36)
Here the two frequencies $`\omega `$ and $`\omega _0`$ define two different time scales, in close analogy with the HMF model. Using the small parameter to renormalize the amplitude of the excitation as $`A=a/\epsilon `$, and choosing the ansatz $`\theta =\theta _0(\tau )+\epsilon \theta _1(t,\epsilon )`$ where $`\tau =\epsilon t`$, the Lagrangian equations for the function $`\theta _1`$, leads to the solution $`\theta _1=A\mathrm{cos}t\mathrm{sin}\theta _0(\tau )`$. This result not only simplifies the above ansatz, but more importantly suggest to average the Lagrangian on the fast variable $`t`$ to obtain an effective Lagrangian $`_{eff}=_t=\frac{1}{2}\left(\frac{d\theta _0}{d\tau }\right)^2+V_{eff}`$ , where the averaged potential is found to be

$$V_{eff}=\mathrm{cos}\theta _0\frac{A^2}{8}\mathrm{cos}2\theta _0+\text{Cste}.$$
(37)
It is now straightforward to show that the inverted position would be stable if $`A^2>2`$, i.e. if $`a\omega >\sqrt{2}\omega _0`$. As the excitation amplitude $`a`$ is usually small, this condition emphasizes that the two time scales should be clearly different, for the inverted position to be stable.

The procedure for the HMF model is analogous, but of course it implies a series of tedious calculations. Since we would like to limit here to a pedagogical presentation, we will skip such details that can be found in Ref. barreepjB . It is however important to emphasize that the potential energy in the HMF model is self-consistently determined and depend on the position of all particles. The magic and the beauty is that, even if this is the potential energy of $`N`$ particles, it is possible to compute the statistical mechanics of the new effective Hamiltonian, derived directly from the effective Lagrangian via the Legendre transform. The main result is that the out-of-equilibrium state, (i.e. the bicluster shown in Fig. 6) corresponds to a statistical equilibrium of the effective mean-field dynamics. No external drive is present in this case, as for the inverted pendulum, but the time dependence of the mean field plays the role of the external drive.

The HMF model represents presumably the simplest $`N`$-body system where out-of-equilibrium dynamically stabilized states can be observed and explained in detail. However, we believe that several systems with long range interactions should exhibit behaviours similar to the ones we have observed here. Moreover, this model represents a paradigmatic example for other systems exhibiting nonlinear interactions of rapid oscillations and a slower global motion. One of this is the piston problem piston : averaging techniques have been applied to the fast motion of gas particles in a piston which itself has a slow motion sinairusse . Examples can also be found in applied physics as for instance wave-particles interaction in plasma physics plasma , or the interaction of fast inertia gravity waves with the vortical motion for the rotating Shallow Water model embidmajda .

## 5 Dynamical Properties II: Lyapunov exponents

In this section, we discuss the chaotic features of the microscopic dynamics of the HMF model. We mainly concentrate on the $`A=0`$ case, presenting in details the behavior of the Lyapunov exponents and the Kolmogorov-Sinai entropy both for ferromagnetic and antiferromagnetic interactions. We also briefly discuss some peculiar mechanisms of chaos in the $`A0`$ case. The original motivation for the study of the chaotic properties of the HMF was to investigate the relation between phase transitions, which are macroscopic phenomena, and microscopic dynamics (see Ref. laporeport  for a review) with the purpose of finding dynamical signatures of phase transitions cmd . Moreover, we wanted to check the scaling properties with the number of particles of the Lyapunov spectrum livipolruffo  in the presence of long-range interactions.

### 5.1 The $`A=0`$ case

In the $`A=0`$ case (see formulae (7)-(20)), the Hamiltonian equations of motion are

$`\dot{\theta }_i`$ $`=`$ $`p_i`$ (38)
$`\dot{p}_i`$ $`=`$ $`M\mathrm{sin}\left(\theta _i\varphi \right).`$ (39)
The Largest Lyapunov Exponent (LLE) is defined as the limit

$$\lambda _1=\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}\frac{d(t)}{d(0)}$$
(40)
where $`d(t)=\sqrt{_{i=1}^N(\delta \theta _i)^2+(\delta p_i)^2}`$ is the Euclidean norm of the infinitesimal disturbance at time $`t`$. Therefore, in order to obtain the time evolution of $`d(t)`$, one has to integrate also the linearized equations of motion along the reference orbit

$$\frac{d}{dt}\delta \theta _i=\delta p_i,\frac{d}{dt}\delta p_i=\underset{j}{}\frac{^2V}{\theta _i\theta _j}\delta q_j,$$
(41)
where the diagonal and off-diagonal terms of the Hessian are

$`{\displaystyle \frac{^2V}{\theta _i^2}}`$ $`=`$ $`M\mathrm{cos}(\theta _i\mathrm{\Phi }){\displaystyle \frac{1}{N}}`$ (42)
$`{\displaystyle \frac{^2V}{\theta _i\theta _j}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{cos}(\theta _i\theta _j),ij.`$ (43)
To calculate the largest Lyapunov exponent we have used the standard method by Benettin et al ben . In Fig. (10), we report the results obtained for four different sizes of the system (ranging from $`N=100`$ to $`N=20,000`$).

In panel (a), we plot the largest Lyapunov exponent as a function of U. As expected, $`\lambda _1`$ vanishes in the limit of very small and very large energies, where the system is quasi-integrable. Indeed, the Hamiltonian reduces to weakly coupled harmonic oscillators in the former case or to free rotators in the latter. For $`U<0.2`$, $`\lambda _1`$ is small and has no $`N`$-dependence. Then it changes abruptly and a region of “strong chaos” begins. It was observed Antoni  that between $`U=0.2`$ and $`U=0.3`$, a different dynamical regime sets in and particles start to evaporate from the main cluster, in analogy with what was reported in other models cmd ; nayak . In the region of strong chaoticity, we observe a pronounced peak already for $`N=100`$ yama . The peak persists and becomes broader for $`N=20,000`$. The location of the peak is slightly below the critical energy and depends weakly on $`N`$.

In panel (b), we report the standard deviation of the kinetic energy per particle $`\mathrm{\Sigma }`$ computed from

$$\mathrm{\Sigma }=\frac{\sigma _K}{\sqrt{N}}=\frac{\sqrt{K^2K^2}}{\sqrt{N}},$$
(44)
where $``$ indicates the time average. The theoretical prediction for $`\mathrm{\Sigma }`$, which is also reported in Fig.10, is firpo ; lat2hmf

$$\mathrm{\Sigma }=\frac{T}{\sqrt{2}}\sqrt{1\left[12M\left(\frac{dM}{dT}\right)\right]^1},$$
(45)
where $`M(T)`$ is computed in the canonical ensemble. Finite size effects are also present for the kinetic energy fluctuations, especially for $`U>U_c`$, but in general there is a good agreement with the theoretical formula, although the experimental points in Fig. 10b lies systematically below it. The figure emphasizes that the behavior of the Lyapunov exponent is strikingly correlated with $`\mathrm{\Sigma }`$: in correspondence to the peak in the LLE, we observe also a sharp maximum of the kinetic energy fluctuations. The relation between the chaotic properties and the thermodynamics of the system, namely the presence of a critical point, can be made even more quantitative. An analytical formula, relating (in the $`A=0`$ model) the LLE to the second order phase transition undergone by the system, has been obtained firpo  by means of the geometrical approach developed in Refs. lapo1 ; laporeport . Using a reformulation of Hamiltonian dynamics in the language of Riemannian geometry, they have found a general analytical expression for the LLE of a Hamiltonian many-body system in terms of two quantities: the average $`\mathrm{\Omega }_0`$ and the variance $`\sigma _\mathrm{\Omega }`$ of the Ricci curvature $`\kappa _R=\mathrm{\Delta }V/N=1/N_{i=1}^N\frac{^2V}{q_i^2}`$, where $`V`$ is the potential energy and $`q_i`$ are the coordinates of the system. Since in the particular case of the HMF model, we have

$$\frac{1}{N}\underset{i=1}{\overset{N}{}}\frac{^2V}{\theta _i^2}=M^2\frac{1}{N}=\frac{2K}{N}+12U\frac{1}{N},$$
(46)
the two quantities $`\mathrm{\Omega }_0`$ and $`\sigma _\mathrm{\Omega }`$ can be expressed in terms of average values and fluctuations either of $`M^2`$ or of the kinetic energy $`K`$

$`\mathrm{\Omega }_0`$ $`=`$ $`M^2{\displaystyle \frac{1}{N}}={\displaystyle \frac{2}{N}}K+(12U){\displaystyle \frac{1}{N}}`$ (47)
$`\sigma _\mathrm{\Omega }^2`$ $`=`$ $`N\sigma _{M^2}^2={\displaystyle \frac{4}{N}}\sigma _K^2.`$ (48)
The formula obtained by Firpo firpo  relates the LLE, a characteristic dynamical quantity, to thermodynamical quantities like $`M^2`$ and $`\sigma _{M^2}`$, or $`K`$ and $`\sigma _K`$, which characterize the macroscopic phase transition. For moderately small values of $`U`$, an approximation of the formula gives

$$\lambda \frac{\sigma _K}{\sqrt{N}}=\mathrm{\Sigma }.$$
(49)
This is in agreement with the proportionality between LLE and fluctuations of the kinetic energy found numerically in Fig. (10). This implies also a connection between the LLE and the specific heat, another quantity which is directly related to the kinetic energy fluctuation. In fact the specific heat can be obtained from $`\mathrm{\Sigma }`$ by means of the Lebowitz-Percus-Verlet formula lpv

$$C_V=\frac{1}{2}\left[12\left(\frac{\mathrm{\Sigma }}{T}\right)^2\right]^1.$$
(50)
In Fig. 11, we report the numerical results for the specific heat as a function of $`U`$ for a system made of $`N=500`$ particles, and we compare them with the theoretical estimate.

In the HMF model and for a rather moderate size of the systems, it is possible to calculate not only the LLE but all the Lyapunov exponents, and from them the Kolmogorov-Sinai entropy. We give first a succinct definition of the spectrum of Lyapunov exponents (for more details see Ruelleeckman ). Once the 2N-dimensional tangent vector $`𝐳=(\delta \theta _1,\mathrm{},\delta \theta _N,\delta p_1,\mathrm{},\delta p_N)`$ is defined, with its dynamics given by Eqs. (41), one can formally integrate the motion in tangent space up to time $`t`$, since the equations are linear,

$$𝐳(t)=J^t𝐳(0),$$
(51)
where $`J^t`$ is a $`2N\times 2N`$ matrix that depends on time through the orbit $`\theta _i(t),p_i(t)`$. The first $`k`$ exponents of the spectrum $`\lambda _1,\mathrm{},\lambda _k`$, which are ordered from the maximal to the minimal, are then given by

$$(\lambda _1+\mathrm{}+\lambda _k)=\underset{t\mathrm{}}{lim}\frac{1}{2t}\mathrm{ln}\text{Tr}J_k^t(J_k^t)^{},$$
(52)
where $`J_k^t`$ is the matrix ($`(J_k^t)^{}`$ its transpose) induced by $`J^t`$ that acts on the exterior product of $`k`$ vectors in the tangent space $`𝐳_1\mathrm{}𝐳_k`$. The spectrum extends up to $`k=2N`$ and in our Hamiltonian system obeys the pairing rule

$$\lambda _i=\lambda _{2N+1i}\text{for}1i2N.$$
(53)
The numerical evaluation of the spectrum of the Lyapunov exponents is a heavy computational task, in particular for the necessity to perform Gram-Schmidt orthonormalizations of the Lyapunov eigenvectors in order to maintain them mutually orthogonal during the time evolution. We have been able to compute the complete Lyapunov spectrum for system sizes up to $`N=100`$. In Fig. 12, we report the positive part of the spectrum for different system sizes and an energy $`U=0.1`$ inside the weakly chaotic region. The negative part of the spectrum is symmetric due to the pairing rule (53). The limit distribution $`\lambda (x)`$, suggested for short range interactions,

$$\lambda (x)=\underset{N\mathrm{}}{lim}\lambda _{xN}(N),$$
(54)
that is obtained by plotting $`\lambda _i`$ vs. $`i/N`$ and letting $`N`$ going to infinity, is found also here for the $`N`$ values that we have been able to explore. At higher energies, this scaling is not valid and a size-dependence is present lat2hmf .

The Kolmogorov-Sinai (K-S) entropy is, according to Pesin’s formula Ruelleeckman , the sum of the positive Lyapunov exponents. In Fig. 13, we plot the entropy density $`S_{KS}/N`$ as a function of $`U`$ for different systems sizes.

As for the LLE, $`S_{KS}/N`$ shows a peak near the critical energy, a fast convergence to a limiting value as $`N`$ increases in the small energy limit, and a slow convergence to zero for $`UU_c`$. A comparison of the ferromagnetic and antiferromagnetic cases is reported in Fig. 14.

Here, for N=100, we plot as a function of $`U`$ the LLE and the Kolmogorov-Sinai entropy per particle $`S_{KS}/N`$. In both the ferromagnetic and antiferromagnetic cases, the system is integrable in the limits of small and large energies. The main difference between the ferromagnetic and the antiferromagnetic model appears at intermediate energies. In fact, although both cases are chaotic (LLE and $`S_{KS}/N`$ are positive), in the ferromagnetic system one observes a well defined peak just below the critical energy, because the dynamics feels the presence of the phase transition. On the other hand, a smoother curve is observed in the antiferromagnetic case.

In the low energy regime, it is possible to work out lat1hmf  a simple estimate $`\lambda _1\sqrt{U}`$, which is fully confirmed for the ferromagnetic case in Fig. 15 for different system sizes in the range $`N=100,20000`$. The same scaling law is also valid in the antiferromagnetic case latoraprogtherophys  and for the $`A=1`$ at  HMF model.

At variance with the $`N`$-independent behavior observed at small energy, strong finite size effects are present above the critical energy in the ferromagnetic case and for all energies for the antiferromagnetic case. In Fig. 16(a), we show that the LLE is positive and $`N`$-independent below the transition (see the values $`U=0.4,0.5`$), while it goes to zero with $`N`$ above. We also report in the same figure a calculation of the LLE using a random distribution of particle positions $`\theta _i`$ on the circle in Eq. (41) for the tangent vector. The agreement between the deterministic estimate and this random matrix calculation is very good. The LLE scales as $`N^{\frac{1}{3}}`$, as indicated by the fit reported in the figure. This agreement can be explained by means of an analytical result obtained for the LLE of product of random matrices pari . If the elements of the symplectic random matrix have zero mean, the LLE scales with the power $`2/3`$ of the perturbation. In our case, the latter condition is satisfied and the perturbation is the magnetization $`M`$. Since $`M`$ scales as $`N^{\frac{1}{2}}`$, we get the right scaling of $`\lambda _1`$ with $`N`$. This proves that the system is integrable for $`UU_c`$ as $`N\mathrm{}`$. This result is also confirmed by the analytical calculations of Ref. firpo  and, more recently, of Ref. Vallejos .

In the antiferromagnetic case, the LLE goes to zero with system size as $`N^{\frac{1}{3}}`$ for all values of $`U`$.

Interesting scaling laws have also been found for the Kolmogorov-Sinai entropy in the ferromagnetic case: at small energies $`S_{KS}/NU^{3/4}`$ with no size dependence, and $`S_{KS}/NN^{1/5}`$ for overcritical energy densities. The latter behavior has been found also in other models tsuchia . Concluding this section we would like to stress that the finite value of chaotic measures close to the critical point is strongly related to kinetic energy fluctuations and can be considered as a microscopic dynamical indication of the macroscopic equilibrium phase transition. This connection has been found also in other models and seems to be quite general laporeport ; cmd ; lapo3 ; barrehmf

The behavior of the HMF model as a function of the range of the interaction celiahmf ; tamahmf ; campa  and the dynamical features before equilibration lat0 ; lat0a  is discussed in a separate chapter of this volume in connection to Tsallis nonextensive thermodynamics tsallisrap .

### 5.2 Mechanisms of chaos in the $`A0`$ case

For $`A=0`$, the origin of chaos is related to the non time-dependence of Eq. (20), since it is obvious that if the phase $`\varphi `$ and the magnetization $`M`$ would become constant the dynamics of the system will reduce to that of an integrable system. There are indeed preliminary indications priv  that in the mean-field limit $`N\mathrm{}`$, $`M`$ and $`\varphi `$ will become constant and $`\lambda 0`$ . It should be noticed that this is true if the mean-field limit is taken before the limit $`t\mathrm{}`$ in the definition of the maximal Lyapunov exponent, and numerical indications were reported in Ref. lat0a . When $`A>0`$ we expect a quite different situation: indeed, even assuming that in the mean-field limit $`M`$ and $`P`$ and their respective phases will become constant, the dynamics will eventually take place in a $`4`$-dimensional phase space and chaos can in principle be observed.

As already shown in at , for $`A=1`$ and $`\epsilon =1`$, two different mechanisms of chaos are present in the system for $`U<U_c`$ : one acting on the particles trapped in the potential and another one felt by the particles moving in proximity of the separatrix. This second mechanism is well known and is related to the presence of a chaotic layer situated around the separatrix. The origin of the first mechanism is less clear, but presumably related to the erratic motion of the minimum of the potential well, i.e. to the time-dependent character of the equations ruling the dynamics of the single particle. Indications in this direction can be found by performing the following numerical experiment. Let us prepare a system with $`N=200`$ and $`U=0.87`$ (the critical energy is in this case $`U_c2`$) with a Maxwellian velocity distribution and with all particles in a single cluster.

For an integration time $`t<2\times 10^6`$, the Lyapunov exponent has a value $`\lambda 0.13`$. But when at time $`t2\times 10^6`$, one particle escapes from the cluster, its value almost doubles (see Fig. 17). The escaping of the particle from the cluster is associated to a decrease of the magnetization $`M`$ and of the kinetic energy $`K`$. This last effect is related to the negative specific heat regime: the potential energy $`V_A`$ is minimal when all the particles are trapped, if one escapes then $`V_A`$ increases and due to the energy conservation $`K`$ decreases. As a matter of fact, we can identify a “strong” chaos felt from the particles approaching the separatrix and a “weak” chaos associated to the orbits trapped in the potential well. We believe that the latter mechanism of chaotization should disappear (in analogy with the $`A=0`$ case) when the mean-field limit is taken before the $`t\mathrm{}`$ limit. Therefore we expect that for $`N\mathrm{}`$ the only source of chaotic behaviour should be related to the chaotic sea located around the separatrix. As already noticed in Ref. at , the degree of chaotization of a given system depends strongly on the initial condition (in particular in the mean-field limit). In the latter limit, for initial condition prepared in a clustered configuration, we expect that $`\lambda =0`$, until one particle will escape from the cluster.

## 6 Conclusions

We have discussed the dynamical properties of the Hamiltonian Mean Field model in connection to ita thermodynamics. This apparently simple class of models has revealed a very rich and interesting variety of behaviours. Inequivalence of ensembles, negative specific heat, metastable dynamical states and chaotic dynamics are only some among them. During the past years these models have been of great help in understanding the connection between dynamics and thermodynamics when long-range interactions are present. Such kind of investigation is of extreme importance for self-gravitating systems and plasmas, but also for phase transitions in finite systems, such as atomic clusters or nuclei, and for the foundation of statistical mechanics. Several progresses have been done during these years. This contribution, although not exhaustive, is an effort to summarize some of the main results achieved so far. We believe that the problems which are still not understood will be hopefully clarified in the near future within a general theoretical framework. We list three important open questions that we believe can be reasonably addresses: the exact solution of the model in the microcanonical ensemble; the full characterization of the out-of-equilibrium states close to phase transitions; the clarification of the scaling laws of the maximal Lyapunov exponent and of the Lyapunov spectrum; the study of the single-particle diffusive motion levyhmf  in the various non-equilibrium and equilibrium regimes.

## Acknowledgements

We would like to warmly thank our collaborators Mickael Antoni, Julien Barré, Freddy Bouchet, Marie-Christine Firpo, François Leyvraz and Constantino Tsallis for fruitful interactions. One of us (A.T.) would also thank Prof. Ing. P. Miraglino for giving him the opportunity to complete this paper. This work has been partially supported by the EU contract No. HPRN-CT-1999-00163 (LOCNET network), the French Ministère de la Recherche grant ACI jeune chercheur-2001 N 21-311. This work is also part of the contract COFIN00 on Chaos and localization in classical and quantum mechanics.