format: text

text:
# Exact coherent states in one-dimensional quantum many-body systems with inverse-square interactions

## I Introduction

The harmonic oscillator (of time-dependent parameter) is a model where the path integral for the kernel (propagator) is Gaussian and thus the kernel can be almost determined by the classical action. From the fact that the kernel should satisfy Schrödinger equation, with a condition for the kernel in the coincident limit, one can obtain the exact expression of the kernel in terms of the solutions of classical equation of motion . This is one of the basic reasons of the fact that wave functions of harmonic oscillators are described by the solutions of classical equation of motion. Two different types of operators have been known for the construction of coherent states from the vacuum state of a simple harmonic oscillator system; the displacement operator and squeeze operator . Recently, it has been shown  that a harmonic oscillator of time-dependent parameters are unitarily equivalent to a simple harmonic oscillator of different timescale (including the case with an inverse-square potential ), where the unitary operator of the relation is again given in terms of the classical solutions of the system of time-dependent parameters. If one considers the set of wave functions whose centers of the probability distribution functions do not move, the operator at a given time corresponds to a squeeze operator. As have been well known, a driven harmonic oscillator (without an inverse-square potential) is unitarily equivalent to a harmonic oscillator without driving force, and this relation can be used in finding the unitary transformation which does not change the form of the Schrödinger equation in a harmonic oscillator. The unitary operator for this transformation corresponds to a displacement operator. By applying this transformation to the eigenstates of stationary probability distribution of the model, one can obtain the wave functions whose centers of the probability distributions move according to the the classical solution.

In this paper, we will consider the models of identical $`N`$-body harmonic oscillators of time-dependent parameters interacting through the potential $`V(x_1,x_2,\mathrm{},x_N)`$ of homogeneous degree -2 which satisfies

$$V(ax_1,ax_2,\mathrm{},ax_N)=a^2V(x_1,x_2,\mathrm{},x_N)$$
(1)
with non-zero constant $`a`$. One of the models which belongs to this category was first solved by Sutherland , based on the earlier work by Calogero . From its inception, this model is closely related to the random matrix model  and has been found relevant for the descriptions of various physical phenomena . The model has been generalized into the several cases , known in the literature as the Calogero-Sutherland models, which have generated wide interest , while some of the generalized models do not belong to the category we will consider (see Sec. V.).

We will show that, there is a unitary transformation which relates the system to the the model of different timescale with unit mass and spring constant. The operator for the transformation is given as a product of the squeeze-type unitary operators in (one-body) harmonic oscillators. If $`V(x_1,x_2,\mathrm{},x_N)`$ is written in terms of the differences between positions of two particles, so that

$$V(x_1+a,x_2+a,\mathrm{},x_N+a)=V(x_1,x_2,\mathrm{},x_N),$$
(2)
there also exists a unitary transformation which does not change the form of the Schrödinger equation, while the operator for the transformation is given as the product of the displacement-type unitary operators in (one-body) harmonic oscillators. Both unitary operators are symmetric under the exchange of any pair of particles. These operators, as in the one-body harmonic oscillator case , thus can be used to find the wave functions of the systems of time-dependent parameters from those of constant parameters, preserving the symmetric property of the wave functions. The unitary transformations will be explicitly applied for three cases in the Calogero-Sutherland models, to find exact wave functions. In a Sutherland model of time-dependent parameters, Sutherland  has found an exact ”coherent” state by directly analyzing the Schrödinger equation; we will show that this ”coherent” state is obtained by applying the squeeze-type operator, while the displacement-type operator can also be applied in this model to give a more general form of the exact wave function. In addition, by extending the definition of the displacement-type operator, we will also show that one can find a unitary relation between the system without external force and the same system with external force. Since the quantum states we will find through unitary transformations can be solely described, up to some parameters of the models, by the classical solutions of a harmonic oscillator, we will call them coherent states.

This paper will be organized as follows; in the next section, we will introduce the squeeze-type unitary operator which relates the interacting $`N`$-body oscillator system of time-dependent parameters and the same system of constant parameters. The displacement-type operator will be also introduced, and it will be shown that, if the potential of interaction satisfies Eq. (2), the unitary transformation by the operator does not change the form of the Schrödinger equation. In Sec. III, the unitary transformations will be explicitly applied for three cases in the Calogero-Sutherland models. In Sec. IV, for the cases that the potential of the interaction satisfies Eq. (2), it will be shown that the displacement-type operator can be extended to give the unitary relation between a system without external force and the same system with external force. The last section will be devoted to a summary and discussions.

## II The unitary transformations

It has been shown that , the harmonic oscillator with inverse-square potential described by the Hamiltonian

$$H_{s,1}=\frac{p^2}{2}+\frac{x^2}{2}+\frac{g}{x^2}$$
(3)
with a coupling constant $`g`$, and the oscillator with time-dependent mass $`M(t)`$, spring constant $`w(t)`$ described by the Hamiltonian

$$H_1=\frac{p^2}{2M(t)}+\frac{1}{2}M(t)w^2(t)x^2+\frac{g}{M(t)}\frac{1}{x^2}$$
(4)
are related through the unitary transformation. For the case of $`g=0`$, the classical equation of motion for the system of the Hamiltonian in Eq. (4) is written as;

$$\frac{d}{dt}(M\dot{x})+M(t)w^2(t)x=0.$$
(5)
If we denote the two linearly independent solutions of Eq. (5) as $`u(t)`$ and $`v(t)`$, the $`\rho (t)`$ defined by $`\rho (t)=\sqrt{u^2+v^2}`$ satisfies

$$\frac{d}{dt}(M\dot{\rho })\frac{\mathrm{\Omega }^2}{M\rho ^3}+Mw^2\rho =0$$
(6)
with a time-constant $`\mathrm{\Omega }`$ ($`M(t)[\dot{v}(t)u(t)\dot{u}(t)v(t)]`$), while the overdots denote differentiation with respect to $`t`$. Without losing generality we assume that $`\mathrm{\Omega }`$ is positive. By defining the operator $`O_{s,1}`$ and $`O_1`$ as

$`O_{s,1}(\tau )=i\mathrm{}{\displaystyle \frac{}{\tau }}+H_{s,1}`$ (7)
$`O_1(t)=i\mathrm{}{\displaystyle \frac{}{t}}+H_1,`$ (8)
if the $`\tau `$, the time of the system of $`H_{s,1}`$, and $`t`$, the time of the system of $`H_1`$, is related as

$$d\tau =\frac{\mathrm{\Omega }}{\rho ^2}dt,$$
(9)
the unitary relation between the two operators has been given in Ref.  as

$$U_1O_{s,1}(\tau )U_1^{}_{\tau =\tau (t)}=\frac{M\rho ^2}{\mathrm{\Omega }}O_1,$$
(10)
where

$$U_1=\mathrm{exp}[\frac{i}{2\mathrm{}}M\frac{\dot{\rho }}{\rho }x^2]\mathrm{exp}[\frac{i}{4\mathrm{}}\mathrm{ln}(\frac{\rho ^2}{\mathrm{\Omega }})(xp+px)].$$
(11)
There is a similar unitary relation between the identical $`N`$-body harmonic oscillators interacting through the potential $`V(x_1,x_2,\mathrm{},x_N)`$ satisfying Eq. (1). The potential $`V`$ may be written as a linear combination of the terms $`1/_{l,m=1}^Na_{lm}x_lx_m`$, where $`a_{lm}=a_{ml}`$. If a system is described by the Hamiltonian

$$H_{s,N}=\underset{i=1}{\overset{N}{}}(\frac{p_i^2}{2}+\frac{x_i^2}{2})+V(x_1,x_2,\mathrm{},x_N)$$
(12)
and another system is described by the Hamiltonian

