text:
# Characterization of the stretched exponential trap-time distributions in one-dimensional coupled map lattices

## I Introduction

The decay of certain quantities characteristic of many complex systems such as glasses Phillips (1996), spin glasses Campbell (1985), quasi-crystals Dzugutov and Phillips (1995), trapping models Donsker and Varadhan (1975); Grassberger and Procaccia (1982); Barkema et al. (2001), coupled non-linear systems Hunt et al. (2002); Roy and Amritkar (1997); Pruthi (1995), turbulence Frisch and Sornette (1997); Gamba and Kolokolov (1997) and others Weibull (1951); Laherrère and Sornette (1998); Jonscher (1999) is often described by a stretched-exponential – or Kohlrausch – functional form:

$$\varphi (t)=\mathrm{exp}\left[(t/\tau )^\beta \right],t0$$
(1)
with $`0<\beta <1`$ 111By analogy, the case $`1<\beta <2`$ is often called “stretched Gaussian” (see e.g. Refs. Kaplan and Heller, 1996; van Zon and MacKintosh, 2002).. Although Eq. (1) provides a good fit to a wide range of experimental and numerical results, in many cases these can also be fitted by power laws with comparable accuracy; while most experimental set-ups can span many decades in time, few have achieved more than 2 or 3 decades in $`\varphi (t)`$. The Kohlrausch form can also be reproduced explicitely by a few theoretical models (see, e.g., Refs. Donsker and Varadhan, 1975; Grassberger and Procaccia, 1982; Campbell, 1985; Gielis and Maes, 1995; Garrhan and Chandler, 2002) on the basis of various assumptions, but it is still unclear whether there can be a unifying solid theoretical justification for it.

Recently, Hunt, Gade and Mousseau Hunt et al. (2002) found that stretched exponentials could fit experimental distributions of trap time — the time a system spends in an uninterrupted state with temporal period two, in a one-dimensional network of coupled diode resonators — over more than 6 decades in the distribution. Numerical models, based on the set-up, could expand this fitting over about 10 orders of magnitude, strongly suggesting that dynamics with underlying stretched exponential distributions could be universal in coupled chaotic systems and ruling out a power-law or any other standard fit.

In this paper, we revisit one such numerical model and provide a detailed characterization of the dynamics of this network as a function of system size and parameters, providing elements of explanation regarding the origin of stretched exponential distributions in this system. In particular, we show that (1) size effects are important only for relatively small systems; (2) the natural invariant density, $`\rho (x)`$, generated by typical orbits has well-defined structure that is self-organized; (3) the structure of $`\rho (x)`$ alone cannot lead to stretched-exponential distributions — the dynamical spatial organization is essential for stabilizing the periodic orbits.

This paper is organized as follows: In the next section, we introduce the mathematical model: a one-dimensional coupled map lattice. Section III discusses the background for the trap-time distributions. In section IV, the results from the simulation are presented, to be discussed in section V. We summarize our conclusions in Sect. VI.

## II Model

Our model is a one-dimensional chain of $`N`$ diffusively coupled nonlinear deterministic maps, $`f(x)`$, with a coupling constant $`\alpha `$ and periodic boundary conditions. The interaction is totalistic and involves only the nearest neighbors. The time evolution of this system is discrete and is described by an iterative equation:

$$x_n(t+1)=(1\alpha )f[x_n(t)]+\frac{\alpha }{2}\left\{f[x_{n1}(t)]+f[x_{n+1}(t)]\right\}.$$
(2)
which, given initial conditions $`x_n(0)`$ for each site $`n=1,2,\mathrm{},N`$, generates a time series $`\{x_n(t)\}`$, for integer values of the time index $`t=0,1,\mathrm{}`$. This model offers a much simplified version of the 1D array of diode resonators studied in Refs. Johnson et al., 1996; Löcher et al., 1998a, b; Löcher et al., 2000. Although the individual diode resonators are best described by an inertial equation Rollins and Hunt (1984), it was shown recently that Eq. 2 captures the dynamics in the regime of interest here Hunt et al. (2002).

The diffusive nature of the coupling in Eq. (2) becomes more apparent if this equation is rewritten in the following form

$`x_n(t+1)`$ $`=`$ $`f[x_n(t)]+`$ (3)
$`{\displaystyle \frac{\alpha }{2}}\left\{f[x_{n1}(t)]2f[x_n(t)]+f[x_{n+1}(t)]\right\}`$
$`=`$ $`f[x_n(t)]+{\displaystyle \frac{\alpha }{2}}D_+D_{}f[x_n(t)]`$
where the central difference operator $`D_+D_{}f(x_n)`$ is the discrete Laplacian of $`f(x_n)`$ on a one-dimensional grid with unit spacing. In the limit $`N\mathrm{}`$ the discrete Laplacian can be substituted by its continuous counterpart and Eq. (3) becomes a nonlinear analogue of a difference-differential diffusion equation, see e.g., Ref. Elmer and Vleck, 2001. For a large $`N`$, $`1/r`$ can be thought of as a “viscosity” and $`\alpha `$ as a “diffusivity” with their usual relation to the “temperature.”

We use the logistic map $`f(x)=rx(1x)`$, Ref. Ott, 1993, in order to describe the dynamics of the basic chaotic elements. Although this differs from the form studied in Ref. Hunt et al., 2002, $`g(y)=1ay^2`$, the results are unaffected; the two maps are conjugate, related by a simple algebraic transformation: $`x=\sqrt{a}y/\sqrt{r}+1/2`$, $`g(y)=f[x(y)]r/4+1`$. We choose the logistic map because it is more studied in the literature. In the simulations presented here, we fix the value of the coupling constant to $`\alpha =0.25`$, and vary $`N`$ and $`r`$, which we refer to as the nonlinearity parameter in the following.

Initial values of $`x_n(0)`$ are taken from a random distribution in the $`[0,1[`$ interval. Runs are then iterated for a few hundred thousands steps in order to avoid any transient effect before starting the accumulation of data. All simulations presented here are done on one-dimensional arrays with variable length and periodic-boundary conditions. Statistics are generally accumulated over 10 million to 10 billion time steps.

Following the experiment Hunt et al. (2002), the analysis of the dynamics is done using coarse-grained variables defined by

$$\sigma _n(t)=\mathrm{sign}[x_n(t)x_{\mathrm{thr}}],t=0,2,4,\mathrm{}$$
(4)
where the quantity $`x_n(t)`$ is defined in Eq. (2) and $`x_{\mathrm{thr}}`$ is a certain threshold value. The results presented here are not very sensitive to the value of the threshold. For simplicity, we select the value of the unstable fixed point, $`x_{}=11/r`$, for the single map. This coarse-graining reduces the problem from continuous to two-state.

The basic dynamics of the coupled oscillator being period two, the analysis is also done only over even (or odd) time steps.

## III Background: the distribution of trap time

We are interested in the statistical distribution of trap time in a coarse-grained state space of a one-dimensional chain of coupled nonlinear maps, as shown in Fig. 1. This quantity is formally equivalent to the distribution of time intervals between zero crossings of renewal processes such as random walks Godrèche and Luck (2001) and has the advantage that it can be measured experimentally and numerically to a high degree of accuracy for this system Hunt et al. (2002).

Generally, however, experimentally measurable relaxation responses are mathematically described by auto-correlation functions of certain dynamical variables. For example, the inverse Fourier transform of the dynamical structure factor, the intermediate scattering function, can be calculated as time auto-correlation function of the microscopic density distribution of particles in a material Hansen and McDonald (1986).

The relation between the auto-correlation function $`C(t)=\sigma (t^{})\sigma (t^{}+t)_t^{}`$ of a renewal process $`\sigma (t)`$ —defined analogously to Eq. (4) — and the distribution of trap time of this process was recently studied for a range of distributions Godrèche and Luck (2001). Using this framework, it is possible to show that $`C(t)`$ corresponding to the stretched exponential distribution of occupation time is also stretched exponential at long time albeit with a different stretching exponent $`\beta `$.

These results are only valid when the traps are uncorrelated in time which is the case for the systems studied here. The details of this work will be presented elsewhere Simdyankin and Mousseau .

This establishes a direct relation between the distribution of trap time, discussed in this paper, and the more standard auto-correlation function, measured in a number of experiments on glasses and other complex systems.

## IV Results

Before discussing the possible origins for the stretched exponential distributions in our model, it is necessary to offer a characterization of its dynamics as a function of the parameters of the model.

First, we consider the size effects on the dynamics of the network, studying experimental set-up varying in length between 15 and 256 oscillators, and simulated arrays of up to tens of thousands of sites.

Then we investigate the dynamics of the system as a function of the nonlinearity parameter, $`r`$, which brings the system through a series of dynamical changes from stable periodic orbits to full chaos. In particular, it is important to assess the parts of the parameter space where the stretched exponential distributions can be observed.

Once the basic phase diagram is established, we discuss the building up of the dynamics on a single site embedded in the network as well as the spatial structure associated with the stretched exponential distribution.

### IV.1 Dependence on $`N`$

As is seen in Fig. 1, the traps are directly associated with a periodic spatial organization, which is controlled by the coupling $`\alpha `$.

As $`\alpha `$ is increased, the spatial organization goes through two rapid transitions, showing qualitatively different features: For $`\alpha <0.1`$, the lattice tends to follow a period-two spatial organization. For intermediate values of $`0.1<\alpha <0.19`$, the lattice immediately freezes into a period-two state in both space and time, for all values of the driving parameter, $`r`$. Above this threshold, $`\alpha >0.19`$, the dynamics becomes stochastic again while dominant spatial period goes to 4 and even longer for large $`\alpha `$. In each of these phases, the variation of $`\alpha `$ affects only minimally the spatial structure of the phase.

In spite of the evident organization seen in Fig. 1, the spatial correlation is short-range. Figure 2 shows that spatial correlations vanish exponentially fast with a typical length scale between about 2 and 8, i.e., the only static spatial correlation appearing in the system is directly associated with the short-range organization.

This sort-range correlation implies that size effects should be very limited. Fig. 3 shows the lattice-size dependence on the bifurcation diagram, for a single site on the lattice.

For small lattices, $`1<N<32`$, the dynamics depends on whether $`N`$ is even or odd. For even $`N`$, arrays synchronize rapidly while the odd sizes continue to display chaotic trajectories. For example, the bifurcation diagram changes qualitatively as one goes from a single isolated site to a chain of 8 or 16 oscillators: the chaotic phase disappears totally and the system remains periodic in the coarse-grained state space until $`r=4.0`$ \[see Fig. 3(b) and (c)\]. While isolated oscillators produce an exponential trap-time distribution, this distribution tends to the stretched exponential form at short times for lattices with odd $`N`$ as small as 15. Size effects are still present for these lattices, however, and the long-time distribution diverges from the stretched exponential \[see Fig. 4 (a)\]. Interestingly, the sign of the deviation for the stretched exponential oscillates as the array is increased in size two by two, indicating a certain spatial frustration in these small systems.

For even $`N`$ that is not a multiple of four, the mismatch with the periodicity results in the presence of inclusions of stable defects, with four sites at a row in the same band. These defects can also be present and stable in lattices where $`N`$ is a multiple of four.

As the number of elements reaches 32, orbits corresponding to larger values of $`r_{\mathrm{min}}r4`$ become chaotic and the bifurcation diagram approaches that of an infinite lattice. $`r_{\mathrm{min}}`$ continues to decrease with growing number of sites, and is well converged for a lattice of a few hundred oscillators. As shown in Fig. 4 (b) the same trend is seen experimentally, although the inherent disorders helps to decrease the finite size effects.

Even though the bifurcation diagrams for an isolated map and the infinite lattice appear similar, their respective orbits are qualitatively different for most values of $`r`$. In particular, most of the chaotic region of the latter shows a stretched exponential trap-time distribution.

### IV.2 Dependence on $`r`$

Fig. 3 also indicates the effect of $`r`$ on the dynamics of a single site in a chain. The coupling stabilizes the orbits, reducing significantly the size of the chaotic region. Its effect is to shift the bifurcation diagram to the right, moving the dynamics towards periodic orbits as is shown in Fig. 5 (Fig. 5 will be discussed in more detail in Sect. IV.3).

The stretched exponential behavior appears, for a large enough lattice, around $`r_0=3.83`$. At this threshold value, the stretching exponent $`\beta 0.33`$, obtained from fitting the long-time part of the trap distribution (Fig. 6) according to Eq. (1). $`0.33`$ is the lowest value of $`\beta `$ we could observe for the explored regions of the parameter space of the model. Below this threshold, chaotic windows are interspersed with periodic ones. Trajectories in these narrow windows tend to fall erratically onto neighboring periodic orbits, generating multiple-step trap distributions with no clear overall time behavior.

As seen in Fig. 6, the value of $`\beta `$ increases with $`r`$. For $`r=3.8888`$, the trap-time distribution is well fitted by a stretched exponential with $`\beta 0.5`$, while at $`r=4`$, the full chaos limit for an isolated logistic map, a fit to $`\rho (t)`$ gives $`\beta =0.70\pm 0.05`$ and $`\tau =3.0\pm 0.5`$.

Because of the short spatial correlation and the excellent quality of the simulation data, the exact value of $`\beta `$ is well defined: increasing the lattice size from 64 to 1000 sites leaves the trap-distribution essentially unchanged. Moreover, the same values of $`\beta `$ are obtained, within the error bars, by changing the length of the interval over which the fitting procedure is performed.

These numbers are also in reasonable agreement with experimental results which show $`\beta `$ varying from $`0.10\pm 0.05`$ to $`0.95\pm 0.05`$. This wider range for the experiment is probably caused by the presence of site disorder.

This confirms that the best functional form for the simulation data is given by the Kohlrausch function.

### IV.3 Single site dynamics.

The iterative rule Eq. (2) can be viewed as an equation describing the dynamics of a single element in the presence of an external additive perturbation $`a(t)`$:

$$x(t+1)=(1\alpha )f[x(t)]+\alpha a(t)$$
(5)
where $`a(t)=\{f[x_{n1}(t)]+f[x_{n+1}(t)]\}/2`$. For simplicity, since all sites are statistically identical, we drop the subscript $`n`$ in Eq. (5).

This form allows us to concentrate on the impact, at the single-site level, of the rest of the network and to try to identify the essential elements for a stretched exponential dynamics. For this purpose, it is convenient to use the natural invariant density $`\rho (x)`$ generated by typical orbits $`\{x(t)\}`$, $`t=1,2,\mathrm{}`$ of the map in Eq. (5). The notion of the natural invariant density is widely used in the studies of chaos Ott (1993). The function $`\rho (x)`$ is defined so that for any interval $`[x,x+dx][0,1]`$ the fraction of the time typical orbits spend in this interval is $`\rho (x)dx`$. In the same way one can define the density $`\rho (a)`$ for the perturbation $`a(t)`$ in Eq. (5).

Fig. 7 shows the behavior of $`\rho (a)`$ and $`\rho (x)`$ for a few sets of the parameters $`\alpha `$, $`N`$, and $`r`$. In the case of an isolated map ($`\alpha =0`$) and in the chaotic regime, $`a(t)`$ is the sum of two identically distributed independent random variables — $`\rho (a)`$ is shown by the dash-dotted line in Fig. 7 (b). According to the central limit theorem, if this sum contained much more than two term, $`\rho (a)`$ would very well agree with the normal (Gaussian) distribution. However, even with only two terms, $`\rho (a)`$ exhibits a maximum at $`a=1/2`$ \[the dash-dotted line in Fig. 7 (a)\].

The presence of one or few peaks in $`\rho (a)`$ makes the orbits $`\{x(t)\}`$ intermittent for some values of the nonlinearity parameter $`r`$. This intermittency means that the (coarse-grained) orbit $`\{\sigma (t)\}`$, see Eq. (4), stays periodic during some time, becomes chaotic for a certain while, and goes periodic again. In order to demonstrate how a periodic attractor develops in a coupled system at a value of $`r`$ for which the isolated map produces a chaotic orbit, we can imagine a density $`\rho (a)=\delta (aa_0)`$, where $`\delta (a)`$ is the delta function and $`a_0[0,1]`$ is a constant. In this case, Eq. (5) becomes $`x(t+1)=(1\alpha )f[x(t)]+\alpha a_0=c_1f[x(t)]+c_2`$, which is the iterative equation for the isolated map scaled by a constant $`c_1`$ and shifted by another constant $`c_2`$. This manipulation with the isolated map results in a bifurcation diagram translated with respect to the original one \[Fig. 3 (a)\] while preserving all its qualitative features, see Fig. 5. In particular, the period-two attractors of the modified map overlap with chaotic orbits of the original map.

As discussed above, in the interval $`0.1<\alpha <0.19`$, the lattice falls into period-two spatial and temporal patterns, i.e., each two neighboring sites oscillate out of phase and each site alternates between two values $`x^{(1)}<x^{(2)}`$, as indicated in Fig. 7 (a) and (b).

For larger values of $`\alpha `$, the natural invariant density becomes more complex. When the lattice size is a small even number, the spatial pattern is periodic with the period which is equal to four (two sites “up” followed by two sites “down”) and each site again undergoes a period-two dynamics \[Fig. 7 (c) and (d)\] which is manifested by the peaks in $`\rho (x)`$. For $`N=16`$, there are narrow peaks in $`\rho (a)`$ and $`\rho (x)`$. These peak broaden and merge rapidly, however, for large $`N`$ and the system displays strong sign of chaoticity.

This chaoticity can be controlled by $`r`$ and we see that the orbits spend more and more time in period-two attractors relative to the chaotic ones as the nonlinearity parameter is decreased, relating directly to a change in the value of the stretching exponent $`\beta `$ \[Fig. 7 (e) and (f)\].

Although providing a clearer picture of the dynamics, knowledge of the invariant density is not sufficient to reproduce the stretched exponential dynamics. For example, choosing randomly the values of $`a(t)`$ in Eq. (5) according to a prescribed $`\rho (a)`$ generates only an exponential trap-time distribution. Variations on this theme, including the introduction of a two-step distribution, which favors staying in the period-two phase once in it, also fail to give a stretched exponential distribution of trap times. Spatial correlations, even though short range, are essential to induce a stable and self-organized stretched exponential distribution.

### IV.4 Stability of spatial period-4 structure.

Fig. 1 plots the time sequence for $`\sigma _n(t)`$ on a small region of a 1000-site chain for $`r=3.83`$ ($`\beta 0.33`$) and $`r=3.8888`$ ($`\beta 0.5`$). These show a remarkable range in the trap size: the larger traps of the first one can be up to 250 times longer than those at $`r=3.8888`$, and corresponding to a time-scale 30,000 times larger than the basic time-step. This time scale is similar to that observed experimentally Hunt et al. (2002). In spite of these very stable traps, for long enough time, all sites display an identical dynamics and no region of the lattice remains frozen.

Looking at the same figure, we also see that, for a given $`r`$, the trap width appears to be uncorrelated with its length: long traps remain narrow and are composed of ordered domains with several occurrences of period 4, 5 or 6. While the period-4 domains can remain totally stable for a long time, the middle site of the 3-up or 3-down segment in period 5 and period 6 domains show some instability. They must therefore be considered more as a defect in the stable phase than as an additional spatial structure. Even so, these larger basic domains, that occur mostly for the lower values of $`r`$, are responsible for the apparent spatial periodicity of 5 in the spatial correlation function at $`r=3.83`$ (see Fig. 2).

In view of this discussion the space-time diagram can be separated in two distinct phases: the ordered spatial period 4 and the chaotic phase. The former can only be destroyed at the boundary or by spatial defects (domains of period 5 or 6).

We verify this observation by two simulations. First, the entire chain is initially set into period-4 state with the period formed by two neighboring sites with $`x_n(0)=x^{(1)}`$, $`n=1,2`$ and the next two sites with $`x_n(0)=x^{(2)}`$, $`n=3,4`$. The values $`x^{(1)}`$ and $`x^{(2)}`$ are taken to be equal to those corresponding to the two peaks in $`\rho (x)`$ in Fig. 7 (f) for $`r=3.83`$ and $`N=1000`$. Simulations starting from this initial conditions remain frozen in this spatial period-4 and temporal period-2 state indefinitely. Although metastable, this state requires the presence of defects to be destroyed.

This can be checked in a second test. We first iterate a 1000-site lattice at $`r=3.83`$, for $`10^6`$ time-steps in order to eliminate any transient effects. At that point, at once, half of the chain sites are set into the same spatial period-4 state as above. As shown in Fig. 8, this phase then disappears gradually, over about 30,000 time steps, invaded from the edges by the chaotic phase.

## V Discussion

The previous sections have shown that the dynamics of the network can be understood in terms of a competition between two dynamical regimes: a stable period-two orbit and a fully chaotic state. This can be deduced from Figs. 5 and 7: when the neighboring sites are in opposite bands, their contribution shifts the bifurcation diagram into a period-two regime. A spatial period 4 is therefore the basic stable motif with this set of parameters.