$`H_N`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}({\displaystyle \frac{p_i^2}{2M(t)}}+M(t)w^2(t){\displaystyle \frac{x_i^2}{2}})`$ (14)
$`+{\displaystyle \frac{1}{M(t)}}V(x_1,x_2,\mathrm{},x_N),`$
from the commutator relation

$`[{\displaystyle \underset{i}{\overset{N}{}}}(x_ip_i+p_ix_i),1/{\displaystyle \underset{l,m=1}{\overset{N}{}}}a_{lm}x_lx_m]`$ (15)
$`=4i\mathrm{}/{\displaystyle \underset{l,m=1}{\overset{N}{}}}a_{lm}x_lx_m,`$ (16)
one may find the unitary relation of the two systems

$$U_NO_{s,N}(\tau )U_N^{}_{\tau =\tau (t)}=\frac{M\rho ^2}{\mathrm{\Omega }}O_N,$$
(17)
where

$`O_{s,N}(\tau )`$ $`=`$ $`i\mathrm{}{\displaystyle \frac{}{\tau }}+H_{s,N}`$ (18)
$`O_N(t)`$ $`=`$ $`i\mathrm{}{\displaystyle \frac{}{t}}+H_N.`$ (19)
This relation has been noticed for a specific case . In Eq. (14), the unitary operator $`U_N`$ is given as

$`U_N`$ (20)
$`={\displaystyle \underset{i=1}{\overset{N}{}}}(\mathrm{exp}[{\displaystyle \frac{i}{2\mathrm{}}}M{\displaystyle \frac{\dot{\rho }}{\rho }}x_i^2]\mathrm{exp}[{\displaystyle \frac{i}{4\mathrm{}}}\mathrm{ln}({\displaystyle \frac{\rho ^2}{\mathrm{\Omega }}})(x_ip_i+p_ix_i)])`$ (21)
$`=({\displaystyle \frac{\mathrm{\Omega }}{\rho ^2}})^{N/4}{\displaystyle \underset{i=1}{\overset{N}{}}}(\mathrm{exp}[{\displaystyle \frac{i}{2\mathrm{}}}M{\displaystyle \frac{\dot{\rho }}{\rho }}x_i^2]\mathrm{exp}[{\displaystyle \frac{1}{2}}\mathrm{ln}({\displaystyle \frac{\rho ^2}{\mathrm{\Omega }}})x_i{\displaystyle \frac{}{x_i}}]).`$ (22)
If $`\varphi _s(x_1,x_2,\mathrm{},x_N)`$ is an eigenstate of Hamiltonian $`H_{s,N}`$ with eigenvalue $`E`$, from Eqs. (14,18), the wave function $`\psi (t;x_1,x_2,\mathrm{},x_N)`$ satisfying $`O_N(t)\psi =0`$ is given as

$`\psi (t;x_1,x_2,\mathrm{},x_N)`$ (23)
$`=e^{iE\tau /\mathrm{}}_{\tau =\tau (t)}U_N\varphi _s(x_1,x_2,\mathrm{},x_N)`$ (24)
$`=({\displaystyle \frac{\mathrm{\Omega }}{\rho ^2}})^{N/4}({\displaystyle \frac{u(t)iv(t)}{\rho (t)}})^{E/\mathrm{}}({\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{exp}[{\displaystyle \frac{i}{2\mathrm{}}}M{\displaystyle \frac{\dot{\rho }}{\rho }}x_i^2])`$ (25)
$`\times \varphi _s(\sqrt{{\displaystyle \frac{\mathrm{\Omega }}{\rho ^2}}}x_1,\sqrt{{\displaystyle \frac{\mathrm{\Omega }}{\rho ^2}}}x_2,\mathrm{},\sqrt{{\displaystyle \frac{\mathrm{\Omega }}{\rho ^2}}}x_N).`$ (26)
If $`V(x_1,x_2,\mathrm{},x_N)`$ is written in terms of the differences of positions of two particles, so that $`V`$ satisfies Eq.(2), there is a unitary operator

$$U_f=e^{\frac{i}{\mathrm{}}N\delta _f}\underset{i=1}{\overset{N}{}}(\mathrm{exp}[\frac{i}{\mathrm{}}M\dot{u}_fx_i]\mathrm{exp}[\frac{i}{\mathrm{}}u_fp_i])$$
(27)
which does not change the $`O_N`$ under a unitary transformation:

$$U_fO_NU_f^{}=O_N.$$
(28)
In Eq. (21), $`u_f`$ is a linear combination of $`u(t),v(t)`$, and $`\delta _f`$ is defined through the relation

$$\dot{\delta }_f=\frac{1}{2}M(w^2u_f^2\dot{u}_f^2).$$
(29)
Therefore, in this case, a coherent wave function from $`\varphi _s(x_1,x_2,\mathrm{},x_N)`$ is given as

$`\psi ^f(t;x_1,,x_2,\mathrm{},x_N)`$ (30)
$`=e^{iE\tau /\mathrm{}}_{\tau =\tau (t)}U_fU_N\varphi _s(x_1,x_2,\mathrm{},x_N)`$ (31)
$`=({\displaystyle \frac{\mathrm{\Omega }}{\rho ^2}})^{N/4}({\displaystyle \frac{uiv}{\rho }})^{E/\mathrm{}}e^{\frac{i}{\mathrm{}}N\delta _f}`$ (32)
$`\times ({\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{exp}[{\displaystyle \frac{i}{2\mathrm{}}}M{\displaystyle \frac{\dot{\rho }}{\rho }}(x_iu_f)^2+{\displaystyle \frac{i}{\mathrm{}}}M\dot{u}_fx_i])`$ (33)
$`\times \varphi _s(\sqrt{{\displaystyle \frac{\mathrm{\Omega }}{\rho ^2}}}(x_1u_f),\mathrm{},\sqrt{{\displaystyle \frac{\mathrm{\Omega }}{\rho ^2}}}(x_Nu_f)).`$ (34)
From the derivations through unitary transformations, it is manifest that

$$\underset{i=1}{\overset{N}{}}dx_i\varphi _s^{}\varphi _s=\underset{i=1}{\overset{N}{}}dx_i\psi ^{}\psi =\underset{i=1}{\overset{N}{}}dx_i\psi _{}^{f}{}_{}{}^{}\psi ^f.$$
(35)
Since the two unitary operators, $`U_N,U_f`$ are invariant under the exchange of a pair of $`i`$-th and $`j`$-th particles, the wave functions $`\varphi _s(x_1,x_2,\mathrm{},x_N)`$, $`\psi (t;x_1,x_2,\mathrm{},x_N)`$ and $`\psi ^f(t;x_1,x_2,\mathrm{},x_N)`$ have the same symmetric property under the exchanges of particles. For the systems of identical particles, one of the quantities of interest which is independent of statistics  is the particle number density defined for $`\varphi _s`$ as

$$\sigma _s(x)=N\frac{_{\mathrm{}}^{\mathrm{}}𝑑x_2\mathrm{}_{\mathrm{}}^{\mathrm{}}𝑑x_N\varphi _s^2(x,x_2,\mathrm{},x_N)}{_{\mathrm{}}^{\mathrm{}}𝑑x_1\mathrm{}_{\mathrm{}}^{\mathrm{}}𝑑x_N\varphi _s^2(x_1,x_2,\mathrm{},x_N)}.$$
(36)
From the Eq. (25), one can easily find the expression of the particle number density for $`\psi ^f`$ as

$$\sigma ^f(x)=\frac{\sqrt{\mathrm{\Omega }}}{\rho }\sigma _s(\frac{\sqrt{\mathrm{\Omega }}}{\rho }(xu_f)).$$
(37)
## III Applications

In this section, the general results of previous section will be explicitly applied for three cases. First, we will consider the Sutherland model of Ref. . We will find a general expression of a coherent state of the system (of time-dependent parameters) and will show that the expression reproduces the known state in the model . Second, we will consider the model of three-body system . Third, we will consider the Calogero model ”in the Jacobi coordinate”  without the degree of freedom of center of mass.

### A Sutherland model

The system described by the Hamiltonian

$$H_{S,s}=\underset{i=1}{\overset{N}{}}(\frac{p_i^2}{2}+\frac{x_i^2}{2})+\underset{i>j=1}{\overset{N}{}}\frac{\mathrm{}^2\lambda (\lambda 1)}{(x_ix_j)^2}$$
(38)
has the (unnormalized) bosonic ground state

$$\varphi _S=(\underset{j>i=1}{\overset{N}{}}|x_jx_i|^\lambda )\underset{i}{\overset{N}{}}e^{x_i^2/2\mathrm{}}$$
(39)
with the energy eigenvalue $`\mathrm{}N[1+\lambda (N1)]/2`$. For the system described by the Hamiltonian

$`H_S`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}({\displaystyle \frac{p_i^2}{2M(t)}}+M(t)w^2(t){\displaystyle \frac{x_i^2}{2}})`$ (41)
$`+{\displaystyle \frac{1}{M(t)}}{\displaystyle \underset{i>j=1}{\overset{N}{}}}{\displaystyle \frac{\mathrm{}^2\lambda (\lambda 1)}{(x_ix_j)^2}},`$
the wave function $`\psi _S^f`$ satisfying

$$i\mathrm{}\frac{\psi _S^f}{t}=H_S\psi _S^f$$
(42)
is given, from the results of previous section, as

$`\psi _S^f`$ $`=`$ $`({\displaystyle \frac{u+iv}{\sqrt{\mathrm{\Omega }}}})^{N(1+\lambda (N1))/2}e^{iN\delta _f/\mathrm{}}`$ (45)
$`\times ({\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{exp}[{\displaystyle \frac{i}{2\mathrm{}}}M{\displaystyle \frac{\dot{\rho }}{\rho }}(x_iu_f)^2+{\displaystyle \frac{i}{\mathrm{}}}M\dot{u}_fx_i])`$
$`\times ({\displaystyle \underset{j>i=1}{\overset{N}{}}}|x_jx_i|^\lambda ){\displaystyle \underset{i=1}{\overset{N}{}}}e^{\mathrm{\Omega }(x_iu_f)^2/2\mathrm{}\rho ^2}.`$
In deriving Eq. (33), we make use of the fact that $`(uiv)/\rho ^2=1/(u+iv)`$. By making use of Eq. (28) and the results in Ref. , one can find the particle density for $`\psi _S^f`$ is give as

$$\sigma _S^f(x)=\frac{\sqrt{2N\mathrm{\Omega }}}{\pi \rho \sqrt{\lambda }}\sqrt{1\frac{\mathrm{\Omega }}{2N\lambda \rho ^2}(xu_f)^2}.$$
(46)
By directly analyzing the Schrödinger equation (Eq. (32)), for the unit mass case, a coherent wave function is given in terms of a complex solution of classical equation of motion, and a hydrodynamic description was shown to hold exactly in the picture that the wave function provides . If we choose $`u_f=0`$, one can easily verify that, for the unit mass case, the $`\psi _S^f`$ and $`\sigma _S^f`$ reduce to the wave function and the density found by Sutherland , respectively.

As in the case of one-body harmonic oscillator , there are, in general, five free parameters in determining $`\psi _S^f`$ or $`\sigma _S^f`$; two of the parameters determine the motion of center of the particle number density, while the other three parameters determine the shape of the density function. To be explicit, we consider the particle density function of the case of unit mass and unit spring constant. In this case, two homogeneous solutions $`u(t),v(t)`$ and the (fictitious) particular solution $`u_f(t)`$ can be taken, without losing generality, as $`\mathrm{cos}(t+t_0),A\mathrm{sin}(t+\alpha +t_0)`$ and $`B\mathrm{cos}(t+\beta )`$, respectively, with real constants $`t_0`$, $`\beta `$, positive constants $`A`$, $`B`$, and a real constant $`\alpha `$ satisfying $`|\alpha |<\pi `$. Then the density function is written as

$$\frac{\sqrt{2NA\mathrm{cos}\alpha }}{\pi \stackrel{~}{\rho }\sqrt{\lambda }}\sqrt{1\frac{A\mathrm{cos}\alpha }{2N\lambda \stackrel{~}{\rho }^2}(xB\mathrm{cos}(t+\beta ))^2},$$
(47)
with

$$\stackrel{~}{\rho }=\sqrt{\mathrm{cos}^2(t+t_0)+A^2\mathrm{sin}^2(t+\alpha +t_0)}.$$
(48)
Due to the time-translational invariance, in this case, one of the parameters is simply related to the time shifting of the density functions.

### B Three-body interaction model

As another example, we consider the three-body interaction model described by the Hamiltonian :

$`H_{3body,s}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{3}{}}}({\displaystyle \frac{p_i^2}{2}}+{\displaystyle \frac{x_i^2}{2}})`$ (50)
$`+{\displaystyle \underset{i>j=1}{\overset{3}{}}}{\displaystyle \frac{\mathrm{}^2\lambda (\lambda 1)}{(x_ix_j)^2}}+3{\displaystyle \underset{i>j=1}{\overset{3}{}}}{\displaystyle \frac{\mathrm{}^2\alpha (\alpha 1)}{y_{ij}^2}},`$
where $`y_{ij}`$ is defined as $`x_i+x_j2x_k`$ ($`ki`$ and $`kj`$). The (unnormalized) bosonic ground state is given as

$$\varphi _{3body}=\underset{i>j=1}{\overset{3}{}}(|x_ix_j|^\lambda |y_{ij}|^\alpha )\underset{i=1}{\overset{3}{}}e^{x_i^2/2\mathrm{}}$$
(51)
with energy eigenvalue $`3\mathrm{}(\frac{1}{2}+(\lambda +\alpha ))`$. One can easily find that the potential of mutual interaction in $`H_{3body,s}`$ is written in terms of the differences between positions of two particles. By applying the formulas in the previous section, one can find that the exact wave function for the system described by the Hamiltonian

$`H_{3body,s}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{3}{}}}({\displaystyle \frac{p_i^2}{2M(t)}}+M(t)w^2(t){\displaystyle \frac{x_i^2}{2}})`$ (54)
$`+{\displaystyle \frac{1}{M(t)}}{\displaystyle \underset{i>j=1}{\overset{3}{}}}{\displaystyle \frac{\mathrm{}^2\lambda (\lambda 1)}{(x_ix_j)^2}}`$
$`+{\displaystyle \frac{1}{M(t)}}{\displaystyle \underset{i>j=1}{\overset{3}{}}}{\displaystyle \frac{\mathrm{}^2\alpha (\alpha 1)}{y_{ij}^2}},`$
is given as

$`\psi _{3body}^f`$ $`=`$ $`({\displaystyle \frac{u+iv}{\sqrt{\mathrm{\Omega }}}})^{3(\lambda +\alpha )3/2}e^{3i\delta _f/\mathrm{}}`$ (57)
$`\times {\displaystyle \underset{i=1}{\overset{3}{}}}\mathrm{exp}[{\displaystyle \frac{i}{2\mathrm{}}}M{\displaystyle \frac{\dot{\rho }}{\rho }}(x_iu_f)^2+{\displaystyle \frac{i}{\mathrm{}}}M\dot{u}_fx_i]`$
$`\times ({\displaystyle \underset{i>j=1}{\overset{3}{}}}|x_ix_j|^\lambda |y_{ij}|^\alpha ){\displaystyle \underset{i=1}{\overset{3}{}}}e^{\mathrm{\Omega }(x_iu_f)^2/2\mathrm{}\rho ^2}.`$
### C Calogero model in Jacobi coordinates

The model described by the Hamiltonian

$$\stackrel{~}{H}_{C,s}=\underset{i=1}{\overset{N}{}}\frac{p_i^2}{2}+\frac{1}{2N}\underset{i>j=1}{\overset{N}{}}(x_ix_j)^2+V_C$$
(58)
has been considered by Calogero , where

$$V_C=\underset{i>j=1}{\overset{N}{}}\frac{\mathrm{}^2\lambda (\lambda 1)}{(x_ix_j)^2}.$$
(59)
If one introduce the ”Jacobi coordinates”

$`y_i={\displaystyle \frac{1}{\sqrt{i(i+1)}}}({\displaystyle \underset{l=1}{\overset{i}{}}}x_lix_{i+1})(i=1,2,\mathrm{},N1),`$ (60)
$`y_N={\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{l=1}{\overset{N}{}}}x_l,`$ (61)
the Hamiltonian in Eq. (41) is written as

$$\stackrel{~}{H}_{C,s}=\frac{p_{y_N}^2}{2}+H_{C,s},$$
(62)
where

$$H_{C,s}=\underset{i=1}{\overset{N1}{}}(\frac{p_{y_i}^2}{2}+\frac{y_i^2}{2})+V_C.$$
(63)
It is easy to see that $`V_C`$ does not depend on $`y_N`$ and the Hamiltonian $`H_{C,s}`$ describes a system of interacting $`N1`$ particles. In fact, Calogero analyzed the Hamiltonian system of $`H_{C,s}`$, and found the (unnormalized) wave functions

$`\varphi _n^C(y_1,y_2,\mathrm{},y_{N1})`$ (64)
$`=({\displaystyle \underset{i>j=1}{\overset{N}{}}}(x_ix_j)^\lambda )\mathrm{exp}({\displaystyle \frac{1}{2\mathrm{}}}{\displaystyle \underset{i=1}{\overset{N1}{}}}y_i^2)L_n^b({\displaystyle \frac{1}{\mathrm{}}}{\displaystyle \underset{i=1}{\overset{N1}{}}}y_i^2)`$ (65)
satisfying

$$H_{C,s}\varphi _n^C=E_n\varphi _n^C(n=0,1,2,\mathrm{}),$$
(66)
where

$`b`$ $`=`$ $`{\displaystyle \frac{1}{2}}(N3)+{\displaystyle \frac{1}{2}}\lambda N(N1),`$ (67)
$`E_n`$ $`=`$ $`\mathrm{}[{\displaystyle \frac{1}{2}}(N1)+{\displaystyle \frac{1}{2}}\lambda N(N1)+2n],`$ (68)
and $`L_n^b`$ is the associated Laguerre polynomials defined through the equation

$$x\frac{d^2L_n^b}{dx^2}+(b+1x)\frac{dL_n^b}{dx}+nL_n^b(x)=0.$$
(69)
If $`y_i`$ is the space coordinate of the $`i`$-th particle, the Hamiltonian $`H_{C,s}`$ does not describe the system of identical particles, as has been explicitly shown in the 3-body system .

Since the $`V_C`$ can not be written in terms of $`y_iy_j`$, only squeeze-type unitary transformation can be applied to give coherent sates. For the system described by the Hamiltonian

$$H_C=\underset{i=1}{\overset{N1}{}}(\frac{p_{y_i}}{2M(t)}+M(t)w^2(t)\frac{y_i^2}{2})+\frac{V_C}{M(t)},$$
(70)
by making use of the unitary relation in Eq. (14), one may find that the wave functions satisfying $`i\mathrm{}(\psi _n^C/t)=H_C\psi _n^C`$ are given as

$`\psi _n^C`$ $`=`$ $`({\displaystyle \frac{u+iv}{\sqrt{\mathrm{\Omega }}}})^{b1}({\displaystyle \frac{uiv}{\rho }})^{2n}{\displaystyle \underset{i>j=1}{\overset{N}{}}}(x_ix_j)^\lambda `$ (72)
$`\times \mathrm{exp}[{\displaystyle \frac{1}{2\mathrm{}}}({\displaystyle \frac{\mathrm{\Omega }}{\rho ^2}}iM{\displaystyle \frac{\dot{\rho }}{\rho }}){\displaystyle \underset{i=1}{\overset{N1}{}}}y_i^2]L_n^b({\displaystyle \frac{\mathrm{\Omega }}{\mathrm{}\rho ^2}}{\displaystyle \underset{i=1}{\overset{N1}{}}}y_i^2).`$
## IV A Generalization to include external force

If $`V`$ can be written in terms of the differences between positions of two particles, by modifying the $`U_f`$, one can find the unitary relation in different Hamiltonian systems, as in one-body harmonic oscillator . By defining the $`U_F`$ as

$$U_F=e^{\frac{i}{\mathrm{}}N\delta _F}\underset{i=1}{\overset{N}{}}(\mathrm{exp}[\frac{i}{\mathrm{}}M\dot{x}_px_i]\mathrm{exp}[\frac{i}{\mathrm{}}x_pp_i])$$
(73)
where $`x_p`$ and $`\delta _F`$ are defined through the relations

$`{\displaystyle \frac{d}{dt}}(M\dot{x}_p)`$ $`+`$ $`M(t)w^2(t)x_p=F(t),`$ (74)
$`\dot{\delta }_F`$ $`=`$ $`{\displaystyle \frac{1}{2}}M(w^2x_p^2\dot{x}_p^2),`$ (75)
one can find the relation

$$U_FO_NU_F^{}=i\mathrm{}\frac{}{t}+H_{N,F}$$
(76)
where

$$H_{N,F}=H_NF\underset{i=1}{\overset{N}{}}x_i.$$
(77)
From Eq. (56), one can find that the wave function

$`\psi ^F`$ $`=`$ $`({\displaystyle \frac{\mathrm{\Omega }}{\rho ^2}})^{N/4}({\displaystyle \frac{uiv}{\rho }})^{E/\mathrm{}}e^{\frac{i}{\mathrm{}}N\delta _F}`$ (80)
$`\times ({\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{exp}[{\displaystyle \frac{i}{2\mathrm{}}}M{\displaystyle \frac{\dot{\rho }}{\rho }}(x_ix_p)^2+{\displaystyle \frac{i}{\mathrm{}}}M\dot{x}_px_i])`$
$`\times \varphi _s({\displaystyle \frac{\sqrt{\mathrm{\Omega }}}{\rho }}(x_1x_p),\mathrm{}{\displaystyle \frac{\sqrt{\mathrm{\Omega }}}{\rho }}(x_Nx_p))`$
satisfies the Schrödinger equation

$$i\mathrm{}\frac{\psi ^F}{t}=H_{N,F}\psi ^F.$$
(81)
## V Summary and discussions

It has been shown that the unitary relations in one-body harmonic oscillator systems can be extended to give the unitary relations in some of the Calogero-Sutherland models. These unitary relations can be used not only in finding coherent states from the given stationary states in a system, but also in finding exact wave functions of the Hamiltonian systems of time-dependent parameters from those of time-independent Hamiltonian systems. If the potential of mutual interactions can be written in terms of the differences between positions of two particles, we have also shown that the wave functions of the system with external force can be found from those of the the same system without the external force. The list of applications given in this paper is not exhaustive.

The unitary relations can be formally extended to the case of identical $`N`$-body free particles interacting through the mutual interaction potential $`V`$; however, in this case, the $`\rho (t)`$ diverges as $`t^2`$ goes to infinity. Even for the case of identical $`N`$-body harmonic oscillators (interacting through the $`V`$), $`\rho (t)`$ could be unbounded in general, as analyzed in detail for the cases of periodic mass and frequencies .

The system of identical $`N`$-body free particles interacting through a potential $`\frac{\mathrm{}\lambda (\lambda 1)\pi ^2}{L^2}_{i>j}^N\frac{1}{\mathrm{sin}^2[\pi (x_ix_j)/L\sqrt{\mathrm{}}]}`$ with a constant $`L`$ , has recently been of great interest , while the potential is not of homogeneous degree -2. Since the potential satisfies Eq. (2), for the unit mass case, if one choose $`u_f`$ as $`at+b`$ with constants $`a,b`$, one can show that the Schrödinger equation is invariant under a unitary transformation, as in Eq. (22) (with $`w=0`$). For this case, the (unnormalized) wave function of the ground state is given as $`\psi _0=_{i>j}(\mathrm{sin}\frac{\pi (x_ix_j)}{L\sqrt{\mathrm{}}})^\lambda `$. By applying the unitary transformation of $`U_f`$ in Eq. (21) to the $`\psi _0`$, one can obtain the wave function $`\psi _a=[_{i=1}^N\mathrm{exp}(iax_i/\mathrm{})]\psi _0`$, up to a purely time-dependent phase. $`\psi _a`$ is an eigenstate of the Hamiltonian. The $`\psi _a`$ has been discussed as an excited state which may be obtained by implementing a Galilei boost to $`\psi _0`$ , while the derivation in this paper clearly supports this interpretation.

###### Acknowledgements.

This work was supported in part by grant No. 2000-1-11200-001-2 from the Basic Research Program of the Korea Science & Engineering Foundation.


text:
# Event-by-event analysis of ultra-relativistic heavy-ion collisions in smoothed particle hydrodynamics

## 1 Introduction

One of the central issues of the high-energy heavy-ion physics is to investigate the phases of the hadronic and quark matter. Hydrodynamical descriptions, which have a rather long history, may be a powerful tool for studying such phase diagrams because we can handle the equation of state directly in the model. However, to extract precisely information about the equation of state, we must be careful in comparing the model predictions with experimental data because procedures of event averaging may make signal of QGP ambiguous.

Event-by-event analysis is one of promising ways to extract clear information from the experimental data, in particular at RHIC and LHC energy regions. The fluctuations at the initial stage (for example, the initial energy density fluctuation) must affect the formation of QGP and the later space-time evolution of the whole system. Fig.1 shows the energy density (counter plot, at $`z=0fm`$ plane) for a Au+Au collision at $`\sqrt{s}=200AGeV`$, $`b=0fm`$, produced by NeXus code. As shown in Fig.1, the fluctuations at the initial stage of relativistic heavy-ion collisions are not negligible. To achieve an event-by-event analysis using relativistic hydrodynamical model, the numerical code must deal with arbitrary initial conditions and equation of state with suitable calculational speed and precision. The method of smoothed particle hydrodynamics (SPH), which was first introduced for astrophysical applications, satisfy such requirements. The main characteristic of SPH is the introduction of “particles” attached to some conserved quantity which are described in terms of discrete Lagrangian coordinates. This feature of the model enables to easily carry out calculations of the ultra-relativistic heavy-ion collisions, which is accompanied by prominent longitudinal expansion. Another advantage of using SPH method is that we can choose a suitable precision in solving hydrodynamical equations because of the smoothing kernel, which allows to smooth out often unnecessary very precise local aspects and then achieve a reduction of computational time. This is a very profitable point in the event-by-event analysis.

In the following section, we briefly formulate the entropy-based SPH, using the variational approach. In Sec.3, to check the performance of SPH, we apply it for solving problems of relativistic hydrodynamics whose solutions are known analytically or numerically. Then, in Sec.4, we demonstrate hydrodynamical evolutions of both resonance gas and quark gluon plasma starting from the same simple initial conditions. In Sec.5, after giving a brief explanation how to take the initial conditions of NeXus into our code, we present typical hydrodynamical evolution of an event at RHIC energy. To introduce freeze-out process into our code, we formulate Cooper-Frye formula in terms of SPH and show several results on single-particle spectrum in Sec.6. Then, we study the event-by-event fluctuation of elliptic flow coefficient $`v_2`$ at RHIC energy in Sec.7, showing that the dispersion is large and it does depend on the equation of state used. Finally we close our discussions with some concluding remarks.

## 2 Relativistic SPH equations

### 2.1 SPH representation

To specify the hydrodynamical state of a fluid, we need to locally know the thermodynamical quantities such as the energy density $`\epsilon `$ and the collective velocity field, besides equation of state of the fluid. If local thermal equilibrium is well established, the thermodynamical variables may be expressed by smooth functions of the coordinates. Let us choose $`N`$ space points to specify the state of the whole system. To know the physical quantities in other space points, we may use an interpolation using the known $`N`$ points. (Hereafter we use the word ‘particles’ instead of ‘points’.)

In SPH, for an arbitrary extensive thermodynamical quantity of $`A`$, the density (in the space-fixed frame) $`a^{}`$ is parametrized in the following way:

$`a^{}(𝒙,t)={\displaystyle \underset{i}{\overset{N}{}}}\nu _iW(𝒙𝒙_i(t);h),`$ (1)
where $`W(𝒙𝒙_i(t);h)`$ is a positive definite kernel function with the properties

$`{\displaystyle 𝑑𝒙^{}W(𝒙^{}𝒙)}1\text{and}W(𝒙^{}𝒙,h)=\delta (𝒙^{}𝒙)\text{when }h0,`$
where $`h`$ is a smoothing-scale parameter. Using the normalization above, we have

$`A_{total}={\displaystyle d^3𝒙a^{}(𝒙,t)}={\displaystyle \underset{i}{\overset{N}{}}}\nu _i.`$ (2)
From eq.(2) the quantity $`\nu _i`$ can be interpreted as a portion of $`A`$ which is carried by the $`i`$th particle. Its velocity, identified as the velocity of the fluid at $`𝒙_i(t)`$, satisfies

$`𝒗_i(t)={\displaystyle \frac{d𝒙_i(t)}{dt}}.`$ (3)
If $`A`$ is a conserved quantity, $`\nu _i`$ must be constant in time. In this case, from eqs.(1) and (3), we obtain

$`{\displaystyle \frac{}{t}}a^{}(𝒙,t)=\mathbf{}𝒋(𝒙,t),`$ (4)
where

$`𝒋(𝒙,t)={\displaystyle \underset{i}{}}\nu _i𝒗_i(t)W(𝒙𝒙_i(t);h)`$ (5)
is interpreted as the current density of the thermodynamical quantity $`A`$ and eq.(4) is then the corresponding continuity equation.

### 2.2 Equations of motion for SPH particles

To know the time evolution of the fluid, we have to find the equation of motion for each SPH particle. It can be obtained by the variational method with the following Lagrangian

$`L_{\text{SPH}}(\{𝒙_i,𝒗_i\})={\displaystyle \underset{i}{}}\nu _i{\displaystyle \frac{\epsilon _i}{a_i^{}}}`$ (6)
It should be noted here that the factor $`\nu _i/a_i^{}`$ can be regarded as the volume $`V_i`$ (in the space-fixed frame) occupied by the $`i`$th ‘particle’. Then each term in eq.(6) is equal to $`E_i/\gamma _i`$, where $`E_i`$ and $`\gamma _i`$ are the energy (in the rest frame) and the gamma factor of the $`i`$th particle, respectively. Then we also write

$`L_{\text{SPH}}(\{𝒙_i,𝒗_i\})={\displaystyle \underset{i}{}}{\displaystyle \frac{E_i(𝒙_i,𝒗_i)}{\gamma _i(𝒗_i)}}.`$ (7)
The equations of motion of SPH particles are obtained minimizing the action $`I_{\text{SPH}}`$ with respect to $`\{𝒙_i\}`$ with the constraint given by eq.(3) taken into account

$`\delta I_{\text{SPH}}(\{𝒙_i,𝒗_i\})=\delta {\displaystyle 𝑑t\frac{E_i}{\gamma _i}}=0.`$ (8)
Using the adiabatic condition

$`\delta E_i=P\delta V_i=P_i\delta \left({\displaystyle \frac{\nu _i\gamma _i}{a_i^{}}}\right),`$ (9)
we obtain

$`{\displaystyle \frac{d}{dt}}(\nu _i{\displaystyle \frac{P_i+\epsilon _i}{a_i}}\gamma _i𝒗_i)+{\displaystyle \underset{j}{}}\nu _j\left[{\displaystyle \frac{P_i}{a_{i}^{}{}_{}{}^{2}}}+{\displaystyle \frac{P_j}{a_{j}^{}{}_{}{}^{2}}}\right]\mathbf{}W(𝒙_i𝒙_j;h)=0,`$ (10)
where $`P_i`$ is the pressure. Because of the adiabatic condition, eqs.(10) are the equations of motion for the ideal fluid. For the extension to the non-adiabatic case, see ref..

### 2.3 Entropy based SPH in general coordinate system

As a possible candidate for the thermodynamical extensive quantity $`a^{}`$, we may consider baryon number density, entropy density, energy density and so on. In the case of energy density, we have (cf. eq.(6) and (7))

$`\nu _i=E_i(𝒙_i,𝒗_i).`$ (11)
However, this is only the internal energy, so $`\nu _i`$ is not constant with respect to the variation of $`𝒙_i`$ in eq.(8), introducing an additional complication without any practical merit. On the other hand, for entropy or baryon number density, for example, $`\nu _i`$ can be kept constant. Since we will apply SPH to ultra-relativistic heavy-ion collisions, baryon number density is not so suitable choice because it is expected to be very small in the central rapidity region. Let us choose the entropy density for $`a^{}`$ in eq.(1).

In relativistic heavy-ion collisions, initial conditions of fluid are often given at constant $`\tau \sqrt{t^2z^2}`$. Then, it is convenient to formulate SPH in the coordinates

$`x^0\tau =\sqrt{t^2z^2},x^1x,x^2y,x^3\eta ={\displaystyle \frac{1}{2}}\mathrm{ln}{\displaystyle \frac{t+z}{tz}}.`$ (12)
The SPH equations in generalized coordinate system are easily derived in a similar way. In the case of the coordinate above, the SPH equations of motion should be written as

$`s(\tau ,𝒙)={\displaystyle \frac{1}{\tau u^0}}{\displaystyle \underset{i}{}}\nu _iW(𝒙𝒙_i:h),`$ (13)
$`{\displaystyle \frac{d}{d\tau }}(\nu _i{\displaystyle \frac{P_i+\epsilon _i}{s_i}}\gamma _i𝒗_i)+{\displaystyle \underset{j}{}}{\displaystyle \frac{\nu _j}{\tau }}\left[{\displaystyle \frac{P_i}{s_{i}^{}{}_{}{}^{2}u_{i}^{0}{}_{}{}^{2}}}+{\displaystyle \frac{P_j}{s_{j}^{}{}_{}{}^{2}u_{j}^{0}{}_{}{}^{2}}}\right]\mathbf{}W[𝒙_i𝒙_j;h]=0.`$ (14)
## 3 Numerical check of SPH

To check the performance of SPH formulated in Sec.2, we applied it to a problem of relativistic hydrodynamics whose solution is known analytically. When the longitudinal rapidity $`\alpha `$ and the transverse one $`\beta `$ are given as

$`\alpha \eta ,\beta {\displaystyle \frac{1}{2}}\mathrm{ln}{\displaystyle \frac{\tau +\sqrt{x^2+y^2}}{\tau \sqrt{x^2+y^2}}},`$ (15)
it is known that the entropy density is given by

$`s={\displaystyle \frac{s_0}{[\tau ^2x^2y^2]^{3/2}}},`$ (16)
where $`s_0`$ is a constant which is determined by the initial conditions. A comparison of the numerical results by SPH with the analytical solution eq.(16) is shown in Fig.2, where the equation of state used was $`c_s^2=1/3`$. It is verified that SPH solves this problem correctly. For a further check, we applied it also to another problem, namely, relativistic hydrodynamics with the longitudinal scaling and the Landau-type initial condition in the transverse directions,

$`\alpha \eta ,\beta 0.`$ (17)
Figure 3 shows a comparison of the SPH results with those obtained by Hama and Pottag. The initial conditions and the equation of state used are exactly the same as in ref.. We found that SPH works well also in this problem.

## 4 Equations of state and hydrodynamical evolutions

To complete the theory of relativistic hydrodynamics, the equation of state (EoS) for the fluid is required. We use two possible EoS in the present study, namely, a resonance gas EoS and an EoS containing both resonance gas and QGP via a first order phase transition. The EoS are parametrized as follows (See also Fig.4.):

$`P/\epsilon =0.20\text{ ( constant for whole }\epsilon \text{ region ),}`$ (18)
for the resonance gas (RG) and

$`P/\epsilon =\{\begin{array}{cc}0.20\hfill & 0.00GeV/fm^3<\epsilon <0.28GeV/fm^3\epsilon _1\hfill \\ 0.056/\epsilon \hfill & 0.28GeV/fm^3<\epsilon <1.45GeV/fm^3\epsilon _2\hfill \\ 1/34B/3\epsilon \hfill & 1.45GeV/fm^3<\epsilon \hfill \end{array}`$ (22)
for the quark gluon plasma (QGP) EoS. The bag constant $`B=0.32GeV/fm^3`$ is determined to give the critical temperature $`T_c=160MeV`$. A remark should be made here that when we call QGP EoS, here and in the following sections, we are actually considering all these three states according to the value of $`\epsilon `$.

We show some results of the hydrodynamical evolution using such RG and QGP EoS in Figs. 5 and 6. For the both cases, the initial ($`\tau _0=1.0fm`$) entropy density $`s`$, longitudinal and transverse rapidity $`\alpha `$, $`\beta `$ are set as

$`s=40[1/fm^3],\alpha =\eta ,\beta =0.\text{for }\sqrt{x^2+y^2}<7.0fm\text{}|\eta |<7.0\text{ .}`$ (23)
As seen in Figs.5 and 6, the space-time evolution of the fluid strongly depends on the EoS. In comparison with the RG case, the temperature of the fluid in QGP EoS drops much slowly due to the existence of the mixed phase, where the pressure is constant.

## 5 Initial conditions produced by NeXus

The most fundamental assumption of hydrodynamical models is the thermalization of partonic or hadronic systems produced after high-energy heavy-ion collisions. In our study, we set the initial conditions using information from Nexus event generator, in the following 2 steps.

1. Align all primary hadrons on the $`\tau =\tau _0`$ hypersurface. Here we assume that the primary hadrons move freely with the momentum $`p^\mu `$. In the present study, we use $`\tau _0=1.0fm`$.
2. The initial energy density and momentum density are estimated using the interpolation kernel $`W`$.

$`\epsilon (\tau _0;𝒙)={\displaystyle \underset{i}{}}p_i^0W(\tau _0,𝒙𝒙_i;h),`$
$`𝒑(\tau _0;𝒙)={\displaystyle \underset{i}{}}𝒑_iW(\tau _0,𝒙𝒙_i;h),`$
where $`𝒙_i=(x_i,y_i,\eta _i)`$ are the coordinates of the $`i`$-th primary hadron after the alignment on $`\tau =\tau _0`$. The smoothing scale parameter $`h=1fm`$ is used, which is the rough size of the hadrons. The collective velocity field is then given by $`𝒗(\tau _0,𝒙)=𝒑(\tau _0,𝒙)/\epsilon (\tau _0,𝒙)`$ .

In Fig.7, we show the hydrodynamical evolution of a typical event of $`Au+Au`$ collision at energy $`\sqrt{s}=130AGeV`$ for the impact parameter $`b=7.0fm`$ given by our hydro code, SPheRIO111 Smoothed Particle hydrodynamical evolution of Relativistic heavy IOn collisions, connected to the NeXus event generator.

## 6 Cooper-Frye formula in SPheRIO

The Cooper-Frye formula is widely used in calculations of single particle spectra using hydrodynamical models because of its simplicity. So, as the first trial to evaluating single particle densities, we shall follow it.

In SPH, using eq.(13), the usual Cooper-Frye formula can be rewritten as

$`E{\displaystyle \frac{d^3N}{d𝒑^3}}`$ $`=`$ $`{\displaystyle \frac{g}{(2\pi )^3}}{\displaystyle \underset{j}{}}\nu _j{\displaystyle _{F.O.}}\tau 𝑑x𝑑y𝑑\eta {\displaystyle \frac{\sigma _\mu p^\mu }{\tau su^0}}W(𝒙𝒙_j(\tau );h){\displaystyle \frac{1}{e^{u_\mu p^\mu /T_f}1}},`$ (24)
where $`\sigma _\mu `$ is defined by

$`\sigma _\mu (1,{\displaystyle \frac{n_x}{n_\tau }},{\displaystyle \frac{n_y}{n_\tau }},{\displaystyle \frac{n_\eta }{n_\tau }})`$ (25)
and $`n_\mu =(n_\tau ,n_x,n_y,n_\eta )`$ is the normal vector of the isothermals:

$`n_\mu ({\displaystyle \frac{T}{\tau }},{\displaystyle \frac{T}{x}},{\displaystyle \frac{T}{y}},{\displaystyle \frac{T}{\eta }}).`$ (26)
The integration is usually done over a constant temperature hypersurface, $`T=T_f`$. When the spline-type kernel function is used, $`W(𝒙𝒙^{};h)=0`$ for $`|𝒙𝒙^{}|>2h`$. Then, within the approximation of small $`h`$,

$`E{\displaystyle \frac{d^3N}{d𝒑^3}}`$ $``$ $`{\displaystyle \frac{g}{(2\pi )^3}}{\displaystyle \underset{j}{}}\nu _j\left[{\displaystyle \frac{p^\mu }{e^{u_\mu p^\mu /T_f}1}}\right]_j{\displaystyle \tau 𝑑x𝑑y𝑑\eta \frac{\sigma _\mu }{\tau su^0}W(𝒙𝒙_j(\tau _\text{F}(𝒙));h)}`$ (27)
$`=`$ $`{\displaystyle \frac{g}{(2\pi )^3}}{\displaystyle \underset{j}{}}\left[{\displaystyle \frac{p^\mu }{e^{u_\mu p^\mu /T_f}1}}\right]_j{\displaystyle 𝑑x𝑑y𝑑\eta \frac{\nu _j}{su^0}\frac{\sigma _\mu W(𝒙𝒙_j;h)}{|1+v_j^x\frac{n_x}{n_\tau }+v_j^y\frac{n_y}{n_\tau }+v_j^\eta \frac{n_\eta }{n_\tau }|}},`$
where $`𝒙_j`$ gives the point where $`j`$-th SPH-particle crosses the freeze-out hypersurface. The result of the single-particle spectra, $`1/m_\text{T}dN/dm_\text{T}`$, $`dN/dy`$ and $`dN/d\phi `$ for Au+Au collisions at $`\sqrt{s}=130AGeV`$, $`b=7.0fm`$ are shown in Figs. 8, 9 and 10, respectively. The freeze-out temperature used is $`T_f=140MeV`$.

## 7 Event-by-event analysis of hadronic observables

### 7.1 Elliptic flow coefficient $`v_2`$

The flow phenomena can be important candidates of QGP signals because they may carry much information about EoS during the expansion of quark or hadronic matter. In particular, the elliptic flow coefficients is one of the interesting observables since it is sensitive to the early stage of high-energy heavy-ion collisions. Because it is expected that different EoS responds to the initial fluctuations in different way, the measure of fluctuation, for example, $`\delta v_2\sqrt{v_2^2v_2^2}`$ is also very important. For this purpose, we investigate $`v_2`$ distribution using different EoS. In Figs.11, 12 and Table 1, we show results of the $`v_2`$ distribution for $`Au+Au`$ collisions at energy $`\sqrt{s}=130AGeV`$ in RG and QGP EoS cases.

As seen in Table 1, the fluctuations of flow coefficient $`v_2`$ in QGP is about 20 $``$ 30% smaller than that in RG (cf. the difference of $`v_2`$ in RG and QGP. It is about a few %.). The $`b`$ dependence of $`\delta v_2`$ is also found. Hence, it is expected that $`\delta v_2`$ also brings us useful information to infer on the phase passed in the early stage of high-energy nuclear collisions.

### 7.2 Fluctuation of the slope parameter

It is also interesting to study the fluctuation of the so-called slope parameter in the transverse momentum spectra. The slope parameter $`\stackrel{~}{T}`$ can be obtained by fitting $`1/m_\text{T}dN/dm_\text{T}`$ with a function $`𝒩\mathrm{exp}[m_\text{T}/\stackrel{~}{T}]`$. In the present study, $`0<m_\text{T}m_\pi <1.0GeV`$ interval is used for the fitting. The results are shown in Table 2. We found $`\stackrel{~}{T}`$ and $`\delta \stackrel{~}{T}`$ of QGP is systematically smaller than that of RG.

### 7.3 Multiplicity fluctuation in the central region

Other interesting measurement may be multiplicity fluctuation. Let us define the fluctuation measure as

$`\delta n(y,\mathrm{\Delta }y)\sqrt{n^2(y,\mathrm{\Delta }y)n(y,\mathrm{\Delta }y)^2}`$ (28)
where

$`n(y,\mathrm{\Delta }y)={\displaystyle _{y\mathrm{\Delta }y}^{y+\mathrm{\Delta }y}}𝑑y{\displaystyle \frac{dN}{dy}}.`$ (29)
We investigate the EoS dependence of the multiplicity fluctuation in two central ($`y=0`$) rapidity windows $`\mathrm{\Delta }y`$ and as function of the impact parameter $`b`$. The results are summarized in Table 3. As seen there, the multiplicities and their fluctuations $`\delta n`$ have a clear EoS dependence. See Figs. 13 and 14. When the initial energy density is bellow $`\epsilon _2`$ the multiplicity fluctuation with QGP EoS is expected to be larger than that of RG without first order phase transition. This is due to the existence of the mixed phase (see Appendix). The values of ratio $`\delta n/n`$ of QGP is slightly larger than those of RG.

## 8 Concluding remarks and discussions

We developed a new hydrodynamical code, SPheRIO, based on the smoothed particle hydrodynamics (SPH), for studying high-energy heavy-ion collisions in event-by-event basis. In this study, we set the initial conditions based on the results of NeXus event generator and used two types of possible equation of states (EoS), $`ie.`$, the resonance gas (RG) and quark gluon plasma (QGP) EoS. In the freeze-out process, single particle spectra are calculated by the Cooper-Frye formula for each set of initial conditions and EoS. Analyzing these single particle spectra, the fluctuations in flow coefficients, slope parameter $`\stackrel{~}{T}`$ and multiplicities in the central rapidity regions are investigated at RICH energy.

As seen in the results of Table 1, we found remarkable difference in $`\delta v_2`$ between RG and QGP EoS predictions. This means that the measure of the fluctuations is also helpful to infer the phase passed in the high-energy heavy-ion collisions. In the QGP EoS case, the fluctuations are about $`2030`$ % smaller than those in RG case. Our predictions of the average value of $`v_2`$ is shown in Fig.15, compared with STAR Collaboration data. It is seen that the main trend of the data is reproduced, but the fluctuation is quite large.

The fluctuation of the slope parameter $`\stackrel{~}{T}`$ is small, a few % of $`\stackrel{~}{T}`$. We found the fluctuation of multiplicity has a clear EoS dependence. The impact parameter dependence of $`\delta n(y,\mathrm{\Delta }y)`$ can be a good measure for the detection of the first order phase transition.

Acknowledgments: The authors acknowledge stimulating discussions with K. Werner, O. Socolwski Jr and H.J. Drescher. This work was partially supported by FAPESP(contract no.s 98/02249-4 and 98/00317-2), PRONEX(contract no. 41.96.9886.00), FAPERJ(contract no. E-26/150.942/99) and CNPq.

## Appendix A Entropy fluctuations in the first order phase transition

### A.1 Mixed phase dominant case:

If in the most part of the fluid the energy density at the initial time $`\tau _0`$ varies only in the interval $`\epsilon _1<\overline{\epsilon }(x)<\epsilon _2`$, where $`\epsilon _1`$ and $`\epsilon _2`$ are the critical densities of the RG and QGP phases, respectively, the produced entropy $`S`$ of the system is

$`S`$ $`=`$ $`{\displaystyle d^3xs_1\left[\frac{\overline{\epsilon }(x)+\delta \epsilon (x)}{\epsilon _1}\right]^{1/[1+c_m^2]}}{\displaystyle d^3xs_1\left(\frac{\overline{\epsilon }(x)}{\epsilon _1}\right)\left[1+\left(\frac{\delta \epsilon (x)}{\overline{\epsilon }(x)}\right)\right]}`$
where $`c_m=0`$(constant) is the sound velocity in the mixed phase(Mix) and $`s_1`$ is the critical entropy density corresponding to $`\epsilon _1`$. $`\overline{\epsilon }(x)`$ and $`\delta \epsilon (x)`$ denote the event averaged energy density and fluctuation, respectively. Here, we assume that $`\delta \epsilon (x)\overline{\epsilon }(x)`$, so that the mixed phase is dominant in all the events. Then, the dispersion of the total entropy $`S`$ is

$`D_{Mix}^2=S^2S^2{\displaystyle d^3xs_1^2\left(\frac{\overline{\epsilon }(x)}{\epsilon _1}\right)^2\left(\frac{\delta \epsilon (x)}{\overline{\epsilon }(x)}\right)^2}`$ (A2)
In the case the phase transition is absent, so that the fluid is described by RG EoS, the same events with the same energy distribution and fluctuation gives

$`S_{RG}`$ $`=`$ $`{\displaystyle d^3xs_1\left[\frac{\overline{\epsilon }(x)+\delta \epsilon (x)}{\epsilon _1}\right]^{1/[1+c_r^2]}}`$ (A3)
$``$ $`{\displaystyle d^3xs_1\left(\frac{\overline{\epsilon }(x)}{\epsilon _1}\right)^{5/6}\left[1+\frac{5}{6}\left(\frac{\delta \epsilon (x)}{\overline{\epsilon }(x)}\right)\frac{5}{72}\left(\frac{\delta \epsilon (x)}{\overline{\epsilon }(x)}\right)^2+\mathrm{}\right]}`$
and

$`D_{RG}^2={\displaystyle d^3xs_1^2\left(\frac{\overline{\epsilon }(x)}{\epsilon _1}\right)^{5/3}\frac{25}{36}\left(\frac{\delta \epsilon (x)}{\overline{\epsilon }(x)}\right)^2},`$ (A4)
respectively. We used the sound velocity $`c_r=1/\sqrt{5}`$(constant) here. So, the appearance of the mixed phase produces

$`\mathrm{\Delta }D^2`$ $``$ $`D_{Mix}^2D_{RG}^2={\displaystyle d^3xs_1^2\left\{\left(\frac{\overline{\epsilon }(x)}{\epsilon _1}\right)^2\frac{25}{36}\left(\frac{\overline{\epsilon }(x)}{\epsilon _1}\right)^{5/3}\right\}\left(\frac{\delta \epsilon (x)}{\overline{\epsilon }(x)}\right)^2}>0`$
Hence,

$`D_{Mix}^2>D_{RG}^2.`$ (A5)
### A.2 QGP phase dominant case:

If the most part of the fluid is in the QGP phase at initial time $`\tau _0`$, the total produced entropy $`S`$ of the system is

$`S_{QGP}`$ $`=`$ $`{\displaystyle d^3xs_2\left[\frac{\overline{\epsilon }(x)+\delta \epsilon (x)}{\epsilon _2}\right]^{1/[1+c_q^2]}}`$ (A6)
$``$ $`{\displaystyle d^3xs_1\left(\frac{\epsilon _2}{\epsilon _1}\right)\left(\frac{\overline{\epsilon }(x)}{\epsilon _2}\right)^{3/4}\left[1+\frac{3}{4}\left(\frac{\delta \epsilon (x)}{\overline{\epsilon }(x)}\right)\frac{3}{32}\left(\frac{\delta \epsilon (x)}{\overline{\epsilon }(x)}\right)^2+\mathrm{}\right]}`$
and

$`D_{QGP}^2={\displaystyle d^3xs_1^2\left(\frac{\epsilon _2}{\epsilon _1}\right)^2\left(\frac{\overline{\epsilon }(x)}{\epsilon _2}\right)^{3/2}\frac{9}{16}\left(\frac{\delta \epsilon (x)}{\overline{\epsilon }(x)}\right)^2},`$ (A7)
where $`c_q=1/\sqrt{3}`$(constant) is assumed. Then, we have

$`\mathrm{\Delta }D^2`$ $``$ $`D_{QGP}^2D_{RG}^2={\displaystyle d^3xs_1^2\left\{\frac{9}{16}\left(\frac{\epsilon _2}{\epsilon _1}\right)^2\left(\frac{\overline{\epsilon }(x)}{\epsilon _2}\right)^{3/2}\frac{25}{36}\left(\frac{\overline{\epsilon }(x)}{\epsilon _1}\right)^{5/3}\right\}\left(\frac{\delta \epsilon (x)}{\overline{\epsilon }(x)}\right)^2}.`$
From this relation, we can conclude that when the energy density $`\overline{\epsilon }(x)`$ becomes larger than a certain value, $`(0.81)^6(\epsilon _2/\epsilon _1)^3\epsilon _111GeV/fm^3`$ for the most part of the fluid, the dispersion of the QGP case becomes smaller than RG case.


text:
# 1 Introduction

## 1 Introduction

The relationship between the integrable systems and the Lie group gauge symmetries is well known ,,,. However, Lie groups by no means exhaust all the symmetries and in fact there are a lot of systems whose gauge symmetries cannot be reduced to the Lie group action  . For this reason it is quite interesting to construct an integrable system using the non Lie group symmetries. The aim of the paper is to present the example of such integrable system. We obtain the system via the symplectic reduction from an extended phase space, where the symmetries are realized by the action of the quadratic Sklyanin algebra . Let us briefly explain the key idea for realizing this action. Note that the well known coadjoint action of a Lie algebra on its dual space can be written in terms of the linear Poisson structure, naturally associated to the Lie algebra. Namely if $`p_\mu `$ are coordinates in the dual space and the respective Poisson tensor is of the form

$$\alpha _{\mu \nu }(p)=f_{\mu \nu }^\lambda p_\lambda ,$$

where $`f_{\mu \nu }^\lambda `$ are the structure constants of the Lie algebra, then the coadjoint action is rewritten as

$$\delta _ϵp_\mu =\alpha _{\mu \nu }(p)ϵ^\nu (p),$$
(1)
where $`ϵ^\nu `$ are the infinitesimal parameters of the action. If we assume now that the Poisson tensor $`\alpha _{\mu \nu }(p)`$ in (1) is no longer linear in $`p_\mu `$ we get the generalization of the coadjoint action to the case of an arbitrary Poisson manifold. Such a generalization of the coadjoint action has been originally proposed by Karasev  and in general it cannot be reduced to the Lie algebra action in a sense that one cannot choose the global basis of the sections $`ϵ^\mu (p)`$ with a Lie algebra commutation relations of the transformations (1). Instead the form of the commutation relations in a general case looks as follows

$$[\delta _\xi ,\delta _\eta ]=\delta _ϵ,$$

where

$$\begin{array}{c}ϵ=[\xi ,\eta ]=d\alpha (\xi ,\eta )+\alpha (d\xi ,\eta )+\alpha (\xi ,d\eta ),\\ \\ ϵ^\mu =^\mu (\alpha _{\nu \lambda }\xi ^\nu \eta ^\lambda )+\alpha _{\nu \lambda }(^\nu \xi ^\mu ^\mu \xi ^\nu )\eta ^\lambda +\alpha _{\nu \lambda }\xi ^\nu (^\lambda \eta ^\mu ^\mu \eta ^\lambda ),\end{array}$$
(2)
$$^\mu =\frac{}{p_\mu }.$$

Thus we get an example of the manifold, equipped with the non Lie algebra symmetries. A mathematical rigor treatment of such symmetries involves a definition of the so-called Lie algebroid ,, which is defined as a bundle $`𝒜`$ over a manifold $`M`$ with a Lie bracket, given in the set of sections $`\mathrm{\Gamma }(𝒜)`$ and an anchor map $`\delta :\mathrm{\Gamma }(𝒜)\mathrm{\Gamma }(TM)`$ , being a homomorphism of the respective Lie algebras. Using this terminology we see that the above construction is an example of the Lie algebroid. Namely, the Lie algebroid is the cotangent bundle over the Poisson manifold. The Lie bracket between the sections $`\xi `$ and $`\eta `$ is given by the equation (2) and the anchor map is defined by the formula (1)aaasuch example of the Lie algebroid, naturally associated to a Poisson manifold has been originally presented in the paper .. Note that the Jacobi identity of the Lie bracket (2) holds in virtue of the respective Jacobi identity for the Poisson tensor $`\alpha _{\mu \nu }(p)`$ . In this paper we start from the Poisson manifold $`𝐑^4`$, equipped with the standard classical Sklyanin bracket  , and thus equipped with the algebroid action, which may be regarded as the action of the Sklyanin algebra. We choose the extended phase space of the integrable system to be the cotangent bundle over $`𝐑^4`$ with the standard symplectic structure. We lift the algebroid transformations from the base manifold to the canonical ones in the cotangent bundle and readily find the respective Hamiltonian generators, which define a natural analogue of the momentum map for the case of the non Lie group symmetries. Then we construct the phase space of the integrable system by performing the symplectic reduction on the surface of the non zero level of the generators. The underlying Hamiltonians are obtained by reduction of the Casimir functions of the Sklyanin bracket. Note that since the Sklyanin bracket  is the deformation of the Poisson bracket, which is equivalent to the linear bracket, associated to the Lie algebra $`𝐮_\mathrm{𝟐}`$ the presented symplectic reduction is similar to the one leading to the rational Calogero-Moser system ,. In particular, one of the exponential degenerations of the Sklyanin bracket leads to the system, whose dynamics is in a sense equivalent to the dynamics of the two-particle rational Calogero-Moser system. As an aside, it is worth noting that although the integrable system is obtained via the symplectic reduction we still cannot find the respective Lax representation because the symmetries, that have been used for constructing the system are not Lie group ones. However, we show that it is possible to resolve the equations of motion for our system via the projection method. So the problem of finding the respective Lax representation is seemed to be rather intriguing and we are going to comment on it in the concluding section. The organization of the paper is as follows. In the second section we present the construction of the integrable system. In the third section we develop the projection method and explicitly resolve the equations of motion. The fourth section is devoted to the quantization. We find that although it is not possible to pose the bound states problem for our system one of the Hamiltonians possesses a discrete spectrum. At last, in the concluding section we discuss some open questions. All through the paper the Greek indices take values from $`0`$ to $`3`$ the Latin indices run from $`1`$ to $`3`$ and whenever the indices $`i,j,k`$ are meat they are assumed to take values $`1,2,3`$ together with their cyclic permutations. $`^\mu `$ denotes $`{\displaystyle \frac{}{p_\mu }}`$ and $`_\nu `$ is reserved for $`{\displaystyle \frac{}{x^\nu }}`$ .

## 2 Construction of the System

We construct the system starting with the extended phase space, which is chosen to be the cotangent bundle $`T^{}𝐑^4`$ , equipped with the canonical symplectic structure

$$\omega =dp_\mu dx^\mu .$$
(3)
Here we use the following unusual conventions. The coordinates on the base $`𝐑^4`$ are denoted by $`p_\mu `$ and the coordinates in the cotangent space are denoted by $`x^\mu `$ . Let us endow the base manifold $`𝐑^4`$ with the quadratic Sklyanin bracket

$$\{p_0,p_i\}=\alpha _{0i}(p)=J_{jk}p_jp_k\{p_i,p_j\}=\alpha _{ij}(p)=p_0p_k,$$
(4)
where the coefficients $`J_{12}`$, $`J_{23}`$ and $`J_{31}`$ satisfy the identity

$$J_{12}+J_{23}+J_{31}=0,$$
(5)
which is equivalent to the Jacobi identity for the bracket (4). The condition (5) also implies that these coefficients can be represented in the form

$$J_{ij}=J_iJ_j,$$
(6)
where the values $`J_i`$ are defined up to the additional constant $`c`$, $`J_iJ_i+c`$ . Following the steps, outlined in the introduction we firstly define the analogue of the coadjoint transformations on the base manifold $`𝐑^4`$ using the Poisson structure (4). Namely, the infinitesimal coadjoint transformations are defined as

$$\delta _ϵp_\mu =\alpha _{\mu \nu }(p)ϵ^\nu (p),$$
(7)
where $`ϵ^\nu `$ are the respective infinitesimal parameters. The transformations (7) can be easily lifted to the canonical ones in the cotangent bundle

$$\delta _ϵp_\mu =\alpha _{\mu \nu }ϵ^\nu (p)\delta _ϵx^\mu =^\mu (\alpha _{\nu \lambda }ϵ^\nu (p)x^\lambda )$$
(8)
and the Hamiltonian generators corresponding to the transformations (8) can be also readily foundbbbthus we have constructed the so-called Hamiltonian algebroid  over $`T^{}𝐑^4`$ ; in the paper  it is shown that the natural construction of the Hamiltonian algebroid exists for an arbitrary Lie algebroid

$$M_ϵ=\alpha (ϵ,x)=\alpha _{\mu \nu }(p)ϵ^\mu x^\nu .$$
(9)
Then in order to obtain a non-trivial dynamical system via the symplectic reduction method we consider the surface of non-zero level of the generators (9) in the extended phase space $`T^{}𝐑^4`$. In analogy to what it is done for case of the rational Calogero-Moser system we set

$$\begin{array}{c}M_0=J_{23}p_2p_3x^1+J_{31}p_3p_1x^2+J_{12}p_1p_2x^3=0\\ \\ M_1=p_0(x^2p_3x^3p_2)J_{23}p_2p_3x^0=\nu \\ \\ M_2=p_0(x^3p_1x^1p_3)J_{31}p_3p_1x^0=0\\ \\ M_3=p_0(x^1p_2x^2p_1)J_{12}p_1p_2x^0=0.\end{array}$$
(10)
The symplectic form (3) being reduced on the constraint surface (10) is of course degenerate and the respective kernel distribution is given by the following vector field, being tangent to the surface (10)

$$\begin{array}{cc}\delta p_0=J_{23}p_2p_3& \delta x^0=(x^2p_3x^3p_2)\\ & \\ \delta p_1=0& \delta x^1=0\\ & \\ \delta p_2=p_0p_3& \delta x^2=(p_0x^3+J_{23}p_3x^0)ϵ^1\\ & \\ \delta p_3=p_0p_2& \delta x^3=(p_0x^2J_{23}p_2x^0).\end{array}$$
(11)
Our system is then presented as a result of the symplectic reduction on the constraint surface (10) with respect to the transformations, corresponding to the vector field (11). In what follows, these transformations are called the gauge ones and the integral curves of the vector field (11) are called the gauge orbits. We perform the symplectic reduction using the gauge fixing condition $`x^2=0`$. In fact it is easy to see that this gauge condition has a horizon, given by the equation

$$p_0x^3J_{32}p_3x^0=0.$$
(12)
Namely, for the points, satisfying (12) the coordinate $`x^2`$ is unchange under the gauge transformations (11), and hence for these points the value of the coordinate $`x^2`$ in general cannot be reduced to zero with the help of the gauge transformations. However, in the paper we do not bother the subtlety and we are not going to discuss it. So taking into account the condition $`x^2=0`$ together with the constraint equations (10) we get the embedding of the reduced phase space into the extended one. Namely, the embedding is given by the following mapping

$$\begin{array}{cccc}x^0=u^1& x^1=0& x^2=0& x^3=u^2\\ p_0=v_1& p_1=0& p_2=\frac{\nu }{v_1u^2J_{32}v_2u^1}& p_3=v_2,\end{array}$$
(13)
where $`u^1,v_1,u^2,v_2`$ are the canonical coordinates parameterizing the reduced phase space. The pair of the Casimir functions of the Sklyanin bracket

$$C_1=\frac{1}{2}\underset{l}{}(p_l)^2C_2=\frac{1}{2}((p_0)^2+\underset{l}{}J_l(p_l)^2)$$
(14)
provides us with the required set of the commuting and independent Hamiltonians of the presented system. Really, the gauge invariance of the Casimir functions and their involution property with respect to the canonical structure on $`T^{}𝐑^4`$ ensures that the Hamiltonians, obtained by reduction of (14) on the surface of the constraint equations and of the gauge fixing condition are commuting and independent. Thus we readily find the Hamiltonians for the our system

$$H_1=\frac{(v_2)^2}{2}+\frac{\nu ^2}{2(v_1u^2J_{32}u^1v_2)^2}H_2=\frac{(v_1)^2}{2}+J_{32}\frac{(v_2)^2}{2},$$
(15)
where we set $`J_2=0`$ using the arbitrariness of the coefficients $`J_i`$ . We have presented the example of the integrable system, obtained via the symplectic reduction, which is performed in the context of the non Lie group symmetries. The most similar symplectic reduction performed in the group theoretical context is the one leading to the two-particle rational Calogero-Moser system ,. This resemblance of the reductions becomes quite clear if we compare the Poisson manifold $`𝐑^4`$ (4) to the dual space of the Lie algebra $`𝐮_\mathrm{𝟐}`$ and the coadjoint transformations (7) to the coadjoint action of the Lie algebra $`𝐮_\mathrm{𝟐}`$ . Note that the exponential degeneration of the Sklyanin bracket corresponding to case $`J_{23}=0`$ leads us to the system with the following Hamiltonians

$$H_1=\frac{(v_2)^2}{2}+\frac{\nu ^2}{2(v_1u^2)^2}H_2=\frac{(v_1)^2}{2},$$
(16)
which describe the dynamics of the two-particle Calogero-Moser system with the coupling constant, depending on the initial value of $`v_1`$ and with $`u^1`$ , being the coordinate of the mass center.

## 3 Projection Method

In the this section we show that in spite of the more complicated nature of the gauge transformations the equations of motion for the presented system (15) can be explicitly resolved via the projection method. Let us briefly recall the main idea of the method. It is well known that each point of the reduced phase space corresponds to a gauge orbit on the constraint surface, where the dynamics is described by the gauge invariant Hamiltonians ((14) in our case). We integrate the equations of motion defined by the Hamiltonians (14) and construct the surface, which consists of the gauge orbits, passing through to all the points of the evolution curves. The intersections of the surface with the level of the gauge fixing condition are just the evolution curves of the reduced system. Thus the explicit solutions of the evolution equations for the reduced system can be derived if the explicit description of the gauge orbits is known. In order to describe the gauge orbits we have to integrate the system of differential equations given by the kernel vector field (11)

$$\begin{array}{cc}\frac{d}{ds}p_0=J_{23}p_2p_3& \frac{d}{ds}x^0=(x^2p_3x^3p_2)\\ & \\ \frac{d}{ds}p_1=0& \frac{d}{ds}x^1=0\\ & \\ \frac{d}{ds}p_2=p_0p_3& \frac{d}{ds}x^2=(p_0x^3+J_{23}p_3x^0)\\ & \\ \frac{d}{ds}p_3=p_0p_2& \frac{d}{ds}x^3=(p_0x^2J_{23}p_2x^0),\end{array}$$
(17)
where $`s`$ parameterizes the gauge orbit. Note that the r.h.s. of the differential equations for the momenta contains the momenta only. Hence we can integrate these equations separately. Then using the fact that the gauge transformations (11) are canonical ones we can resolve the rest equations for the coordinates with the help of the following formula

$$x^\mu (s)=\frac{p_\nu ^0}{p_\mu }x_0^\nu ,$$
(18)
where $`x_0^\mu ,p_\nu ^0`$ denote the initial conditions for the system (17). Note also that the phase coordinates $`x^1,p_1`$ are not changed under the gauge transformations so in what follows we are going to omit them. It is easy to see that the system of differential equation for the momenta is similar to the system of equations for the ordinary $`\mathrm{𝐬𝐨}_\mathrm{𝟑}`$ top . Using the analogy we can readily find two independent integrals of motion for the system (17)

$$A=\sqrt{(p_2)^2+(p_3)^2}B=\sqrt{(p_0)^2+J(p_3)^2},$$
(19)
where $`J=J_{32}`$ (in the rest of the paper we assume that $`J>0`$). Let us consider the integrals (19) as the new independent momenta. Then the rest third coordinate in the momentum space can be chosen in two different ways. First, one can introduce the following angle variable

$$\phi =arctg(p_3/p_2),$$
(20)
for which the momenta $`p_0,p_2,p_3`$ take the following form

$$p_0=\sqrt{B^2JA^2sin^2\phi }p_2=Acos\phi p_3=Asin\phi $$
(21)
Second, the third coordinate can be chosen as

$$\phi _2=arctg(\sqrt{J}p_3/p_0).$$
(22)
And in this case we get

$$p_0=Bcos\phi _2p_2=\sqrt{A^2\frac{B^2}{J}sin^2\phi _2}p_3=\frac{B}{\sqrt{J}}sin\phi _2$$
(23)
The first way to choose the new momenta is applicable when

$$k=\frac{\sqrt{J}A}{B}<1,$$
(24)
because it is the case when the expression $`B^2JA^2sin^2\phi `$ do not change the sign. Otherwise, one should use the second way of choosing coordinates in the momentum space. In both cases one can develop the projection method and find the evolution equations of the reduced system proceeding similarly at each stage. So in our paper we are going to discuss the projection method only for the first case (24) omitting another possibility. In what follows, it is more convenient to use instead of $`x^0,x^2,x^3`$ the coordinates, which are canonically conjugated to $`A,B`$ and $`\phi `$. They can be readily found in the following form

$$\begin{array}{c}x_A=x^2cos\phi +x^3sin\phi \frac{JAsin^2\phi }{\sqrt{B^2JA^2sin^2\phi }}x^0\hfill \\ x_B=\frac{Bx^0}{\sqrt{B^2JA^2sin^2\phi }}\hfill \\ x_\phi =Ax^2sin\phi +Ax^3cos\phi \frac{JA^2sin\phi cos\phi }{\sqrt{B^2JA^2sin^2\phi }}x^0.\hfill \end{array}$$
(25)
In terms of $`A,B,\phi ,x_A,x_B`$ and $`x_\phi `$ the solution for the system (17) looks as follows

$$\begin{array}{c}\phi (s)=am(Bs+F(\phi _0;k);k)A(s)=A_0B(s)=B_0\\ x_\phi (s)=\sqrt{\frac{1k^2sin^2\phi _0}{1k^2sin^2\phi }}x_\phi ^0\\ x_A(s)=x_A^0+\frac{\sqrt{J}}{B}\sqrt{1k^2sin^2\phi _0}(F_k(\phi ;k)F_k(\phi _0;k))x_\phi ^0\\ x_B(s)=x_B^0\sqrt{1k^2sin^2\phi _0}(s+\frac{A\sqrt{J}}{B^2}[F_k(\phi ;k)F_k(\phi _0;k)])x_\phi ^0,\end{array}$$
(26)
where $`u=am(t;k)`$ is the elliptic Jacobi function being the inverse function to the first order elliptic integral

$$F(u;k)=_0^u\frac{dt}{\sqrt{1k^2sin^2t}},$$

$`F_k(u;k)`$ denotes the corresponding derivative

$$F_k(u;k)=_0^u\frac{ksin^2tdt}{(1k^2sin^2t)^{3/2}}$$

and $`A_0,B_0,\phi _0,x_A^0,x_B^0x_\phi ^0`$ are the initial conditions of the respective phase coordinates. Now suppose we are given the initial conditions $`\overline{u}^1,\overline{v}_1,\overline{u}^2,\overline{v}_2`$ for the equations of motion of the reduced system (15)

$$\begin{array}{cc}_{t_1}u^1=\frac{\nu ^2u^2}{(v_1u^2Jv_2u^1)^3}& _{t_2}u^1=v_1\\ _{t_1}u^2=v_2+\frac{J\nu ^2u^1}{(v_1u^2Jv_2u^1)^3}& _{t_2}u^2=Jv_2\\ _{t_1}v_1=\frac{J\nu ^2v_2}{(v_1u^2Jv_2u^1)^3}& _{t_2}v_1=0\\ _{t_1}v_2=\frac{\nu ^2v_1}{(v_1u^2Jv_2u^1)^3}& _{t_2}v_2=0.\end{array}$$
(27)
The initial point of the evolution curve (27) is mapped via the embedding (13) to the point in the extended phase space with the phase coordinates $`\overline{x}^\mu ,\overline{p}_\nu `$ , such that $`\overline{x}^1=\overline{p}_1=0`$ and the gauge fixing condition $`\overline{x}^2=0`$ is satisfied. Here we assume for simplicity that the inequality (24) is also satisfied. The gauge fixing condition and the gauge invariant Hamiltonians (14) in terms of the phase coordinates (19),(20),(25) take the following form

$$Ax_Acos\phi =x_\phi sin\phi ,$$
(28)
$$C_1=\frac{(p_1)^2}{2}+\frac{A^2}{2}C_2=J_1\frac{(p_1)^2}{2}+\frac{B^2}{2}.$$
(29)
Solving the equations of motion with respect to the Hamiltonians (29) we get that the momenta (19,20) are conserved $`x^1`$ and $`p_1`$ do not change their zero values and the coordinates (25) are linear functions of the times $`t_1`$ and $`t_2`$, corresponding to the Hamiltonians $`C_1`$ and $`C_2`$ (29) respectively, namely

$$\begin{array}{cc}x_0^1=0& p_1^0=0\\ x_\phi ^0=\overline{x}_\phi & \phi _0=\overline{\phi }\\ x_A^0=\overline{x}_A+\overline{A}t_1& A_0=\overline{A}\\ x_B^0=\overline{x}_B+\overline{B}t_2& B_0=\overline{B}.\end{array}$$
(30)
Here we suspend the phase space coordinates by $`0`$ and in order to show that they should be considered as the initial conditions for the system of differential equations (17) because in general a point of the phase curve (30) does not satisfy the gauge fixing condition (28) but just define the orbit that should intersect the level of the gauge fixing condition in another point. Substituting now the phase space coordinates of the curve (26) with the initial data (30) into the gauge fixing condition (28) we get the following equation

$$\overline{A}(\overline{x}_A+\overline{A}t_1)=\overline{x}_\phi (\sqrt{\frac{1k^2sin^2\overline{\phi }}{1k^2sin^2\phi (s)}}tg\phi (s)$$

$$k\sqrt{1k^2sin^2\overline{\phi }}(F_k(\phi (s);k)F_k(\overline{\phi };k))),$$
(31)
where $`k={\displaystyle \frac{\sqrt{J}\overline{A}}{\overline{B}}}`$ . Since the r.h.s. of the equation (31) does not explicitly depend on the parameter $`s`$ we can forget about it and treat (31) as the equation, defining $`\phi `$ as an implicit function of $`t_1`$. Really if we extract the dependence on $`\phi `$ in the r.h.s. of (31) we get the function

$$f(\phi )=\frac{tg\phi }{\sqrt{1k^2sin^2\phi }}kF_k(\phi ;k)$$

with the positive derivative

$$f^{}(\phi )=\frac{1}{cos^2\phi \sqrt{1k^2sin^2\phi }}>0.$$

Thus the intersection of the surface, which consists of the gauge orbits, passing through the points of the phase curve (30), with the level of the gauge fixing condition (28) is given by the following equations

$$\begin{array}{c}x_\phi (t_1,t_2)=\overline{x}_\phi \sqrt{\frac{1k^2sin^2\overline{\phi }}{1k^2sin^2\phi (t_1)}}x_A(t_1,t_2)=\frac{tg\phi (t_1)}{\overline{A}}\overline{x}_\phi \sqrt{\frac{1k^2sin^2\overline{\phi }}{1k^2sin^2\phi (t_1)}}\\ x_B(t_1,t_2)=\overline{x}_B+\overline{B}t_2\frac{\overline{x}_\phi }{\overline{B}}\sqrt{1k^2sin^2\overline{\phi }}(F(\phi (t_1);k)F(\overline{\phi };k))+\\ +\frac{\overline{A}\sqrt{J}}{\overline{B}^2}(F_k(\phi (t_1);k)F_k(\overline{\phi };k))\\ A=\overline{A}B=\overline{B}\phi =\phi (t_1),\end{array}$$
(32)
where $`\phi (t_1)`$ is the implicit function given by (31). Using the equations (32) we get the explicit solutions of the differential equations (27) in the following form

$$\begin{array}{c}u^1(t_1,t_2)=\frac{\sqrt{\overline{B}^2J\overline{A}^2sin^2\phi (t_1)}}{\overline{B}}x_B(t_1,t_2)\\ u^2(t_1,t_2)=x_A(t_1,t_2)sin\phi (t_1)+\frac{x_\phi (t_1,t_2)}{\overline{A}}\mathrm{cos}\phi (t_1)+\frac{J\overline{A}}{\overline{B}}x_B(t_1,t_2)sin\phi (t_1)\\ v_1(t_1,t_2)=\sqrt{\overline{B}^2J\overline{A}^2sin^2\phi (t_1)}v_2(t_1,t_2)=\overline{A}sin\phi (t_1),\end{array}$$
(33)
where the initial data in the extended phase space can be obtained with the help of the following formulas

$$\begin{array}{c}\overline{A}=\sqrt{(\overline{v}_2)^2+\frac{\nu ^2}{(\overline{v}_1\overline{u}^2J\overline{v}_2\overline{u}^1)^2}}\overline{B}=\sqrt{(\overline{v}_1)^2+J(\overline{v}_2)^2}\\ \overline{\phi }=arctg\frac{\overline{v}_2(\overline{v}_1\overline{u}^2J\overline{v}_2\overline{u}^1)}{\nu }\\ \overline{x}_A=\overline{u}^2sin\overline{\phi }\frac{J\overline{A}sin^2\overline{\phi }}{\sqrt{\overline{B}^2J\overline{A}^2sin^2\overline{\phi }}}\overline{u}^1\\ \overline{x}_B=\frac{\overline{B}\overline{u}^1}{\sqrt{\overline{B}^2J\overline{A}^2sin^2\overline{\phi }}}\\ \overline{x}_\phi =\overline{A}\overline{u}^2cos\overline{\phi }\frac{J\overline{A}^2sin\overline{\phi }cos\overline{\phi }}{\sqrt{\overline{B}^2J\overline{A}^2sin^2\overline{\phi }}}\overline{u}^1.\end{array}$$
(34)
## 4 Quantization

In this section we consider the canonical quantization of the presented system restricting ourselves to the case of the positive $`J=J_{32}>0`$ , that is to the case when the quantum Hamiltonian $`\widehat{H}_2`$ is an elliptic operator. We realize the Hamiltonians (15) at the quantum level in such a way that the strong involution relation between them is not destroyed by quantum corrections. Then, given these commuting Hamiltonians, we solve the respective eigenvalue problem by requiring that the wave functions have no branch points. The surprising thing is that the spectrum of the Hamiltonian $`\widehat{H}_2`$ turns out to be discrete and the coupling constant $`\nu ^2`$ turns out to be negative. Our starting point is the classical Hamiltonians

$$H_1=\frac{(v_2)^2}{2J}+\frac{\nu ^2}{2J(v_1u^2v_2u^1)^2}H_2=\frac{(v_1)^2}{2}+\frac{(v_2)^2}{2},$$
(35)
which are achieved from (15) by the obvious rescaling

$$v_1v_1u^1u^1$$

$$v_2\frac{v_2}{\sqrt{J}}u^2\sqrt{J}u^2.$$

We perform the quantization of the system (35) in the momentum representation using the polar coordinates on the momentum plane

$$v_1=rsin\varphi v_2=rcos\varphi .$$
(36)
In this setting

$$\widehat{u}_1=i\mathrm{}(sin\varphi _r+\frac{cos\varphi }{r}_\varphi )\widehat{u}_2=i\mathrm{}(cos\varphi _r\frac{sin\varphi }{r}_\varphi )$$

$$\widehat{v}_1=rsin\varphi \widehat{v}_2=rcos\varphi ,$$
(37)
and the quantum Hamiltonians take the following form

$$\widehat{H}_1=\frac{r^2cos^2\varphi }{2J}+\frac{\nu ^2}{2J\widehat{M}^2}\widehat{H}_2=\frac{r^2}{2},$$
(38)
where the operator $`\widehat{M}`$ is defined as

$$\widehat{M}=\widehat{v}_1\widehat{u}_2\widehat{v}_2\widehat{u}_1=i\mathrm{}_\varphi .$$

In order to determine the inverse operator to $`\widehat{M}^2`$ we use the following anzats for the basis wave functions

$$\psi (r,\varphi )=R(r)\mathrm{\Phi }(\varphi )$$
(39)
and represent the factor $`\mathrm{\Phi }(\varphi )`$ as the respective Fourier series

$$\mathrm{\Phi }(\varphi )=\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{\Phi }_me^{im\varphi },$$
(40)
$$\mathrm{\Phi }_m=\frac{1}{2\pi }_0^{2\pi }\mathrm{\Phi }(\varphi )e^{im\varphi }.$$

Then the operator is naturally defined on the functions (39) by linearity if we set that the action of $`\widehat{M}^2`$ on the eigenfunctions $`e^{im\varphi }`$ of $`\widehat{M}`$ looks as follows

$$\widehat{M}^2e^{im\varphi }=\frac{1}{(\mathrm{}m)^2}e^{im\varphi }.$$
(41)
Of course we assume here that the zero mode in the expansion (40) is vanishing $`\mathrm{\Phi }_0=0`$ . It is easy to see that thus defined Hamiltonians (38) are commuting operators and one may pose the respective eigenvalue problem

$$\widehat{H}_1\psi =E_1\psi \widehat{H}_2\psi =E_2\psi .$$
(42)
Using the same anzats (39) for the wave function we see that this problem (42) reduces to a single iteration equation for the coefficients $`\mathrm{\Phi }_m`$

$$m^2(\mathrm{\Phi }_{m+2}+\mathrm{\Phi }_{m2}c\mathrm{\Phi }_m)+b\mathrm{\Phi }_m=0$$
(43)
since the dependence of the wave function (39) on $`r`$ is obviously defined by the second equation in (42)

$$R(r)=\delta (r\sqrt{2E_2}).$$

The iteration equation (43) is equivalent to the following differential equation for the generating function $`\mathrm{\Phi }(z)={\displaystyle \underset{m}{}}\mathrm{\Phi }_mz^m`$

$$(z_z)^2((z^2+\frac{1}{z^2}c)\mathrm{\Phi }(z))+b\mathrm{\Phi }(z)=0,$$
(44)
where

$$c=2(2\frac{JE_1}{E_2}1)b=\frac{2\nu ^2}{E_2\mathrm{}^2}.$$

In what follows it is more convenient to work with the unfolded form of this equation (44), namely

$$z^2(z^4cz^2+1)\mathrm{\Phi }^{\prime \prime }(z)+z(5z^4cz^23)\mathrm{\Phi }^{}(z)+(4z^4+4+b)\mathrm{\Phi }(z)=0.$$
(45)
It goes with out saying that one can hardly find the solutions for the equation (45) in the explicit form. However, we may use the standard results of the complex analysis and conclude that the solutions are analytical functions on $`𝐂`$ besides may be the points where the polynomial $`P(z)=z^2(z^4cz^2+1)`$ is vanishing. Hence, we can qualitatively analyze asymptotical the behaviour of the solutions of (45) in the vicinity of the roots of the polynomial $`P(z)=z^2(z^4cz^2+1)`$ . And thus, we can derive the properties of the desired functions $`\mathrm{\Phi }(\varphi )`$, obtained by restricting the solutions to the unit circle $`\{z=e^{i\varphi }:\varphi 𝐑\}`$ . In fact we have to require that the restriction $`\mathrm{\Phi }(\varphi )`$ is the function on the circle. It means that the respective analytic function $`\mathrm{\Phi }(z)`$ have no branch points in the interior of the unit disc $`\{z𝐂:|z|<1\}`$ or the factors that arise from these branches should combine into the trivial factor $`1`$ if the branch point is not the only one. In addition, we claim that the function $`\mathrm{\Phi }(\varphi )`$ belongs to $`L^1(S^1)`$ . The roots of the polynomial $`P(z)=z^2(z^4cz^2+1)`$ that are worthy of notice are $`z_1=0`$ for arbitrary values of the parameters $`c`$ and $`b`$ $`z_{2,3}=\pm 1`$ for $`c=2,b8`$ and $`z_{4,5}=\pm i`$ for $`c=2`$ and $`b8`$ . It is easy to show that the solutions in the vicinity of the point $`z_1=0`$ behaves like $`\mathrm{\Phi }(z)=z^\mu ,\mu =2\pm \sqrt{b}`$ if $`b4`$ and behaves analytically if $`b=4`$ . Hence $`z_1=0`$ is the branch point for our solutions if $`bn^2,n𝐍`$ . It turns out that there are no other branch points in the interior of the unit disc and thus we conclude that the admissible values for $`b`$ are $`n^2`$ with natural numbers $`n`$ . Considering the points $`z_{2,3}=\pm 1`$ for $`c=2,b8`$ and $`z_{4,5}=\pm i`$ for $`c=2`$ and $`b8`$ it is easy to verify that the solutions of (45) behaves like $`\mathrm{\Phi }(z)=(zw)^\beta `$ where $`w=z_{2,3}`$ or $`w=z_{4,5}`$ respectively and $`\beta ={\displaystyle \frac{3\pm \sqrt{1b}}{2}}`$ . Since $`n𝐍`$

$$\frac{\sqrt{n^2+1}3}{2}>1$$

we conclude that for an arbitrary $`b=n^2,n𝐍`$ there exists a solution $`\mathrm{\Phi }(\varphi )L^1(S^1)`$ . Note that in virtue of the relation

$$b=\frac{2\nu ^2}{E_2\mathrm{}^2}=n^2$$

the coupling constant $`\nu ^2`$ turns out to be negative since $`E_2`$ , being the eigenvalue of the elliptic operator $`\widehat{H}_2`$ is necessarily positive. Thus, we get the following spectrum for the Hamiltonians (38)

$$E_1>0E_2=\frac{2\nu ^2}{\mathrm{}^2n^2},n𝐍\nu ^2<0.$$
(46)
The discrete spectrum of the Hamiltonian $`\widehat{H}_2`$ is not so surprising because if we consider the exponential degeneration of our system (16) and quantize it in the representation of the wave functions $`\psi =\psi (u^2,v_1)`$ we note that the momentum $`v_1`$ defining the coupling constant of the two-particle Calogero-Moser system is quantized and thus the respective quantum Hamiltonian $`\widehat{H}_2`$ also has a discrete spectrum.

## 5 Concluding Remarks

The main lesson we have learned from the above considerations is that we may not restrict ourselves to the Lie group symmetries in constructing the integrable systems via the symplectic reduction method. Of course if we go far beyond the Lie group symmetries the explicit description of the gauge orbits cannot be guaranteed and we cannot be sure that we are able to resolve the equations of motion via the projection method. Nevertheless we suspect that the non Lie group symmetries may be useful for constructing the mechanisms of the symplectic reduction for well known integrable systems. For example, although the mechanisms of the symplectic reduction for the Calogero-Moser systems associated to the root systems of the simple Lie algebras  are known  the projection method for these system has not been developed. We suppose that the projection method for the systems can be presented if we construct the mechanisms of the symplectic reduction in the context of the non Lie group symmetries. It should be also very interesting to perform analogous symplectic reduction starting with the Odesskii bracket, which generalizes the Sklyanin bracket (4) to the case of the arbitrary elliptic $`r`$-matrix  . We suppose that the projection method for the respective systems should be related to solving the equations of motion for some higher dimensional tops  as is the projection method for our system related to the equations of motion for the ordinary $`\mathrm{𝐬𝐨}_\mathrm{𝟑}`$ top. As it has been noted in the introduction the problem of finding the Lax representation for our system seems to be rather intriguing. Here we mention that this problem raises the question of whether it is possible to get the quadratic Sklyanin bracket (4) via the reduction from the linear Poisson bracket. As far as we know the latter question is not answered yet even for the quadratic Poisson brackets, related to the constant triangular $`r`$-matrices . Finally, it is worth mentioning that in our paper we quantize the system starting with the reduced phase space of the model. However, it would be very interesting to obtain the same quantum mechanics via the respective quantum reduction, which may be performed with the help of the BRST approach . We suspect that this quantum reduction should be some how related to the representations of the quantum Sklyanin algebra  as are the quantum Calogero-Moser systems related to the representation of the respective Lie algebras .
Acknowledgments. I am grateful to M.A. Olshanetsky and A.M. Levin for many useful discussions, suggestions and valuable comments. I also acknowledge M.A. Olshanetsky for a careful reading of the manuscript and constructive criticisms concerning the first version of this article. I am indebted to A. Gorsky and A. Zabrodin for useful remarks. I would like to express my sincere thanks to S. Lyakhovich for support and A. Morozov for hospitality and support during my stays in Moscow. It is my pleasure to thank for hospitality Tamm Theory Division of the Lebedev Physical Institute where a part of this work has been done. I also acknowledge discussions with Yu. Chernyakov, I. Gordeli, A. Kotov, S. Loktev, D. Talalaev and A. Zotov. My work is partially supported by RFBR grant 00-02-17-956 and Grant for Support of Scientific Schools 00-15-96557.