As mentioned above, we cannot reconstruct a stretched exponential dynamics simply using the natural invariant density inserted into Eq. (5). We have found that the only way to obtain a stretched exponential trap-time distribution in this situation is to impose a stretched exponential trap-time distribution on $`a(t)`$ for the values of $`a`$ within the peak region in $`\rho (a)`$ (see Fig. 7). More interestingly, connecting a single site, with a unidirectional coupling, to two sites selected at random on a lattice, is sufficient to induce a stretched exponential trap distribution on this single site, albeit with a larger $`\beta `$ than that for the elements which belong to the lattice.

Thus once an external perturbation, following a stretched exponential dynamics, is imposed on a chaotic element, the latter will immediately adopt a similar dynamics. However, stretched exponential distributions cannot be observed in a self-organized process without spatial organization. For example, connecting the nodes of a balanced binary tree unidirectionally towards its root results in the densities $`\rho (x)`$ and $`\rho (a)`$ similar to those depicted in Fig. 7, (e) and (f) already for two levels in the tree. Nevertheless, the trap-time distribution measured at the root node is a pure exponential. This suggests that the stretched exponential distributions requires some spatial organization allowing time-limited spatial period-4 structures to occur and to be stabilized.

The nature of this spatial organization is somewhat paradoxical. Correlations are very short range and the width of these domain is only weakly varying with the stretching parameter. While traps at $`r=3.83`$ are up to 250 times longer than those at $`r=3.8888`$, their width is only 4 to 5 times larger, as displayed in Fig. 1. This decoupling between spatial and temporal correlation is reminiscent of dynamical heterogeneities observed in in glass-forming systems Garrhan and Chandler (2002); Vollmay-Lee et al. (2002); Dzugutov et al. (2001). Associating a “temperature” with the control parameter $`r`$, we see that the spatial size of the traps increases with decreasing $`r`$ without diverging. In the same way, the size of dynamical heterogeneities increases and the stretching exponent $`\beta `$ decreases under cooling in glass-forming systems.

As we showed also, once in a perfect spatial period-4 regime, the system will never become chaotic. The finite life-time arises from defects in the period-4 phase or from the chaotic boundaries. This behavior is similar to what is seen in the spin-system models with effective constraints on the dynamics Garrhan and Chandler (2002) and the state-space partition model Nadler et al. (1996). Interestingly, both models show a stretched exponential decay of certain statistical quantities. The distinctive feature of the present model, Sect. II is that the effective dynamical constraints arise in the course of the evolution described by the dynamical equation (2), unlike the above cited models where the fixed constrains are imposed on the neighboring sites and the dynamics is due to a Monte Carlo procedure.

## VI Conclusion

The coupled array of chaotic oscillators presented here has a number of properties that make it an important model. As was shown above and in Ref. Hunt et al., 2002, it is a faithful representation of the dynamics of an experimental set-up of coupled diode-resonators in chaotic regime characterized by a stretched exponential distribution of trap time. Moreover, the low cost associated with solving numerically the model of Eq. (2) allows us to study the system on time scale unreachable with atomic models; it takes less than a day, on a fast processor, to iterate a lattice of 1000 sites over $`10^9`$ steps. With traps extending to more than 30,000 times the basic time step, long simulations are absolutely necessary to establish the nature of the dynamics in these systems.

Results of these simulations demonstrate that the stretched exponential distributions arise from the competition between a chaotic and a period-two regimes. The stretched-exponential requires some spatial organization to appear but does no imply diverging length scale as the traps become longer and longer: the period-4 structure is sufficient to stabilize a site onto a periodic orbit and the space of the spatial correlation is only weakly related to the length of the traps appearing in these systems.

There are strong similarities between this system and the configurational glasses and we suggest that the understanding gained here could be extended directly to these important materials. We are currently pursuing this avenue of research.

## Acknowledgments

We thank C. Godrèche and I.A. Campbell for valuable comments. This work in supported in part by the Office of Naval Research (ERH), and by the Natural Sciences and Engineering Council of Canada and the NATEQ fund of Québec (NM and SIS). Most of the calculations were one on the computers of the Réseau québécois de calcul de haute performance (RQCHP). NM is a Cottrell Scholar of the Research Corporation.


format: text

text:
# Bianchi type I cosmology and the Euler-Calogero-Sutherland model

## 1 Introduction

Dealing with cosmological studies, the conceptual problem of the identification of true degrees of freedom of gravity, the “golden fleece” of canonical gravity , transforms into a real practical problem. In the framework of the conventional formulation of Hamiltonian dynamics of the cosmological models (see, e.g., ), based on the Dirac  and the Arnowitt-Deser-Misner (ADM) treatment of general relativity , the solution of the problem involves two ingredients: choosing a certain coordinate fixing condition (gauge) and solving the constraints. In the application to the minisuperspace models of Bianchi type, the standard choice of the intrinsic time gauge leads to a reduced Hamiltonian which is in general time-dependent and represents a square-root expression . In spite of the fact that a “relativistic” physicist takes delight in this square-root expression as a reminder of the relativistic particle energy, the obtained result is not quite satisfactory. Particularly, passing to the quantum theory, it brings the unpleasant difficulties of defining the quantized square-root Hamiltonian as an operator in Hilbert space.

In the present paper, we would like to draw attention to the fact that for the simplest Bianchi type I cosmological model, it is possible to use an alternative more attractive generator of evolution, free from the above-mentioned difficulties. Below it will be shown that one can arrive at a reduced, quadratic in momenta, time-independent Hamiltonian, in contrast to the above-mentioned textbook result. Moreover, a transparent correspondence between the resulting deparametrized theory and a well-known nonrelativistic many-body integrable system naturally arises. Namely, it turns out that the expression for the generator of evolution of the unconstrained Bianchi type I cosmology coincides with the Hamiltonian describing the relative motion of the particles in the three-body hyperbolic Euler-Calogero-Sutherland model. 111 This system is a generalization of the many-body Calogero-Sutherland model  (see also the useful reviews ), describing particles on a line and interacting via pairwise hyperbolic potential $`1/\mathrm{sinh}^2x`$, by introduction of spin variables of the particles . For the other types of generalizations by inclusion of spinlike variables in the Calogero-Sutherland model, see, e.g., .

The plan of the paper is as follows. The first part offers a short review of Bianchi cosmological models including as well geometrical aspects as its Dirac generalized dynamics . Then, we restrict our consideration to the simplest Bianchi type I cosmology. Analyzing the energy constraint and the extended Poincaré-Cartan 1-form, we show how to deparametrize the theory, identify the physical degrees of freedom, and as a result establish the correspondence with the three-body Euler-Calogero-Sutherland model. Finally, using the relation between the dynamics of the Euler-Calogero-Sutherland system and the geodesic motion on the space of symmetric matrices, we discuss the general solution of the Bianchi type I model.

## 2 Bianchi cosmology: geometry and dynamics

The canonical formulation of the cosmological models intensively uses the geometrical symmetry to restrict the gravitational configuration space, the space of all possible pseudo-Riemannian metrics defined on a given manifold. In the case of large symmetry of the space-time manifold, the gravitational degrees of freedom are effectively reduced to finite number and this circumstance essentially relieves the analysis of the theory. Below we describe very briefly this effective gravitational configuration space for the homogeneous Bianchi type models. For details, we refer to the comprehensive reviews .

### 2.1 Geometrical settings

The conventional Hamiltonian analysis requires the existence of a preferred timelike variable and therefore views the Universe in terms of space plus time. Below in quick review of the $`3+1`$ geometry of gravitation and its canonical formulation, we follow the definitions and notations from the very transparent review by Isenberg and Nester .

So, the space-time is supposed to be a smooth manifold, $`(,𝐠)`$ endowed with a metric $`𝐠`$ of signature $`(,+,+,+)`$, metric-compatible connection, and foliated by a family of one-parameter three-dimensional nonintersecting surfaces $`\mathrm{\Sigma }_t`$ on which the induced metric has three positive eigenvalues, i.e., $`=\times \mathrm{\Sigma }_t`$. According to such a space-time foliation, it is useful to choose on $``$ a surface-compatible frame of vector fields 222 Following , the boldface notation is used to distinguish four-dimensional quantities from three-dimensional ones. $`(\text{e}_{},e_1,e_2,e_3)`$ with timelike unit-length vector field $`\text{e}_{}`$ orthogonal to $`\mathrm{\Sigma }_t`$ and three spacelike vector fields 333 Hereafter, the latin indices run over 1,2,3. $`e_a=(e_1,e_2,e_3)`$tangent to it. In the corresponding frame of forms $`(𝜽^{},\theta ^a)`$, the metric g reads

$$\text{g}=𝜽^{}𝜽^{}+\gamma _{ab}\theta ^a\theta ^b,$$
(1)
with a spatial metric $`\gamma `$ induced on $`\mathrm{\Sigma }_t`$. According to this representation, the Lie derivative $`_\text{e}_{}`$, derivative with respect to the proper time along the normal to the hypersurface $`\mathrm{\Sigma }_t`$, describes the evolution of the dynamical variables. Note that the frame $`(\text{e}_{},e_a)`$ is not a coordinate frame and usually, instead of dealing with Eq. (1), the so-called ADM metric  is exploited. The latter is based on using the coordinate vector fields $`\text{e}_0=\mathbf{}\mathbf{/}\mathbf{}𝒕`$ and $`e_a=/x^a`$, where the vector field $`\text{e}_0`$ is expressed in terms of the normal vector field $`\text{e}_{}`$ and spatial vector field $`N^ae_a`$ tangent to the hypersurface $`\mathrm{\Sigma }_t`$ as

$$𝐞_0=N𝐞_{}+N^ae_a.$$
(2)
However, when investigating the homogeneous Bianchi cosmological models, it is more convenient to deal with a special noncoordinate frame instead of the ADM coordinate frame for the 3-space. Its construction is dictated by geometrical considerations. By definition, in a spatially homogeneous space-time, a three-dimensional Lie group $`G_3`$ acts on the space-time as a group of isometries, such that each orbit on which $`G_3`$ acts simply transitively is a spacelike hypersurface. The advantage of considering the simply transitive action is that we can put the element of $`G_3`$ into one-to-one correspondence with the points of $`\mathrm{\Sigma }_t`$ and thus the space-time is considered topologically as the product space $`=\times G_3`$. Based on this observation, it is clear that instead of the usual coordinate frame of the spatial vector fields we need to choose a new space basis $`e_a`$, adapted to the Lie group structure of the three-dimensional hypersurface $`\mathrm{\Sigma }_t`$. Namely, the algebra of the infinitesimal generators of isometries, i.e., the Killing vector fields $`\xi _a`$,

$$[\xi _a,\xi _b]=C_{ab}^c\xi _c,$$
(3)
dictates what kind of basis should be chosen. Indeed, the vector fields $`\xi _a`$ provide a basis $`(\text{e}_0,e_a)`$ invariant under the isometries

$$_{\xi _a}\text{e}_0=0,_{\xi _a}e_a=0.$$
(4)
In this case, one can specify the form of the three-dimensional part of the Bianchi metric as

$$𝜸=\gamma _{ab}\omega ^a\omega ^b,$$
(5)
where $`\omega ^a`$ are group-invariant 1-forms dual to the vector fields $`e_a`$ with structure coefficients $`C_{bc}^a=2d\omega ^a(e_b,e_c)`$, defined by the structure constants of the homogeneity group $`C_{bc}^a`$.

Thus, instead of the ADM coordinate frame, the adapted representation for the Bianchi-type metrics in noncoordinate basis is

$$\text{g}=(N^2N^aN_a)\text{dt}\text{dt}+2N_a\text{dt}\omega ^a+\gamma _{ab}\omega ^a\omega ^b.$$
(6)
The preferable role of this choice for the frame is transparent: from the Killing equation

$$_{\xi _a}\text{g}=0$$
(7)
it follows that the functions $`N,N^a`$, and $`\gamma _{ab}`$, entering in the expression for the Bianchi metric (6), depend only on the time parameter $`t`$.

### 2.2 Generalized dynamics

As a result of this geometrical analysis, in an adapted basis, all gravitational field configuration variables, the lapse function $`N`$, the shift vector $`N^a`$, and the spatial metric $`\gamma _{ab}`$ can be identified with a set of ten Lagrangian coordinates of a certain “mechanical system.”

By definition, the dynamics of this system is determined from the Hilbert action

$$𝒜=_{}𝝈\text{R},$$
(8)
where the ansatz (6) is plugged into the expression for the space-time scalar curvature R and for the four-dimensional volume element $`𝝈=\sqrt{𝒈}\omega ^0\omega ^1\omega ^2\omega ^3.`$ As a result the restricted variational problem (for Bianchi class A models) 444 Bianchi class A models are defined writing the structure constants of the isometry Lie group as $`C_{ab}^d=ϵ_{lab}S^{ld}+A_{[d}\delta _{b]}^d`$ and supposing that $`C_{ad}^d=A_a=0`$. Here $`S^{ab}`$ is a three-dimensional second-rank symmetric tensor and $`A_a`$ is a three-dimensional vector.   is

$$𝒜[N,N_a,\gamma _{ab},\dot{\gamma }_{ab}]=\underset{t_1}{\overset{t_2}{}}𝑑t\sqrt{\gamma }N\left({}_{}{}^{3}RK_a^aK_b^b+K_{ab}K^{ab}\right)+2\underset{t_1}{\overset{t_2}{}}𝑑t\frac{d}{dt}\left(\sqrt{\gamma }K_a^a\right),$$
(9)
where $`{}_{}{}^{3}R`$ is the scalar curvature formed from the spatial metric $`\gamma `$,

$${}_{}{}^{3}R=\frac{1}{2}\gamma ^{ab}C_{da}^cC_{cb}^d\frac{1}{4}\gamma ^{ab}\gamma ^{cd}\gamma _{ij}C_{ac}^iC_{bd}^j,$$
(10)
and

$$K_{ab}=\frac{1}{2N}\left[(\gamma _{ad}C_{bc}^d+\gamma _{bd}C_{ac}^d)N^c+\dot{\gamma }_{ab}\right]$$
(11)
is the extrinsic curvature of the slice $`\mathrm{\Sigma }_t`$ defined as $`K_{ab}=\frac{1}{2}_e_{}\gamma _{ab}`$.

The Lagrangian (9) belongs to the class of so-called degenerate ones; its Hessian is zero. To deal with its Hamiltonian description, we shall follow the Dirac generalization of Hamiltonian dynamics .

Implementing the Legendre transformation on variables $`N,N^a`$, and $`\gamma _{ab}`$, we get the canonical Hamiltonian

$$H_C=N+N^a_a$$
(12)
the primary

$$P=0,P_a=0$$
(13)
and secondary constraints

$`={\displaystyle \frac{1}{\sqrt{\gamma }}}\left[\pi ^{ab}\pi _{ab}{\displaystyle \frac{1}{2}}\pi _a^a\pi _b^b\right]\sqrt{\gamma }{}_{}{}^{3}R=0,`$ (14)
$`_a=2C_{ab}^d\pi ^{bc}\gamma _{cd}=0.`$ (15)
Here the canonical variables $`(N,P),(N^a,P_a),`$ and $`(\gamma _{ab},\pi ^{ab})`$ obey the following nonvanishing fundamental Poisson brackets relations:

$`\{N,P\}=1,`$ (16)
$`\{N^a,P_b\}=\delta _b^a,`$ (17)
$`\{\gamma _{ab},\pi ^{cd}\}={\displaystyle \frac{1}{2}}(\delta _a^c\delta _b^d+\delta _a^d\delta _b^c).`$ (18)
Due to the reparametrization symmetry of Eq. (9) inherited from the diffeomorphism invariance of the initial Hilbert action, the evolution of the system is unambiguous and it is governed by the total Hamiltonian

$$H_T=N+N^a_a+uP+u^aP_a,$$
(19)
with four arbitrary functions $`u(t)`$ and $`u^a(t)`$. One can verify that the secondary constraints are first class and obey the algebra

$`\{,_b\}=0,`$ (20)
$`\{_a,_b\}=C_{ab}^d_d.`$ (21)
Thus the dynamics of system (9) in the extended phase space, spanned by the canonical variables $`(N,P),(N^a,P_a)`$, and $`(\gamma _{ab},\pi ^{ab})`$, is described by the action rewritten in the Hamiltonian form as

$$𝒜=\mathrm{\Theta }+d(\pi ^{ab}\gamma _{ab}),$$
(22)
where $`\mathrm{\Theta }`$ is the Poincaré-Cartan 1-form

$$\mathrm{\Theta }=\pi ^{ab}d\gamma _{ab}+PdN+P_adN^aH_Tdt.$$
(23)
## 3 Deparametrized version of the model

### 3.1 Misner-Ryan decomposition

Let us now specialize to the Bianchi type I model. In this case, the three-dimensional group of isometries is an Abelian one, i.e., $`C_{ab}^d=0`$, and thus the momentum constraints (15) are identically satisfied, while the energy constraint reduces to the simplest form

$$=\pi ^{ab}\pi _{ab}\frac{1}{2}\pi _a^a\pi _b^b.$$
(24)
To find the corresponding unconstrained Hamiltonian system with a certain observable time parameter, we shall use the variables introduced by Misner  and modified by Ryan . In Misner’s representation, the spatial metric $`\gamma `$ is given by

$$\gamma _{ij}=r_0^2e^{2\mathrm{\Omega }}e_{ij}^{2\beta },$$
(25)
where $`\beta _{ij}`$ is a $`3\times 3`$ symmetric traceless matrix and $`\mathrm{\Omega }`$ is a scalar, both being functions of the time parameter only, and $`r_0`$ is a constant. 555 In Eq. (25), the constant $`r_0`$ was introduced for convenience in choosing units. For notational convenience, in this paper the constant $`r_0`$ is scaled to unity hereafter. For our purposes it is convenient, following , to pass to the main-axes decomposition for the nondegenerate symmetric matrix $`\beta `$ and finally write the metric $`\gamma `$ in the Misner-Ryan form as

$$\gamma =R^T(\chi )\left(\begin{array}{ccc}e^{\beta _++\sqrt{3}\beta _{}}& 0& 0\\ 0& e^{\beta _+\sqrt{3}\beta _{}}& 0\\ 0& 0& e^{2\beta _+}\end{array}\right)R(\chi ).$$
(26)
Here $`R(\chi )`$ is an orthogonal $`SO(3,)`$ matrix parametrized with the three Euler angles $`\chi _i`$. From the Jacobian of the transformation (26),

$$J\left(\frac{\gamma _{ab}(\mathrm{\Omega },\beta ,\chi )}{\mathrm{\Omega },\beta ,\chi }\right)e^{6\mathrm{\Omega }}\mathrm{sinh}(2\sqrt{3}\beta _{})\mathrm{sinh}(3\beta _++\sqrt{3}\beta _{})\mathrm{sinh}(3\beta _+\sqrt{3}\beta _{}),$$
(27)
it follows that Eq. (26) can be used as a definition of the new configuration variables, ($`\mathrm{\Omega },\beta _\pm )`$, and the three angles $`(\chi _1,\chi _2,\chi _3`$) only if all eigenvalues of the matrix $`\gamma `$ are different. To have a uniqueness of the inverse transformation, we are forced to treat only this type of configurations, the so-called principle orbits of the action of the $`SO(3,)`$ group, while the analysis of the orbits with coinciding eigenvalues of the matrix $`S`$ (singular orbits) requires special consideration.

The point transformation (26) induces the canonical transformation from the coordinates $`(\gamma _{ab},\pi ^{ab})`$ to the six new canonical pairs $`(\mathrm{\Omega },p_\mathrm{\Omega })`$, $`(\beta _\pm ,p_\pm )`$, and $`(\chi _i,p_{\chi _j})`$,

$`\{\mathrm{\Omega },p_\mathrm{\Omega }\}=1,\{\beta _\pm ,p_\pm \}=1,\{\chi _i,p_{\chi _j}\}=\delta _{ij}.`$ (28)
Using the requirement of invariance of the symplectic 1-form,

$$\pi ^{ab}d\gamma _{ab}=p_\mathrm{\Omega }d\mathrm{\Omega }+p_+d\beta _++p_{}d\beta _{}+\underset{i=1}{\overset{3}{}}p_{\chi _i}dp_{\chi _i},$$
(29)
one can find the explicit formulas for the new variables as functions of the old ones and finally rewrite the energy constraint as

$`2={\displaystyle \frac{1}{12}}(p_\mathrm{\Omega }^2+p_+^2+p_{}^2)`$
$`+{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{\eta _1^2}{\mathrm{sinh}^2(3\beta _+\sqrt{3}\beta _{})}}+{\displaystyle \frac{\eta _2^2}{\mathrm{sinh}^2(3\beta _++\sqrt{3}\beta _{})}}+{\displaystyle \frac{\eta _3^2}{\mathrm{sinh}^2(2\sqrt{3}\beta _{})}}\right).`$ (30)
Note that in this expression all angular variables are gathered into the three right-invariant $`SO(3,)`$ Killing vector fields $`\eta _i`$,

$`\eta _1=\mathrm{sin}\chi _1\mathrm{cot}\chi _2p_{\chi _1}+\mathrm{cos}\chi _1p_{\chi _2}+{\displaystyle \frac{\mathrm{sin}\chi _1}{\mathrm{sin}\chi _2}}p_{\chi _3},`$ (31)
$`\eta _2=\mathrm{cos}\chi _1\mathrm{cot}\chi _2p_{\chi _1}+\mathrm{sin}\chi _1p_{\chi _2}{\displaystyle \frac{\mathrm{cos}\chi _1}{\mathrm{sin}\chi _2}}p_{\chi _3},`$ (32)
$`\eta _3=p_{\chi _1},`$ (33)
satisfying the $`SO(3,)`$ Poisson brackets algebra

$$\{\eta _i,\eta _j\}=ϵ_{ijk}\eta _k.$$
(34)
### 3.2  Deparametrization and correspondence to the hyperbolic three-particle Euler-Calogero-Sutherland model

In this section, it will be shown that the unconstrained Hamiltonian, describing the evolution of observables in Bianchi type I cosmology, can be identified with the Hamiltonian of the relative motion in the three-particle hyperbolic Euler-Calogero-Sutherland system.

The direct method to “deparametrize” the initial reparametrization-invariant theory with first-class constraints consists of two steps, namely fixation of the gauge and projection of the initial degenerate action to the constraint shell. Then the expression for the reduced, unconstrained Hamiltonian can be extracted from the Poincaré-Cartan 1-form. Here it is in order to make one general comment on the gauge-fixing when some preferred set of variables is chosen. It turns out that in this case the action projected on the constraint shell becomes independent of the part of the variables without imposing some gauge-fixing condition. These eliminated variables should naturally be treated as unphysical ones. Such a phenomenon is well-known from the theory of Hamiltonian systems possessing a Lie group symmetry; when the cordinates adapted to the symmetry group action are found, then the so-called cyclic (or ignorable) coordinates disappear from the effective Hamiltonian. In the reparametrization-invariant theories, there is a certain peculiarity of this effect. Namely, when such a set of well-defined variables exists, then after reduction, one of these cyclic variables takes the role of time in the corresponding deparametrized system. The existence of such an exceptional variable in the reparametrization-invariant theories means that the equivalent unconstrained system can be written in an autonomous form.

Below we shall see that the model under consideration can be brought into a class of such Hamiltonian systems. 666 The corresponding representation for the Friedmann cosmological model with scalar and spinor fields has been discussed recently in . To see how the proper variable which takes the role of time in the unconstrained system appears, let us consider the energy constraint (3.1) in terms of the Misner-Ryan coordinates. One can easily verify that performing the canonical transformation

$$T=6\frac{\mathrm{\Omega }}{P_\mathrm{\Omega }},\mathrm{\Pi }_T=\frac{1}{12}P_\mathrm{\Omega }^2,$$
(35)
the energy constraint (3.1) linearizes with respect to the momentum $`\mathrm{\Pi }_T`$, 777 Such a linearity in one of the momenta forms of the constraints in accordance with Kuchař’s observation  shows the possibility of a global deparametrization of the degenerate theory.

$$2=\mathrm{\Pi }_T+h(\beta _\pm ,p_\pm ,\eta _i),$$
(36)
where $`h(\beta _\pm ,p_\pm ,\eta _i)`$ is given by

$`h(\beta _\pm ,p_\pm ,\eta _i)={\displaystyle \frac{1}{12}}(p_+^2+p_{}^2)`$
$`+{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{\eta _1^2}{\mathrm{sinh}^2(3\beta _+\sqrt{3}\beta _{})}}+{\displaystyle \frac{\eta _2^2}{\mathrm{sinh}^2(3\beta _++\sqrt{3}\beta _{})}}+{\displaystyle \frac{\eta _3^2}{\mathrm{sinh}^2(2\sqrt{3}\beta _{})}}\right).`$ (37)
Note that after all of these transformations, the Poincaré-Cartan 1-form (23) changes by a total differential,

$$\mathrm{\Theta }=\mathrm{\Pi }_TdT+p_+d\beta _++p_{}d\beta _{}+\underset{i=1}{\overset{3}{}}p_{\chi _i}d\chi _i+PdN+P_adN^aH_Tdt\frac{1}{2}d(P_\mathrm{\Omega }\mathrm{\Omega }).$$
(38)
Now to find the evolution parameter and the Hamiltonian of the deparametrized theory, we project Eq. (38) to the constraint shell $`P=P_a==_a=0`$, 888 We omit all differentials assuming that proper boundary conditions are taken to derive the classical equation of motion.

$$\mathrm{\Theta }|_{constraints}=p_+d\beta _++p_{}d\beta _{}+\underset{i=1}{\overset{3}{}}p_{\chi _i}d\chi _i+h(\beta _\pm ,p_\pm ,\eta _i)dT.$$
(39)
The expression for the projected 1-form (39) prompts the following interpretation. Supposing that the deparametrized version of the Bianchi type I model has a conventional Poincaré-Cartan 1-form,

$$\mathrm{\Theta }_{BI}^{}=p_+d\beta _++p_{}d\beta _{}+\underset{i=1}{\overset{3}{}}p_{\chi _i}d\chi _iH_{BI}(\beta _\pm ,p_\pm ,\eta _i)d\tau ,$$
(40)
we make the following identifications: the initial dynamical variable $`T(t)`$ “metamorphosis” to the evolution parameter $`\tau =T`$ of the corresponding deparametrized system and the Hamiltonian, as a time-independent generator of evolution, is given by $`H_{BI}=h(\beta _\pm ,p_\pm ,\eta _i)`$. Thus to summarize, we arrive at the unconstrained Hamiltonian for the Bianchi type I model in the form

$`H_{BI}={\displaystyle \frac{1}{12}}(p_+^2+p_{}^2)`$
$`+{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{\eta _1^2}{\mathrm{sinh}^2(3\beta _+\sqrt{3}\beta _{})}}+{\displaystyle \frac{\eta _2^2}{\mathrm{sinh}^2(3\beta _++\sqrt{3}\beta _{})}}+{\displaystyle \frac{\eta _3^2}{\mathrm{sinh}^2(2\sqrt{3}\beta _{})}}\right).`$ (41)
One final comment is in order concerning the form of the derived Hamiltonian. As was mentioned in the Introduction, the appearance of the square-root expression for the unconstraned Hamiltonian of the Bianchi cosmologies is usually interpreted in a certain analogy to relativistic particles . In contrast to this picture, our analysis points to the existence of a certain correspondence with a system of nonrelativistic particles.

To describe this correspondence, let us recall the formulation of the integrable three-particle hyperbolic Euler-Calogero-Sutherland (ECS) model . The particles are characterized by two types of variables. The first ones are the canonical coordinates $`(\beta _i,p_i)`$ describing the position and momenta,

$$\{\beta _i,p_j\}=\delta _{ij},$$
(42)
and the other are the so-called “internal” coordinates $`l_{ab}=l_{ba}`$, obeying the $`SO(3,)`$ Poisson brackets algebra,

$$\{l_{ab},l_{cd}\}=\delta _{ac}l_{bd}\delta _{ad}l_{bc}+\delta _{bd}l_{ac}\delta _{bc}l_{ad}.$$
(43)
The dynamics of both variables is determined by the Hamiltonian

$$H_{ECS}=\frac{1}{2}\underset{i=1}{\overset{3}{}}p_i^2+\frac{1}{8}\underset{(ij)}{}\frac{l_{ij}^2}{\mathrm{sinh}^2(\beta _i\beta _j)}.$$
(44)
Now it is easy to find the relation between the derived unconstrained Hamiltonian (3.2) and the Euler-Calogero-Sutherland Hamiltonian (44). To achieve this, we identify the internal variables $`l_{ij}`$ with the Killing vectors $`\eta _i`$ of the Bianchi type I model,

$$l_{ij}=ϵ_{ijk}\eta _k,$$
(45)
while the variables $`(\beta _\pm ,p_\pm )`$ are identified with the relative variables in the Jacobi system for $`(\beta _i,p_i)`$,

$`\beta _1=\beta _++\sqrt{3}\beta _{}X,`$ $`p_1={\displaystyle \frac{1}{6}}p_++{\displaystyle \frac{1}{2\sqrt{3}}}p_{}{\displaystyle \frac{1}{3}}P,`$ (46)
$`\beta _2=\beta _+\sqrt{3}\beta _{}X,`$ $`p_2={\displaystyle \frac{1}{6}}p_+{\displaystyle \frac{1}{2\sqrt{3}}}p_{}{\displaystyle \frac{1}{3}}P,`$ (47)
$`\beta _3=2\beta _+X,`$ $`p_3={\displaystyle \frac{1}{3}}p_+{\displaystyle \frac{1}{3}}P.`$ (48)
Here $`(X,P)`$ are the corresponding center-of-mass coordinates. Now one can verify that the expression (44) takes the following separable form:

$$H_{ECS}=\frac{1}{6}P^2+H_{BI}$$
(49)
and thus one arrives at the above-stated correspondence between the Hamiltonian, describing the relative motion in a three-particle hyperbolic Euler-Calogero-Sutherland system, and the deparametrized version of the Bianchi type I model Hamiltonian $`H_{BI}`$.

### 3.3 Solution of the unconstrained equations of motion

The established relation to the many-body integrable model allows us to write down the general solution for the Bianchi type I metric. Here it is in order to note that the solution of the initial Einstein equations for the Bianchi-type metrics contain arbitrary functions, reflecting the reparametrization invariance, of the theory and the solution should be classified according to their equivalence. Passing to the reduced system, when the constraints generating this reparametrization are eliminated, this freedom is fixed. However, if the system has the additional rigid symmetry, not all of solutions describe the nonequivalent metrics and now the problem of classification reduces to the identification of the metrics, with a fixed number of parameters. Below we shall present the general solution for the reduced Bianchi type I model, while the detailed analysis of the essential parameters and their origin will be done elsewhere.

To write down an explicit form of the solution for the general metric of Bianchi type I space-time, we recall the roots of the integrability of the Euler-Calogero-Sutherland model. It is known (see, e.g.,  and references therein) that the integral curves of the Hamiltonian (44) are in one-to-one correspondence with the geodesics of the space of symmetric matrices $`S`$ endowed with the bi-invariant metric,

$$L_{ECS}=\frac{1}{2}\text{Tr}(S^1dS)^2.$$
(50)
To prove this dynamical equivalence, it is enough to use again the main-axes decomposition for the symmetric matrix $`S`$ of the form

$$S=R^Te^{2X}R$$
(51)
with the diagonal matrix $`X=\text{diag}\beta _1,\beta _2,\beta _3`$. Assuming that $`S_{ij}`$ are six independent Lagrange coordinates, from the action (50) it follows that the equation of motion can be written in the form of conservation law,

$$\frac{d}{dt}\left(S^1\frac{d}{dt}S\right)=0,$$
(52)
and therefore the general solution for the matrix $`S`$ can be represented in a simple form,

$$S(t)=S_0e^{Jt}.$$
(53)
Here $`S_0`$ and $`J`$ are arbitrary constant matrices, $`S_0`$ is a symmetric matrix, the value of the matrix $`S(t)`$ at an initial moment $`t=0`$, while $`J`$ is an integral of motion, $`J=S_0^1\dot{S}_0`$. 999 One can easily check that the solution (53) is a symmetric matrix if its initial value $`S_0`$ ware a symmetric one. Comparison of the two expressions (51) and (53) gives the solution of the Euler-Calogero-Sutherland model and therefore of the Bianchi type I metric due to the correspondence derived in the previous section.

## 4 Concluding remarks

Following the Dirac generalized Hamiltonian approach we derive the deparametrized version of the Bianchi type I cosmology in the form of the integrable three-particle Euler-Calogero-Sutherland model. It is worth posing the question of the possibility to extend the result presented here and to find the quadratic time-independent unconstrained Hamiltonians for the other more complicated Bianchi systems. In seeking the answer, first of all it is necessary to investigate the rigid symmetries that Bianchi models possess, because the additional rigid symmetry as a rule leads to the conception of preferred time choice . More precisely, the existence of a timelike Killing vector field allows us to define an intrinsic time and construct the observables, whose dynamics is governed by a time-independent Hamiltonian. 101010 Note that sometimes it is enough to have even a conformal timelike Killing vector field; see, e.g., .

Another application of the discussed correspondence is to analyze the metric of Bianchi type I cosmology with respect to its dependence on the essential parameters. In a forthcoming publication, we plan to discuss in detail a certain parametrization of the metric using the solution (53) of the three-body Euler-Calogero-Sutherland problem and we shall attempt to clarify how the topology of the flat Bianchi type I cosmology emerges. One can find the recent discussions of the topological characteristics of the three-dimensional space in the cosmology in review .

## Acknowledgments

Helpful discussions during the work on the paper with H. Braden, R. Camassa, G. Dvali, P. Hajicek, G. Lavrelashvili, M.D. Mateev, V.P. Pavlov, O. Ragnisco, S. Ruijsenaars, and Y. Suris are acknowledged. A.K. would like to thank the Abdus Salam International Centre for Theoretical Physics (ICTP), Trieste, Italy for support and kind hospitality extended to him at ICTP High Energy Section. D.M. is grateful to the organizers, with special thanks to Orlando Ragnisco, of the workshop ”The Calogero-Moser System Thirty Years Later,” Roma, and to the participants for many interesting conversations and discussions. This work was supported in part by the RFBR grant No. 01-01-00708 and INTAS grant No. 00-00561.