Xavier Leoncini

Alberto D Verga
email: verga@marius.univ-mrs.fr

Stefano Ruffo

Hamiltonian Dynamics and the Phase Transition of the XY Model

Keywords: y, 64.60-i, 64.60.Cn

A Hamiltonian dynamics is defined for the XY model by adding a kinetic energy term. Thermodynamical properties (total energy, magnetization, vorticity) derived from microcanonical simulations of this model are found to be in agreement with canonical Monte-Carlo results in the explored temperature region. The behavior of the magnetization and the energy as functions of the temperature are thoroughly investigated, taking into account finite size effects. By representing the spin field as a superposition of random phased waves, we derive a nonlinear dispersion relation whose solutions allow the computation of thermodynamical quantities, which agree quantitatively with those obtained in numerical experiments, up to temperatures close to the transition. At low temperatures the propagation of phonons is the dominant phenomenon, while above the phase transition the system splits into ordered domains separated by interfaces populated by topological defects. In the high temperature phase, spins rotate, and an analogy with an Ising-like system can be established, leading to a theoretical prediction of the critical temperature T KT ≈ 0.855.

I. INTRODUCTION

The two dimensional XY model, also known as the planar spin model, presents many interesting behaviors. Despite the presence of a continuous symmetry group, a particular form of phase transition exists [START_REF] Itzykson | Théorie Statistique des Champs[END_REF][START_REF] Bellac | Des Phénomènes Critiques aux Champs de Jauge[END_REF], which can be characterized by the change in the behavior of the correlation functions. In the low temperature phase these latter have power law decay, showing that the system is in a long range order state; while they decay exponentially at high temperatures, the long range order is broken, even though thermodynamic quantities remain smooth across the transition [START_REF] Berezinskii | [END_REF]. These observations have been interpreted by Kosterlitz and Thouless [4] using an analogy with the transition of a Coulomb gas from a dielectric phase, where charges are bounded into dipoles, to a plasma (conducting) phase where temperature fluctuations destroy the dipoles, and the charges become free. In the XY system the charges are replaced by topological excitations called vortices.

From the analytical point of view, beyond the spin-wave approximation [5], and the Villain model [6], the use of renormalization group in the critical region has been the main issue [7,8]. In order to confirm the analytical results and to satisfy the need for a better understanding of the transition, many numerical studies have been performed. Different dynamics, like Monte-Carlo [9] or Langevin [10], have been introduced. The system has been confirmed to be dominated by spin wave excitations in the low temperatures region. The transition region is then determined using the information on the correlation functions, and is interpreted in terms of dipoles unbinding. However, there is a number of observations in the literature that seems to complicate this simple picture of the phase transition mechanism. Among the different results those worth mentioning are: the visualization of the vortex distribution [9], the presence of domains delimited by topological defects [10], the precise determination of the transition temperature T KT ≈ 0.89 [11], much lower than the transition temperature of the Villain model. Results are also found on the interaction between vortices [10], or the vortex interaction energy [12].

Although the basic mechanism of the Kosterlitz-Thouless transition, in terms of the breaking of vortex dipoles associated to the emergence of a disordered state, is well understood, the observation of the spatial distribution of defects, which is not uniform (defects tend to appear organized into clusters at temperatures slighty larger than the transition temperature), and the presence of large ordered domains where the spins are almost parallel, seem to indicate that the physics of the phase transition is not exhausted by this unbinding process but that some kind of partial local order is present even beyond the transition temperature. The investigation of the system properties near the transition is one of the points addressed in this paper.

On the other hand, the question whether statistical physics is able to describe over a wide range of temperatures the behavior of a classical Hamiltonian system with many degrees of freedom, still remains an open one (see e.g. the literature on the Fermi-Pasta-Ulam model quoted in [13]). Therefore, it seemed interesting to consider the XY model from the point of view of Hamiltonian dynamical systems by adding a kinetic energy term to the XY Hamiltonian. Such kind of approach has proved to be fruitful of interesting information on typical relaxation time scales and collective behaviors in the 1D case [14] and in the mean field approximation [15]. More recently Lyapunov exponents have been computed, confirming the presence of long relaxation times to equilibrium both in the very low and in the very high temperature limits in one dimension [16] and the presence of a change in slope of the maximal Lyapunov exponent vs the energy near the transition temperature in 2D [17] (see also Ref. [18] for a preliminary study of this latter phenomenon).

The present model is one of coupled rotators (spins) sitting on the sites of a square lattice interacting with near neighbors, whose statistical properties are described by the microcanonical ensemble, the total energy being set by the initial conditions. Convergence to a Gibbsian equilibrium distribution is not given for granted in the whole temperature range on accessible time scales. Therefore, a particular attention must be devoted to the temporal behavior of the different quantities which characterize macroscopic thermodynamical properties. One notices, however, that the spatial topology is unchanged and that the existence of a continuous symmetry group of rotation is also maintained in the dynamics. This implies that a Kosterlitz-Thouless type transition must be observed (in the thermodynamic limit) with typically strong finite size effects, like the existence of a non-zero magnetization [5,19].

In this paper we concentrate on a study of those properties of the dynamics of the XY model which reproduce equilibrium features, postponing to a future work the study of non equilibrium effects. We anticipate that we are able to reproduce most of the equilibrium behaviors of macroscopic quantities, which makes microcanonical dynamics competitive with canonical Monte Carlo, since we have not here to extract random numbers, stochasticity being supplied by the intrinsic chaoticity of the model.

The fact that we actually deal with a dynamical system, allows us to develop an original analytical approach to the study of the thermodynamics, which is based on the approximate solution of the equations of motion. This method is based upon the ergodic properties of the dynamics and on the separation of temporal scales present in the spin motion.

In Section II we present the Hamiltonian model, the basic aspects of the numerical computations and we introduce the thermodynamic and dynamical quantities that characterize the state of the system. In Section III we compare the statistical properties of our model to the usual ones, using both analytical and numerical approaches. We derive the thermodynamical properties using the hypothesis that the dynamics is dominated by random phased waves. We then propose, in Section IV, an Ising-like model based on the observation that, above the transition, synchronized regions of spins appear. The relation of the XY with an Ising-like model allows us to describe the high temperature phase and to derive an approximate value of the critical temperature. A brief summary of the main results and conclusions is presented in Section V.

II. THE HAMILTONIAN MODEL, BASIC PROPERTIES AND NUMERICAL COMPUTATIONS

The XY model was introduced in statistical mechanics as a two-dimensional version of the Heisenberg Hamiltonian. The spins are fixed on the sites of a square lattice, and are characterized by a rotation angle

θ i ∈ [-π, π] H XY = J N (i,j) [1 -cos(θ i -θ j )] , (2.1)
where J is the coupling constant (with J > 0 corresponding to the ferromagnetic case, that we study here), i and j label the N sites of a square lattice of side √ N , i.e. i = (i x , i y ) with 1 ≤ i x , i y ≤ √ N of coordinates (x, y). The summation is extended over all i and its neighboring sites j. In the following, without loss of generality, we set J = 1 and the lattice step to unity.

The spins evolve in time, θ i = θ i (t), after adding a kinetic energy term to the XY Hamiltonian,

H = N i=1 p 2 i 2 + H XY , (2.2)
where p i = θi is the spin momentum. The choice J = 1 is equivalent to set time units and to rescale accordingly momentum (in these units the "inertia" is also unity). With this kinetic energy term the spins become in fact rotators, and the XY model becomes a system of coupled rotators. The equations of motion are

θi (t) = - 4 j(i) sin[θ i (t) -θ j (t)] , (2.3)
where the summation is over the four neighbors j of site i. In addition to the energy H = E, there exists a second constant of the motion, the total angular momentum P = i p i , which can be chosen to be zero. We choose periodic boundary conditions in both the x and the y directions. Numerical integration of Eq. (2.3) is performed using the Verlet algorithm, which conserves the energy O(∆t 2 ), ∆t being the time step, but exactly preserves momentum and the symplectic structure. Thermodynamical quantities are computed by averaging over time and over the sites a single orbit (the evolution of the system from a given initial condition). Typically, the system is started with a Gaussian distribution of momenta and with all the spins pointing in the same direction θ i = θ 0 , for reasons that are clarified in the following. No strong dependence on the chosen initial condition in this class was observed, but one could improve statistically our results by averaging over many orbits with different initial conditions and the same energy.

The temperature is computed through the average squared momentum per spin

T = 1 N N i=1 p i (t) 2 , (2.4)
where the overbar stays for temporal averaging.

The thermodynamical state can be characterized by several macroscopic variables: the internal energy per spin h = h(T ) = E(T, N)/N (E being the constant total energy of the system); the magnetization M = M(T, N), which, as mentioned in the Introduction, vanishes for N → ∞, but is sizeable for any finite N; the density of topological defects ρ v , or vortices, which is intimately related to the mechanism underlying the phase transition.

The magnetization M = (M x , M y ) is given by

M ≡ M(t) = 1 N N i=1 cos θ i , sin θ i , (2.5)
and describe the mean orientation of the spin field. A topological defect is identified, as usual, by computing the total angle circulation on a given plaquette (the sum of the four plaquette angles mod[-π, π]); when it equals ±2π this quantity identifies a positive (negative) unitary vortex on the plaquette. The total density of vortices, or vorticity, which is an intensive quantity, is then given by

ρ v = 1 N [i,j] (θ j+1 -θ j ) mod [-π, π] , (2.6)
where [i, j] denotes the sites j of the i-th plaquette. We have studied the temporal evolution and the spatial distribution of the vortices as a function of the temperature. Because of the periodic boundary conditions and since P = 0, the number of positive vortices equals the number of negative ones.

We have performed simulations using various sets of parameters, to study the system behavior depending on the temperature, the number of spins, and the total time to test the stationarity of the relaxed state. Having chosen parallel spins, one has an initially vanishing potential energy, allowing the exploration of the low temperature region through the reduction of kinetic energy. Random angles in the [-π, π] interval would give an energy per spin of 2, which would then remain fixed producing a high temperature configuration, as can be also deduced from the h(T ) function we compute below (see Fig. 1). Therefore, initially, the magnetization is |M| = 1. Other initial conditions, with |M| = 0 were investigated, for instance with half the spins oriented in the x direction, and the other half in the -x direction. After a transient, the system relaxes to the same thermodynamical state. In Fig. 1, we plot the total energy per spin h(T ). The low (T ≪ 1, T ≈ 1 corresponds to the value for which the kinetic energy is of the same order of the potential energy) and high (T ≫ 1) temperature behaviors are easy to understand. At low temperature equipartition of kinetic and potential energies gives h(T ) ≈ T (this is the linear regime, where angles between neighboring spins are small). At high temperature, angles are uniformly distributed in [-π, π], the cosine interaction in (2.2) is negligible with respect to the kinetic energy, then h(T ) ≈ T /2 + 2. These simple arguments suggest that for T ≪ 1, the spin field is almost linear and can be represented as a superposition of waves, which correspond to phonons, while for T ≫ 1, the potential energy is negligeable, and this field becomes a set of almost free fast rotators. The change in the behavior of the h(T ) curve starts around the Kosterlitz-Thouless critical temperature T KT . The peak in the specific heat, related to the second derivative of the energy, occurs instead at a somewhat higher temperature, a fact well documented in Monte-Carlo computations. Figure 2, shows the absolute value of the magnetization as measured using formula (2.5). Although in the thermodynamic limit the magnetization must vanish, in a finite system and for low temperature we expect an observable macroscopic magnetization, which decreases logarithmically with the number of spins as noted by Berezinskii [3]. This slow decrease with system size is hardly observable in Fig. 2, but points corresponding to larger sizes give systematically a smaller magnetization. In the high temperature region, since angles are randomly distributed in space at any time, the magnetization vanishes algebraically with the number of spins.

The behavior of vorticity with temperature, plotted in Fig. 3, is related to the form of the h(T ) curve. Indeed, at low temperature aligned angle configurations are typical and consequently vorticity vanishes, ρ v → 0. In the high temperature regime, spins are randomly distributed and the vortex density approaches the asymptotic value ρ v → 1/3. A dramatic growth of the vorticity occurs in the phase transition region (the temperature in the range T ≈ 0.8 to 1.5). All our results are in perfect agreement with those of Monte-Carlo simulations (see e.g. [20]), some of them are plotted in the Figs. 1 and 2 for comparison. Therefore, microcanonical (Hamiltonian dynamics) and canonical (Monte-Carlo) computations give, at least for the class of initial conditions studied here, the same thermodynamical equilibrium states.

III. DYNAMICAL AND STATISTICAL PROPERTIES

We know that at low temperatures, the main contribution to the canonical partition function comes from the configurations in which the spins are almost aligned, configurations where the angle differences θ i -θ j are small. In such a situation, the equations of motion (2.3) can be linearized, and the spin field is therefore represented as a superposition of linear waves. We will here determine the effective dispersion relation for these waves at low temperature using consistency relations with temperature and internal energy. Let us introduce a representation of the spin angles in the form of a random phased field

θ i = k α k cos(ψ i k ) , ψ i k = k • x i -ω k t + φ k , (3.1)
where the summation is over the wave-

vectors k = (k x , k y ) = 2π(n x , n y )/ √ N , with n x , n y = 1, . . . , √ N
integers; the wave spectrum is given by the α k , and the phases φ k are supposed to be random, uniformly distributed in the circle. We also should have added to (3.1) a term of the form Ω i t reflecting the individual rotation of spins with some frequency Ω i , which in fact can be considered a function of the temperature. However, at low temperature most of these frequencies must vanish, because the energy needed to trigger the free rotation is only reached at energies of the order of h = 2 (or temperatures of the order of T = 1.5, as can be seen in Fig. 1), when a spin, considered as a perturbed pendulum, crosses the separatrix. Then, in the following we put Ω i = 0 in the low temperature (this term becomes important at high temperatures).

Although the main assumption in (3.1), that is the random character of the wave phases, may only be justified a posteriori, we performed some numerical tests which showed that this ansatz is consistent with the properties of the system (at low and also at finite temperatures). We found that the probability distribution of the spin velocities is almost uniform, and that moreover, the motion of one single spin is ergodic (in the sense that its temporal mean and variance coincide with the spatial ones). In addition, if the system were linear this assumption is equivalent to the assumption of thermodynamic equilibrium (with a bath at temperature T ). In fact this is the case at very low and very high T ; of course, in the intermediate range the influence of the nolinearity becomes important, for instance, establishing non-local interactions or some type of self-organization, in which case this assumption would not be necessarly valid.

In formula (3.1) there are two unknown functions of the wave-number which remain to be determined, the wave frequencies ω k and the spectrum α k . Note, in addition, that (3.1) is in fact a change of variables of the angles θ i to the amplitudes α k , where implicitly, a slow temporal dependence may be included, the fast time dependence being assured by the phase ω k t. In this sense, formula (3.1) is rather general and can take into account a large variety of field states.

The definition of the temperature imposes a constraint on the form of the wave spectrum. Indeed, replacing formula (3.1) into the definition of the temperature (2.4) we obtain

T = 1 N N i θ2 i = 1 2 k α 2 k ω 2 k , (3.2)
where • • • stays for averaging over the random phases φ k , and where we used the identity

cos(ψ i k ) cos(ψ i k ′ ) = 1 2 δ k,k ′ ,
being δ k,k ′ the Kronecker symbol. For instance, if equipartition of the energy among different degrees of freedom is assumed, one obtains the usual Jeans spectrum, given by

α 2 k = 2T Nω 2 k , (3.3)
which means that each degree of freedom takes a fraction T /2 of the total kinetic energy.

In fact, as we will show below, the detailed knowledge of the spectrum turns out not to be necessary to derive the dispersion relation and compute thermodynamical quantities. Before proceeding to the computation of the dispersion relation, let us derive a general expression of the energy per spin h,

h = T 2 + 2 - 1 N N <i,j> cos(θ i -θ j ) , (3.4)
where the first term comes from the definition of the temperature, the second term is the constant added for convenience to the Hamiltonian to make the energy vanish at zero temperature, and the last term includes, for each lattice site, a summation over two neighbors, for instance west (in the x direction) and south (in the y direction). When the expansion (3.1) is introduced into Eq. (3.4) the typical contribution to the sum in the third term on the r.h.s. is

U (i,j) N = cos N n=1 A (i,j)
kn sin(B kn (i, j) + φ kn ) ,

where n relabel the modes and

A (i,j) kn = 2 sin k • (x i -x j ) 2
B (i,j) kn = k • (x i + x j ) 2 + ω k t
are local in n and independent on the phases φ kn . Let us split the summation over the wavenumbers into two terms, one containing the phase, φ k N , and the other with the remaining n = 1, . . . , N -1 phases. Since cos A (i,j)

k N sin(B (i,j) k N + φ k N ) = J 0 (A (i,j) k N ) sin A (i,j) k N sin(B (i,j) k N + φ k N ) = 0 , (3.5)
where J 0 is the Bessel function of zero order, we obtain the recursion relation

U (i,j) N = J 0 (A (i,j) k N ) U (i,j) N -1 ,
which readily gives the result

U (i,j) N = n J 0 (A (i,j) kn ) .
We finally get

h = T 2 + 2 - k J 0 (β x ) - k J 0 (β y ) , (3.6)
where β x = 2α k sin k x /2 and β y = 2α k sin k y /2 and the last two products come from north/south and east/west neighbours on the lattice. Eq. (3.6) is the general expression for the energy per spin, as long as the system is dominated by random phased waves. We note that in the case of a symmetric spectrum α kx,ky = α ky,kx (this is the case for the Jeans spectrum, when the dispersion relation satisfies ω kx,ky = ω ky,kx ), the products of zero order Bessel functions are equal,

v 2 = k J 0 (β x ) = k J 0 (β y ) , (3.7)
where v is related to h by h(T ) = T /2 + 2 -v(T ), and is the average potential energy. In order to obtain an explicit form of the energy per spin as a function of the temperature, h = h(T ), we must determine in a self-consistent way the dispersion relation; this may be done by using the equation of motion in its full nonlinear form. But let us begin with the linear case, which will serve as a guide to the self-consistent computation. The linear equation of motion reads θi =

<j> θ i -θ j , (3.8)
Substituting (3.1), one obtains

- k ω 2 k α k cos ψ i k = - k β x [sin(ψ i k + k x 2 ) -sin(ψ i k - k x 2 )] + - k β y [sin(ψ i k + k y 2 ) -sin(ψ i k - k y 2 )] (3.9) = -4 k α k [sin 2 ( k x 2 ) + sin 2 ( k y 2 )] cos ψ i k ,
although the linearity of the equation allows an identification term by term in the summation over k, a different strategy for the solution would be to average (3.9) over all the phases but one, e.g. over the phases φ k ′ = φ k after multiplying both sides by cos ψ i k , in order to isolate only one term in the summation on the left hand side. One finally obtains the linear dispersion relation

ω k = ω 0k = 4(sin 2 k x 2 + sin 2 k y 2 ) , (3.10)
The frequency spectrum is symmetric with respect to the exchange k x ↔ k y ; in the following we assume, without loss of generality, that the spectrum is symmetric with respect to this transformation. This procedure can be generalized to the nonlinear case, although the r.h.s. of Eq. ( 2.3) is no more local in φ k . After substituting the expansion (3.1) into Eq. (2.3), we must average over the phases

φ k ′ = φ k , which we denote < • • • > ′ , terms of the form sin(θ -θ E ) ′ = -sin k β x sin(ψ k + k x 2 ) ′ , (3.11)
where θ E is e.g. the east neighbor of angle θ, having a phase

ψ E = ψ k +k x = k•x+k x -ω k t+φ k .
We now split the summation in (3.11) into a term containing the phase φ k and the others, and develop the sine of the summation of two terms into a product. The averaging of these two terms reduces to cos

k ′ β x sin(ψ ′ k ± k ′ x 2 ) ′ = v 2J 0 (β x ) , (3.12)
because the averaging of the sine term is zero. The two neighbors in the x direction (west and east sites) give the following expression

sin(θ -θ E ) ′ + sin(θ -θ W ) ′ = - v 2J 0 (β x ) sin β x sin(ψ k + k x 2 ) -sin β x sin(ψ k - k x 2 ) . (3.13)
Those in the y direction give exactly the same contribution with β x → β y . On the other hand, the left hand side of Eq. (2.3), the temporal part of the equation of motion, is reduced, after substitution of the random waves expression, to

θi = -ω 2 k α k cos ψ i k . (3.14)
As before, in the linear computation, we multiply both sides of the equation of motion (3.14) and (3.13) (including the similar terms for the north and south neighbors), by cos(ψ i k ) and average over the phase φ k , noting that this average involves first order Bessel functions J 1 . We finally find the desired dispersion relation

ω 2 k α k = 2v sin( k x 2 ) J 1 (β x ) J 0 (β x ) + sin( k y 2 ) J 1 (β y ) J 0 (β y ) . (3.15)
This dispersion relation is nonlinear, i.e. the frequency depends on the spectrum amplitudes, in two ways: implicitly through the α k s in v and explicitly in the Bessel functions. To go further, let us investigate the dependence of the arguments of the Bessel functions on the parameters T and N, and the related limiting form of the total energy per spin h. We know that v, being an intensive thermodynamic quantity, does not depend on the number of spins (the size of the system); moreover, the linear frequency (3.10) is bounded from below, ω 2 0k ∼ sin 2 (k x /2) > O(1/N). We also note, as it is natural for a system near an equilibrium state, that a large number of k modes must contribute to the energy of the system, and then in general we have α 2 k sin k x /2 → 0 when N → ∞, for a large range of temperatures. If this were not the case, the energy would be concentrated in a few high k modes, which is clearly in contrast with the observations. For the Jeans spectrum, one finds specifically

α 2 k sin k x /2 ≈ O(T / √ N )
. Using these approximations we can now develop the logarithm of the product of the Bessel functions in (3.7) to obtain

v = 2 exp - 1 8 k α 2 k ω 2 0k . (3.16)
In the same approximation, for fixed temperatures and large lattices, the dispersion relation (3.15) reduces to

ω 2 k = vω 2 0k 2 . (3.17)
Introducing this expression into the formula (3.16) for v and using the definition of the temperature in terms of random phased waves (3.2) we obtain an implicit equation for the potential energy,

v = 2 exp - T 2v . (3.18)
We notice that the nonlinear dispersion relation reduces to the linear one for T = 0. The nonlinear effects, taken into account in (3.17), appear as a renormalization of the phonon frequency (energy) due to the coupling of the phonons with a thermal bath created by the other phonons like in a mean field. We also note that the factor 2 found in the right hand side of (3.18) and in the exponent can be related to the lattice dimensionality (it results from the addition of the x and y terms). Moreover, as we anticipated, formula (3.18) does not depend explicitly on the form of the spectrum α k . Let us now investigate the implicit equation (3.18) by solving it numerically. In Fig. 4, where the energy density h versus the temperature T is plotted using formulas h = T /2 + 2 -v(T ) and (3.18), we see three branches according to the solutions of (3.18). Although the three branches are in principle acceptable solutions, they would have different "weights". For instance, if one associates to each branch a Gibbsian probability, proportional to exp(-h i /T ), where h i denotes the branch energy i = 0, 1, 2, the lower branch has, at a given temperature, the higher probability. The choice of a Gibbsian probablity distribution of states is not in contradiction with our microcanonical approach. Indeed, the energies h i are related to macroscopic equilibrium states, and the Gibbs probabilities represent the best choice, taking into account the usual constraints of maximal entropy and normalization.

The straight line corresponds to the solution v = 0, which represents a random spin field. Indeed, in the high temperature region, since the potential energy is bounded, we can expect, that the rotators are freely rotating without any order, i.e. that (θ i -θ j ) is random with constant distribution in [0, 2π). This gives us an energy density h = T /2 + 2 and v(T ) = 0. The implicit solution of (3.18) extrapolates thus the good asymptotic behavior of the energy for high temperatures to the low temperature region. This solution exists for any temperature, but in the low temperature region spin configurations associated to this branch should be unstable (in the sense that for general initial conditions, as in our numerical computations, the system cannot evolve to this state).

Below a certain temperature T = T w we also see that two other branches appear. As the lower branch has the higher probability, in the low temperature region, we may expect the lower energy branch, which we call the phonon branch, to be the physical relevant solution, and the third branch to have no physical meaning at all. However, if we now consider that in fact the lattice is populated with two types of species of different energies h 1 (T ) and h 2 (T ) corresponding to the lower (phonons) and upper (which we call vortex) branch respectively, we obtain

h = h 1 e -2h 1 /T + h 2 e -2h 2 /T e -2h 1 /T + e -2h 2 /T , (3.19)
which is almost in total agreement with the numerical data up to T w , as we can see in Fig. 5. The upper branch can be related to "dipole vortices" since it starts at an energy equal to 2 (see Fig. 4) corresponding to a lattice composed by topological defects forming a perfect crystal (θ i -θ j = ±π/2 and j(i) sin(θ i -θ j ) = 0 overall the lattice). This crystal state, which is in fact a periodic array of dipoles, is a stationary solution of the equations of motion at T = 0, around which no linear development can be made (the linear term is zero, vortices can be considered as nonlinear particles). The factor 2 in the exponents of (3.19) takes into account that vortices come in pairs (dipoles at low temperature) on the lattice in order to conserve the total circulation.

The temperature T w determines the upper border of the wave dominated regime, above this temperature the wave branch disappears (as well as the vortex dipole branch). In order to determine precisely T w it is useful to rewrite (3.18), as a formula for the temperature in terms of v, T = -2v log (v/2) .

(3.20)

If we now calculate the extremum of T (v) for 0 < v < 2, we obtain for v = 2/e = T w /2, T = T w = 4/e ≈ 1.47, which is exactly the critical temperature found in [21] using a Hartree-Fock approximation. This temperature also appears with the use of renormalization techniques in [22]. In this case, we note that the equation for v is formally identical to the one for the renormalization constant in this renormalization group calculation. However, the two methods are basically different, and in particular the relation between v and the dispersion relation (3.17), allows us to get more precise informations on the different thermodynamical quantities of the system. For instance, as mentioned above, the relation between v and h, which is different from the relation between the renormalization constant and the energy, leads to a very good description of the h(T ) curve up to T w . Let us now evaluate the magnetization for a finite size system with periodic boundaries using the same approach as before: we assume that the spin field can be represented by a superposition of random phased waves (3.1). We call θ 0 , the average over the lattice of the θ i , which is a constant of motion, since the total momentum is conserved. An average of the magnetization over the random phases brings up to a similar calculation as the one for the energy per spin and leads to,

M = k J 0 (α k ) (cos θ 0 , sin θ 0 ) (3.21)
This expression remains exact as long as the random phases approximation is valid. In order to develop the logarithm of this expression, we have to take into account that α 2 k < O(T ), which implies that the development is only valid in the low temperature regimes T ≪ 1. Moreover a detailed knowledge of the spectrum is also required. Therefore, considering the observed almost flat spectrum of the momentum, we assume the equipartition of the kinetic energy among the modes and use the Jeans spectrum given in (3.3). The absolute value of the magnetization, |M| , is given by the product over the different k of the Bessel functions. In the low temperature regime using (3.17) and (3.3), we obtain for its logarithm,

log ( |M| ) = - k α 2 k 4 = - T 2N k 1 ω 2 k = - T v G (0) , (3.22)
where G is the Green function of the linear wave equation, G(r) = k exp(ik • r)/ω 2 0k , with G(0) = (1/4π) log(2N). The expression for the magnetization is then

M = 1 2N T 4πv
(cos θ 0 , sin θ 0 ) .

(3.23)

A plot of this expression, scaled in order to obtain a function of the temperature, -4πv log(M)/ log(2N), is shown in Fig. 6 (top). Within the errors of the numerical data, the points collapse, at low temperature, to a unique curve. The result (3.23), obtained including the correction due to the nonlinear contribution to the dispersion relation, substantially improves the usual estimation based on the linear wave approximation (see the bottom of the figure). The agreement of the theoretical results and the numerical data can be made more precise by taking into account the different energy branches. Indeed, if we now take into account that the lattice is populated by two types of species, we can consider that only the phonons contribute to the magnetization, vortex dipoles having total magnetization zero. As the density of phonons is given by n 1 = e -2h 1 /T /(e -2h 1 /T + e -2h 2 /T ), the observed magnetization must then be simply n 1 M. The absolute value of this quantity is plotted in Fig. 6 (bottom). The agreement with numerical data is valid up to the critical temperature.

In the same way we compute the averaged absolute value of the magnetization

|M| 2 = 1 N 2 i,j cos(θ i -θ j ) ,
and we obtain the expression

|M| 2 = 1 N 2 i,j k J 0 2α k sin k • x i -x j 2 , (3.24)
which, in the low temperature regime, assuming a Jeans spectrum, and making the same expansion as in (3.22), leads to

|M| 2 = 1 N 2 i,j exp - 2T v G (r ij ) ,
where r ij = x i -x j . A first order calculation leads to the same result as (3.23), implying that the variance Var(M) ≡ |M| 2 -| M | 2 vanishes at leading order in the temperature. However, at finite temperature, the two values can differ, and the variance of the magnetization can have a nonzero value. This variance can moreover be used to characterize the phase transition, since it is the average over the lattice of the correlation function e θ i -θ j ; it is also related to the susceptibility χ ≡ (N/T )Var(M). Indeed, in Fig. 7 we plot the magnetization variance as a function of the temperature, and find that up to temperatures of the order of the Kosterlitz-Thouless T KT , the variance remains very small, and suddenly it grows around T KT . This behavior suggests that the magnetization is distributed almost like a delta function, in the low temperatures regime (all the spins are pointing to the same direction), and that near the transition, the amount of randomly oriented spins increases dramatically. Higher order effects on the variance of the magnetization can be computed using the identity

J 0 (2z sin α) = J 2 0 (z) + ∞ m=1 J 2 m (z) cos(2mα) ,
which gives, after replacing it into (3.24),

Var(M) = | M | 2   1 N 2 i,j k [1 + M k (r ij )] -1   ,
where

M k (r ij ) = m J 2 m (α k ) J 2 0 (α k ) cos (mk • r ij ) .
To the first nonvanishing order (m = 1) we obtain

Var(M) = 2 T vN 2 | M | 2 k 1 ω 4 0k .
This allows us to compute an approximate value of the suceptibility χ, χ ≈ 3.87 10

-3 T v 2 (2N) (1-T /2πv) .
A similar result can be found in [23], in the limit v → 2. Concerning the density of vortices, we can see in Fig. 3 that around T ≈ 1.3 the curve has an inflexion point. This change in the number of defects may be explained as follows. Let us consider a single plaquette, a defect may appear when the four angles are in increasing (decreasing) order, and the last angle is larger (smaller) than π (-π), if the first angle is set to zero. In our model, due to the continuous symmetry group, no particular direction is favored, and therefore the successive differences between the angles θ i -θ j (discrete gradients) are all equivalent. A defect is then obtained when these gradients satisfy (in average) |θ i -θ j | > π/3. At low temperatures the amplitudes of the local gradients are determined by the phonons and become steeper at higher temperatures, reaching a point where they are large enough to generate the vortices. Using the previous results, we calculate the temperature for which this condition is reached,

(θ i -θ j ) 2 = 2 k α 2 k sin 2 k x 2 = T v(T ) = π 2 9 . (3.25)
If we substitute now this result into equation (3.18), we find

T = 2π 2 9 exp - π 2 18 ≈ 1.27 , (3.26)
which is in very good agreement with the position of the inflexion point in the numerical data. One may consider that around this temperature a proliferation of vortices should occur.

IV. THE PHASE TRANSITION

The usual physical picture of the Kosterlitz-Thouless transition is based on the unbinding of vortex pairs: below the critical temperature long range correlations are established by spin waves, the phonons in our Hamiltonian model; above the critical temperature long range order is destroyed by the proliferation of free vortices. This simple physical picture of the mechanism of the transition does not exhaust the complex processes related to the appearence of a kind of self-organization in the system: the vortex distribution is not uniform, vortices form clusters which separate domains of relatively well ordered spins. In this Section we study, using numerical simulations, the spatial distribution of topological defects and the equilibrium configurations of the system, using an Ising-like model, in order to describe the partial order which is still present above T KT .

In the previous Section we showed that the ordered state of the system is dominated by phonons, and that a second type of excitations, the vortex dipoles, is also present. The energies of the two branches meet at a temperature T w , above which both modes disappear, suggesting that a phase transition to a disordered state, might occur. Within this picture, T w would be thus the temperature where the unbinding of vortex dipoles appears, destroying the long range order created by the waves. However, the actual transition temperature T KT is much lower than T w , indicating that other processes, not included in the representation (3.1) of the spin motion, take place.

Indeed, in this low temperature analysis we have omitted the rotation of the spin, excluding terms of the form Ω i t. Moreover, most significantly, we have assumed that the phases of the spin field were random and uncorrelated from site to site, thus neglecting the influence of organized spin motion. Although the proliferation of vortices breaks the long range order, the spin field can still be organized into domains (where the spins rotate synchronously) separated by lines of vortices. The phase transition would be associated in this case, with the appearence of a self-organized state, where, although long range order is absent, a kind of partial local order, established by separated domains of coherent spins, sets in.

We consider this possibility and study the spatial distribution of defects around the critical temperature. By direct visualization of the spin field m i = (cos θ i , sin θ i ), we noted that the number of isolated defects is negligible, and that they often appear to be in small clusters along the domain borders (this may already be seen in Fig. 7 of Ref. [9]). The problem one encounters with direct visualization is that the spins move too fast to allow the observation of coherent structures. To gain some further insight we introduce then a new diagnostic to correctly visualize the spatial structures. We update each spin of the lattice with one fifth of the sum of itself and its four nearest neighbors and repeat this operation a few times, to obtain an effective local magnetization centered on the considered spin. This i is smoother than the initial field m

i ; at these temperatures the amplitude of spin motion is large, and it defines a direction at a given site, weighted by the orientation of the surrounding spins.

The effective spin field is shown in Fig. 8. The upper plots are linear grey scale images of the quantity sin 2 2θ (l) i (see also [10]), where θ (l) i is the angle of the local magnetization, and the bottom ones are images of the quantity |m (l) i |. In the upper plots we can easily locate the vortices by looking at pinching of the darker and brighter areas. They appear to be all bounded in dipoles or chains of dipoles. The presence of these chains can be interpreted as the birth of interfaces separating domains with local order. These latter are more visible in the bottom images, where the brighter regions characterize strong local magnetization, and therefore locally aligned spins; while the darker regions, where the orientation of spins change rapidly with position, indicate the presence of vortex defects and disorder. In the left images, representing a system of N = 128 2 spins at a temperature slightly above the critical temperature, T = 1.02, we see that the disordered regions are highly concentrated and that they tend to connect themselves along lines. This tendency is even more clear in the right images, where we show a system at T = 1.29. We also note that the size of the ordered regions, the domains, are much smaller at higher temperatures. At T = 1.02 the domains occupy connected regions of a size comparable to that of the whole system, at higher temperature the domains are instead confined in localized regions. showing that its fluctuations are fast compared to the evolution of θ i .

At higher temperature the weight of random regions increases, and the size of domains is reduced. In this regime, since the density of potential energy is bounded, most of the energy is kinetic and therefore most of the spins are fully rotating, with similar amplitudes. This is illustrated in Fig. 9, where the typical temporal behavior of one spin in a high temperature field is shown. On the contrary, at low temperatures, angles remain bounded in time. On the other hand the spin momentum takes essentially constant values (positive or negative ones, depending on the sense of rotation). The individual behavior of the spins also reflects into the statistical properties of the system. The angle variance is bounded at low temperatures, in agreement with the random phase approximation, while it steadily increases in time at high temperatures, as we can see in Fig. 10. In analogy with the properties of a perturbed pendulum, the change of the topology of trajectories is related to the crossing of a separatrix. Below a certain value of the energy density h c (h c = h(T KT ) ≈ 1 in our model), the spins are collectively trapped and after the transition, as the behavior of the topological defects suggests, they start to rotate in organized domains, whose size progressively decreases as the energy density increases.

An important consequence of these results is that the system, at these temperatures, can be represented by local ordered regions of synchronized spin motion, and random interfaces where the defects accumulate. To these rotating spin domains one may add some level of fast (with respect to the time scale of rotation) fluctuations characterized by an effective low temperature T , which takes into account the potential energy. In this context we may introduce relevant variables s i , which take the values ±1 according to the sign of spin rotation, θi . It is important to keep in mind that the system is considered to be in its thermodynamical state, and that the temporal average of spin rotation velocity θ2 i = T is, using ergodicity, independent of the spin site. We can then assume that the angle velocity is of the form θi = Ω i ≈ s i √ T , and write

θ i = s i Ω t + θi (t/ǫ) , (4.1)
where Ω = √ T , and ǫ characterize the fast temporal fluctuations of the function θi , i.e. |ǫ| ≪ 1 and | θi (t/ǫ)| ≪ s i Ω t (see Fig. 9). The field θi is in fact similar to the θ i used in (3.1), but associated to another temperature T : in a rotating reference system, attached to each domain (where all the s i 's are equal), it is an effective "low temperature" XY spin field. At the bottom we have the plot for T = 0.4783, while on top it is for T = 2.7905. The numerical data is obtained for a lattice size of N = 256 2 . We clearly see that the variance of θ is an increasing function of time only for T > T KT .

The introduction of the new field s i = ±1, which label each lattice site by the sign of the spin rotation, suggests an analogy between the high temperature XY system and the Ising model. We will exploit this analogy to construct a model aimed at describing the complex properties of the system near the Kosterlitz-Thouless transition, which do not resume to the simple unbinding mechanism of dipole vortices.

In order to map the XY system into an Ising system, we take advantage of the time scale separation between the relative slow rotation and the fast thermal fluctuations. Let us now average the potential energy over a few rotation periods during which the approximation made in (4.1) remains valid. For that purpose, we are reduced to compute the average of terms of the form:

cos(θ i -θ j ) = cos[(s i -s j )Ω t] cos( θi -θj ) -sin[(s i -s j )Ω t] sin( θi -θj ) , (4.2)
where i and j denote two neighbors. The fast temporal behavior of the θi 's assures its fast thermalization at a temperature T and decorrelates the two terms in each product. This allows us to split the time average in two steps, first on the fast time dependence and then on the slower time scale. The first average leads to

cos( θi -θj ) = v( T ) 2 , sin( θi -θj ) = 0 ,
where we use the results of Section III, due to the effective low temperature T within a domain. Then, the mean of cos[(s i -s j )Ωt] is zero if the two spins do not rotate in the same direction (do not belong to the same domain) and one if the two spins corotate, which means that the average of this terms can be represented by a Heaviside function Θ(s i s j ).

We then obtain the expression

v = 1 2N <i,j> cos (θ i -θ j ) = v T 4N <i,j> Θ(s i s j ) , (4.3)
for the mean interaction energy, where the summation is over all neighbors, and the Heaviside function just states that only synchronized, corotating spins contribute to the averaged potential energy. Using the identity Θ(s i s j ) = (1 + s i s j )/2, we finally obtain a new effective Hamiltonian,

H e = N T 2 + 2 - v( T ) 2 - v( T ) 4 <i,j> s i s j . (4.4)
We notice that this effective Hamiltonian is the sum of two different terms, which we denote H T and H I . The first term H T depends on the temperature, the size of the lattice and v( T ), and can therefore be considered as an energy reference term. The second term H I introduces a coupling energy between the spins, and can be recognized as a ferromagnetic Ising-like Hamiltonian whose coupling constant J I (T ) = v( T )/4 is a function of the temperature. Using an Ising terminology, we have now a new population of spins s i = ±1, interacting with their close neighbors on a square lattice. Therefore, under the assumption of ergodicity and taking into account that spins rotate on a well separated time scale with respect to their fluctuations, we have mapped the XY system into a Ising-like model; these two models are linked through a coupling constant dependent on a temperature related to fluctuations. We shall note, that the time dependence of the hamiltonian is masked, but the s i 's are still functions of time, and we still are within the microcanonical ensemble. Invoking again, the thermodynamical equilibrium, and the fact that both the microcanonial and canonical approaches leed to the same thermodynamic limit, we can now continue our study within the canonical ensemble.

As it is well known, the Ising model undergoes a phase transition and generate a spontaneous magnetization M I when J I /T > β Ic ≈ 0.44 [START_REF] Landau | Physique Statistique[END_REF]. In order to know if the present Ising system reaches the transition, we have to investigate the behavior of v( T ). The previously observed domains (Fig. 8), are in agreement with the phenomenology of the Ising model above its transition, the maximum value of the coupling constant being v( T )/4 < 1/2 implies that both models are surely in their high temperature states for T > 1/2β Ic ≈ 1.14.

In the low temperature regime, none of the spins are rotating, T = T and v I = v, all the potential energy is due to the waves. We know that the energy is continuous through the Kosterlitz-Thouless transition and also through the Ising transition. As a consequence, we expect that v(T KT ) = v( T ) at the transition temperature. Moreover, the effective domain temperature T has to be an increasing function of T , and then, v I is a decreasing function of T , in order to obtain h = T /2 + 2 at high temperature, v( t) has to vanish as T → ∞. The equation

J I /T = v( T )/(4T ) = β Ic (4.5)
has therefore a unique solution. A calculation of T KT , using equations (4.5) and (3.18) gives T KT ≈ 0.855, which is in good agreement with the numerical value of T KT ≈ 0.898 found in the literature.

There is however, an important difference between the usual Ising model, and the present one, derived from the XY system. It is the global constraint associated to the conservation of the total momentum. The mean value of the total momentum P = 0, is given by 0 = θi = s i . But this last expression is precisely the Ising magnetization:

M I = 1 N i s i = 0 , (4.6)
which means that only symmetric distribution of positive and negative spins are allowed. This constraint prevents the system from undergoing, as in the Ising model, a second order phase transition, with spontaneous appearence of a macroscopic magnetization below the critical temperature. In the present case, this would mean that a significant fraction of the spins would rotate in a preferred direction, but (4.6) forbids this kind of phase transition, and forces the system to accomodate to a vanishing magnetization. In fact what happens is that, below the transition temperature, the spins cannot fully rotate, domains dissapear, and the long range order is established by a reorganization of spin motions in the form of spin waves, which can still generate a nonintensive finite-size magnetization.

V. CONCLUSION

In this paper we took a dynamical point of view to analyze the statistical properties of the XY model. This approach has the advantage of offering a natural physical framework. From the structure of the evolution equations of the spins (which become coupled oscillators in the linear approximation), one is prompted to consider the phonons (propagating waves in the lattice) as the basic excitation at low temperature. In the thermodynamic limit, and assuming that the system reaches an ergodic equilibrium state, it is reasonable to introduce random phases in the waves. Using this basic mechanism, we obtained a nonlinear dispersion relation containing the fundamental physics of the system. Macroscopic quantities, such as the energy and the (finite size) magnetization, are well described by this method up to temperatures close to the Kosterlitz-Thouless critical temperature.

On the other hand, for high temperatures, our dynamical approach allowed us to map the XY model into an Ising model with a coupling constant depending on the temperature. This relation to an exact model, is beneficial to understand the physics near the transition, and moreover, permits the analytical computation of the critical temperature in excellent agreement with Monte Carlo simulations.

These results lead to a physical description of the Kosterlitz-Thouless transition in terms of the change from an ordered state (long range correlation are established by the wave excitations) to a local ordered state where macroscopic domains are separated by interfaces populated with topological defects. Topological defects have then a tendency to bind themselves in clusters. The generation of disorder by the unbinding of dipoles, is accompanied by the fragmentation of the low energy unique domain of non rotating spins into separated regions of synchronized rotating spins.

FIG. 1 .

 1 FIG. 1. Total energy per spin as a function of the temperature. At low temperature the energy grows as T and at high temperature the energy tends to T /2 + 2. The Kosterlitz-Thouless phase transition occurs at T KT ≈ 0.89. +, x, and * signs refer to the Hamiltonian dynamical simulations for different lattice sizes, N = 64 2 , N = 128 2 , and N = 256 2 respectively; circles refer to canonical Monte-Carlo simulations for a lattice size of N = 100 2 .

FIG. 2 .

 2 FIG. 2. Absolute value of the magnetization as a function of the temperature. Finite size effects give a finite magnetization at low temperatures. At low temperatures the magnetization decreases only logarithmically with the number of spins. +, x, and * signs refer to the Hamiltonian dynamical simulations for different lattice sizes, N = 64 2 , N = 128 2 , and N = 256 2 respectively; circles refer to Monte-Carlo simulations for a lattice size of N = 100 2 .

FIG. 3 .

 3 FIG. 3. Number density of vortices as a function of the temperature normalized to the total number of spins. This quantity is intensive. +, x, and * signs refer to the Hamiltonian dynamical simulations for different lattice sizes, N = 64 2 , N = 128 2 , and N = 256 2 respectively.

FIG. 4 .

 4 FIG.4. Plot of the energy density versus temperature computed from Eq. (3.18). The straight line characterize the random field state, the lower branch, the phonons (spin waves), and the upper branch, lattice excitations of the vortex dipole type.

FIG. 5 .

 5 FIG.5. Plot of the energy density versus the temperature. The circles refer to the numerical data, the dashed line to the phonon branch of Eq.(3.18), and the solid line to the theoretical averaged density of energy, Eq.(3.19).

FIG. 6 .

 6 FIG.6. Finite size effects on the magnetization.The upper plot shows -4πv log(M )/ log(2N ) as a function of the temperature; the line refers to formula (3.23); +, x, and o refer to numerical data for lattice sizes of N = 64 2 , N = 128 2 , and N = 256 2 , respectively. The lower plot shows M (T ), circles are for a N = 256 2 lattice; the thin line is the analytical result and the dotted line is the result using only the contribution of linear spin waves: M = exp(-T log(2N )/8π).

2 FIG. 7 .

 27 FIG.7. Plot of the variance of the magnetization versus the temperature (N = 128 2 ). The transition region (T ≈ 1) is clearly identified.

FIG. 8 .

 8 FIG. 8. Spatial distribution of the local magnetization, showing the domains and the interfaces (see text). (a) and (c) orientation (sin 2 (θ l i /2)), one passes from white to black by a rotation of π/4; (b) and (d) intensity (|m (l) i |), black are disordered regions. (a) and (b): T = 1.02, (c) and (d): T = 1.29.

FIG. 9 .

 9 FIG.9. Typical temporal behavior of a rotator at high temperature. Top, spin θ i (t) at the central site for T = 4.100, measured in units of 2π. The spin rotates in the same direction and with an almost steady speed for long periods of time. Bottom, plot of θi (t) showing that its fluctuations are fast compared to the evolution of θ i .

FIG. 10 .

 10 FIG.10. The spatial variance of θ versus time is plotted for different temperatures. At the bottom we have the plot for T = 0.4783, while on top it is for T = 2.7905. The numerical data is obtained for a lattice size of N = 256 2 . We clearly see that the variance of θ is an increasing function of time only for T > T KT .

ACKNOWLEDGMENTS

We greatfully acknowledge useful discussions with O. Agullo and D. Benisti. We thank P. C. W. Holdsworth for giving us his Monte Carlo data used in Figs. 1 and2. Part of the numerical simulations were performed at the Centre de Calcul Régional, Région Provence-Alpes-Côte d'Azur. A.V. acknowledges INFM-FORUM for financially supporting his visits to Florence. This work is also part of the EC network No. ERBCHRX-CT94-0460.


J Kowalski-Glikman
email: jkowalskiglikman@ift.uni.wroc.pl

A Starodubtsev

Effective particle kinematics from Quantum Gravity

Particles propagating in de Sitter spacetime can be described by the topological BF SO(4, 1) theory coupled to point charges. Gravitational interaction between them can be introduced by adding to the action a symmetry breaking term, which reduces the local gauge symmetry down to SO(3, 1), and which can be treated as a perturbation. In this paper we focus solely on topological interactions which corresponds to zeroth order in this perturbative expansion. We show that in this approximation the system is effectively described by the SO(4, 1) Chern-Simons theory coupled to particles and living on the 3 dimensional boundary of spacetime. Then, using Alekseev-Malkin construction we find the effective theory of particles kinematics. We show that the particles action contains standard kinetic terms and the deformation shows up in the presence of interaction terms. The strength of the interactions is proportional to deformation parameter, identified with Planck mass scale.

Introduction

It is well known that in 3 dimensions gravity is described by a topological field theory, and therefore has a finite number of topological degrees of freedom reflecting the topology of spacetime [START_REF] Witten | 2+1)-Dimensional Gravity as an Exactly Soluble System[END_REF]. If 3d gravity is coupled to point particles, which can be modeled as charged punctures of space manifold, these topological degrees of freedom can be "integrated out" leading to effective, deformed particle kinematics [START_REF] Matschull | Quantum mechanics of a point particle in 2+1 dimensional gravity[END_REF], [START_REF] Meusburger | Phase space structure of Chern-Simons theory with a non-standard puncture[END_REF], [START_REF] Schroers | Lessons from (2+1)-dimensional quantum gravity[END_REF]. This effective system of deformed particles is of Doubly Special Relativity type (see [START_REF] Amelino-Camelia | Testable scenario for relativity with minimumlength[END_REF], [START_REF] Amelino-Camelia | Relativity in space-times with short-distance structure governed by an observer-independent (Planckian) length scale[END_REF], [START_REF] Bruno | Deformed boost transformations that saturate at the Planck scale[END_REF] for the original proposal and [START_REF] Kowalski-Glikman | Introduction to doubly special relativity[END_REF], [START_REF] Kowalski-Glikman | Doubly special relativity: Facts and prospects[END_REF] for reviews), being characterized by two scales, velocity of light c and Planck mass κ. Similarly it turns out that integrating out topological gravitational degrees of freedom in the case of gravity coupled to fields leads to effective field theory on non-commutative spacetime [START_REF] Freidel | Effective 3d quantum gravity and noncommutative quantum field theory[END_REF], [START_REF] Freidel | Ponzano-Regge model revisited. III: Feynman diagrams and effective field theory[END_REF]. As a result spacetime symmetries of deformed systems become quantum symmetries, being described by quantum groups, instead of the standard Lie groups.

The question arises if something similar could happen in the case of physical, 4-dimensional gravity coupled to particles and/or fields? Certainly, in 4d gravity possesses local degrees of freedom exhibited by Newtonian interactions and gravitational waves, for example. However, as stressed in the context of recent investigations in quantum gravity phenomenology (see [START_REF] Amelino-Camelia | Quantum Gravity Phenomenology[END_REF] for recent review and references to earlier works) our best hope to see signals of quantum gravity is to look for high-energetic events events (like scattering of ultra high energy cosmic rays) in description of whose local degrees of freedom of gravity play no role. If we expect to see some quantum gravity imprints there, to describe them we must look for "no local gravity limit of gravity", hoping that in this limit the effective theory behaves like gravity in 3d, effectively deforming particles and fields kinematics. This is why we make use of DSR-like test theories, that predict deformation of spacetime-symmetries characterized by Planck-scale deformation parameter to describe quantum gravity signal that may be detected in foreseeable experiments.

It is therefore of interest to look for a formulation of gravity such that one keeps good control over the limit, in which local gravitational degrees of freedom are not present. In this limit gravity is described locally by its maximally symmetric vacuum state: Minkowski space for zero cosmological constant, and (anti) de Sitter space in the case of the (negative) positive one. Coupling such a theory to point particles and then taking the limit makes it possible investigate the effective behavior of the particles. One then can ask the question if the quantum gravity scale is still present in this effective theory, acting as a deformation? Naively the answer would be in the negative, and after taking a limit we would end up with the standard theory of particles moving on Minkowski space of special relativity (or on (anti) de Sitter space when the cosmological constant is non-zero.) On the other hand, the experience of 3d gravity suggests that even in the limit the effective theory might be deformed by the presence of Planck scale. This is the problem that we would like to investigate in this paper.

As said above, in 4 dimensions gravity is certainly not a topological theory. However, it can be nevertheless described by a BF topological theory with de Sitter SO(4, 1) gauge group, appended by a small term that breaks gauge symmetry down to Lorentz SO [START_REF] Schroers | Lessons from (2+1)-dimensional quantum gravity[END_REF][START_REF] Witten | 2+1)-Dimensional Gravity as an Exactly Soluble System[END_REF]. The presence of this symmetry breaking term switches on the local degrees of freedom of gravity, while in the limit, in which this term vanishes the theory becomes topological, as in 3d. Such formulation of gravitational field has its roots in the proposal of MacDowell and Mansouri [START_REF] Macdowell | Unified Geometric Theory Of Gravity And Supergravity[END_REF] and has been recently investigated in depth in [START_REF] Smolin | General relativity with a topological phase: An action principle[END_REF] and [START_REF] Freidel | Quantum gravity in terms of topological observables[END_REF].

The construction presented in these works is as follows. The building blocks are the SO(4, 1) connection A IJ with curvature F IJ = dA IJ + A I K ∧ A KJ and the so(4, 1) valued two-form field B IJ (I, J = 0, . . . , 4), in terms of whose the action takes the form

S = B IJ ∧ F IJ - β 2 B IJ ∧ B IJ - 1 2 B IJ ∧ B KL IJKLM V M (1)
In this action a constant algebra element V M enforces breaking the gauge symmetry group down to SO(3, 1) (the subgroup of SO(4, 1) that leaves V M invariant.) Without loss of generality we can take this vector to be V M = (0, 0, 0, 0, V 4 ) = (0, 0, 0, 0, α 2 ) so that the action takes the form

S = B IJ ∧ F IJ - β 2 B IJ ∧ B IJ - α 4 B IJ ∧ B KL IJKL4 (2)
The first two terms in the action describe the topological field theory, and in the limit α → 0 the only solution of equations of motion is de Sitter space, the topological vacuum of the full theory. Notice that both the topological lagrangian and the gauge breaking one are manifestly diffeomorphism invariant, and thus the perturbation theory in α, around topological vacuum, corresponding to topological theory at α = 0, is going to be manifestly diffeomorphism invariant as well.

Remarkably, it can be shown [START_REF] Freidel | Quantum gravity in terms of topological observables[END_REF] that if one decomposes the connection A a4 =1 e a , with being a constant of dimension of length, A ab = ω ab , a, b = 0 . . . 3, where e a and ω ab are tetrad and Lorentz connection one forms, respectively, after solving for B this action reproduces, up to topological terms, to the standard Einstein-Cartan one

S P = - 1 2G (R ab ∧ e c ∧ e d - Λ 6 e a ∧ e b ∧ e c ∧ e d ) abcd - 2 γ R ab ∧ e a ∧ e b (3)
The last term in this action does not modify field equations if torsion vanishes. The physical constants: Newton's constant G, cosmological constant Λ and Immirzi parameter γ are related to the dimensional parameters of the original action (2) and the scale as follows

γ = β α , 1 2 = Λ 3 , G = 3α(1 -γ 2 ) Λ = 3β(1 -γ 2 ) Λγ (4)
One can couple gravity described by the action (2) to particles in a rather straightforward way [START_REF] Freidel | Particles as Wilson lines of gravitational field[END_REF] (see also [START_REF] Fairbairn | On gravitational defects, particles and strings[END_REF] for recent discussion.) Each particle with (dimensionless) mass µ = m and spin s at rest is described by an appropriate element of the so(4, 1) algebra 1 D = µT 04 + sT 23 (5

)
where T 04 and T 23 are "translational" "rotational" generators of so(4, 1) algebra, respectively. Then the lagrangian describing the particle at rest is simply

S rest = dτ Tr(DA τ ), A τ ≡ A µ żµ (τ ) (6)
where z µ (τ ) is the particle worldline, and τ is the affine parameter. The action [START_REF] Amelino-Camelia | Relativity in space-times with short-distance structure governed by an observer-independent (Planckian) length scale[END_REF] breaks gauge invariance at the particle worldline. However SO(4, 1) gauge transformations are just translations and Lorentz transformations, so acting by gauge transformations on the particle at rest just makes the particle moving. In this way the gauge degrees of freedom of gravity on worldline become dynamical degrees of freedom of the particle and the action describing arbitrarily moving particle has the form

S particle = dτ Tr(DA h τ ), A h τ ≡ h -1 A τ h + h -1 ∂ τ h, h ∈ SO(4, 1) (7)
It can be shown [START_REF] Freidel | Particles as Wilson lines of gravitational field[END_REF] that this action leads to correct equations of motion for the particle and generalized Einstein equations with mass and spin of the particle being the source of curvature and torsion, respectively. Notice that the action (7) easily generalizes to the case of a finite number of particles in which case it reads

S particles = i dτ Tr(D i A h τ ), (8)
For further discussion of this formulation of gravity coupled to particles see [START_REF] Freidel | Particles as Wilson lines of gravitational field[END_REF] and [START_REF] Kowalski-Glikman | Can we see gravitational collapse in (quantum) gravity perturbation theory?[END_REF].

The action of gravity (2) coupled to particles ( 8) is a convenient starting point to address the question raised above. The point is that the presence of the parameter α in the action (2) makes the "no gravity limit" easy to control, and to set up a perturbative theory. In the next section we consider the quantum perturbative expansion in parameter α and we argue that the zeroth order of this expansion is described by a holographic theory living on the boundary of spacetime. Then, in the following two sections, we will show how such boundary theory reduces to effective particles dynamics. In the final section we discuss the obtained results.

2 From gravity in 4d to 3d Chern-Simons theory

The purpose of this section is to demonstrate that the zeroth order approximation of perturbation theory of quantum gravity with particles around BF topological quantum field theory is described by Chern-Simons theory coupled to point sources. We will not study higher order corrections here.

For shortness let us write the action (2) as

S = L BF + α 4 L I , (9)
where now L BF is a lagrangian of topological field theory possibly coupled to particles, and L I is a symmetry breaking interaction term. The general perturbative expression for the partition function coupled to arbitrary finite number of particles (where we neglect all the interactions except gravitational) looks like

Z({g pi }, {g p f }) = DADB n (iα) n n! L I (x) n exp i 1 β M L BF = n (iα) n n! Z n ({g pi }, {g p f }), (10)
Here g pi and g p f are SO(4, 1) group elements labeling initial and final positions and orientations of p-th particle with respect to the selected reference point. They can be obtained as holonomies of connection A between corresponding points (see [START_REF] Kowalski-Glikman | Can we see gravitational collapse in (quantum) gravity perturbation theory?[END_REF] for detailed discussion.)

The expression [START_REF] Kowalski-Glikman | Doubly special relativity: Facts and prospects[END_REF] is formal, of course, and we have to define it precisely. First we must specify the measure DA in the path integral. Since the action describes a system with gauge symmetry the conventional way to proceed would be to introduce a gauge-fixing term. An immediate problem with this approach has been pointed out in [START_REF] Cattaneo | Four-dimensional Yang-Mills theory as a deformation of topological BF theory[END_REF] for Yang-Mills theory. Namely, since the interaction term breaks the symmetry, higher order terms in perturbative expansion have less gauge symmetry than the free action. On the other hand, if we do not fix all the gauge symmetries of the free action we cannot construct a propagator because of non-invertibility of the quadratic form. The approach of [START_REF] Cattaneo | Four-dimensional Yang-Mills theory as a deformation of topological BF theory[END_REF] has been to introduce an auxiliary field that turns gauge degrees of freedom into physical ones already at zero order. Such procedure however trivially reduces the perturbative expansion to the standard one around fixed background. This is appropriate in the case of Yang-Mills theory, but cannot be applied in the case of a theory that is supposed to be background independent.

In this paper we consider a different approach which does not use a gauge fixing. To define a path integral without gauge fixing one has to explicitly construct the reduced phase space spanned by the complete set of gauge invariant observables and define a measure on it. In most situations the later is not possible for technical reasons, e.g. it is certainly impossible to construct all the diffeomorphism invariant observables of four dimensional General Relativity. Fortunately, this turns out to be possible for any finite order terms of the perturbation theory considered here. This happens because the starting point for the expansion is a topological field theory whose reduced phase space (moduli space) is finite dimensional. Then, in finite order of perturbative expansion the dimensionality of the moduli space is getting larger, but remains finite. Therefore, in any finite order of expansion the moduli space can be, explicitly constructed and the measure on it can be defined.

To show this let us replace the integral over the symmetry breaking lagrangian in [START_REF] Kowalski-Glikman | Doubly special relativity: Facts and prospects[END_REF] by its Riemann definition. Divide the manifold M into N cells M i where each cell is sufficiently small so that every field φ on which the lagrangian depends can be considered constant within each cell

φ(x) Mi = φ i
The integral is a sum of contributions from every cell

L Int (φ(x)) = N i S I (φ i ), where S I (φ i ) = Mi L I (φ i ).
Here we took into account that the lagrangian L I is a density and the volume element is already contained in it.

Consider the contribution of a cell M i to the partition function at the first order Z 1i = Dφe iS T S I (φ i ) and compare it with the contribution Z 1i from another cell M j . One can always find a diffeomorphism x → x such that M i → M j . The fields also transform under this diffeomorphism φ → φ . Due to diffeomorpfism invariance of the interaction term

S I (φ i ) = S I (φ j ).
The free action is also diffeomorphism invariant and we assume that we can define a diffeomorphism invariant measure of the path integral

S T (φ) = S T (φ ), Dφ = Dφ .
As a result, the contribution to the path integral from different cells is equal, Z 1i = Z 1j , and the sum over cells becomes trivial.

At higher order one has to distinguish the situations when the cells on which the interaction term is applied are the same or distinct.

Z n = N n Dφ i S I (φ i )e iS T + ... + N DφS I (φ i ) n e iS T
In the limit of large number of cells the first term, where interaction is applied to the distinct cells, will be dominating.

The above argument is analogous to that of coordinates-independence of npoint functions of diffeomorphism invariant theories (see e.g. [START_REF] Rovelli | Graviton propagator from background-independent quantum gravity[END_REF] and references therein).

So far we were using only diffeomorphism invariance of the free action. Its topological invariance allows us to give a complete definition of every term in the expansion of the the path integral. The free action can be replaced by its discretized version S T → i S T (φ i ), and the measure can be specified as Dφ → i dφ i . The resulting path integral is finite dimensional and due to topological invariance of the free action is independent of the discretization.

Let us now specify the model in [START_REF] Kowalski-Glikman | Introduction to doubly special relativity[END_REF] to be (2) with particles coupling in the α → 0 limit, i.e. being defined, after integrating out B-field, by

L BF = 1 β F IJ ∧ F IJ + i D IJ i A 0IJ δ 3 (x -x i (τ )) (11)
In the formula above D i = µ i T 04 + s i T 23 is an algebra element defining mass and spin of i-th particle (both µ i = m , where m is the physical mass and is the length scale of the model (4), and s i are dimensionless here), x i (τ ) is a timelike particle trajectory which can be taken arbitrary due to topological invariance of the model. In the present paper we are studying only the zeroth order contribution, and therefore we do not consider the interaction term L Int here.

Due to Bianchi identity the integral of the first term in [START_REF] Freidel | Effective 3d quantum gravity and noncommutative quantum field theory[END_REF] after integration by parts reduces to a Chern-Simons action on the boundary of the original manifold M.

1 β M F AB ∧ F AB = 1 β ∂M Y CS (A) (12)
Below we consider the boundary ∂M to be the direct product of a (punctured) sphere S 2 n with the real line R. The second term in [START_REF] Freidel | Effective 3d quantum gravity and noncommutative quantum field theory[END_REF] breaks the gauge symmetry at the location of the particles thus promoting some of the gauge degrees of freedom to the physical ones. We can include the latter explicitly by substituting into the action the connection in an explicit gauge transformed form:

A → h -1 dh + h -1 Ah (13)
At the location of the particles h(x i ) ≡ h i become physical degrees of freedom.

To make a link with canonical formulation we will rewrite the resulting action explicitly, decomposing the connection into .spacelike and timelike components.

S BF [A S ,A 0 , h i ] = R dx 0 S 2 n k 4π ∂ 0 A S ∧ A S - R dx 0 n p=1 D i , h -1 i ∂ 0 h i (14)
+ M d 4 x A 0 , k 2π F S δ(∂M) - n p=1 T i δ (3) (x -x (i) )dx 1 ∧ dx 2 ∧ dx 3
Let us explain the notation used in the formula above. The Chern-Simon connection one-form A is decomposed into time and space part, to wit

A = A 0 dx 0 + A S , F = dx 0 ∧ (∂ 0 A S -d S A 0 + [A 0 , A S ]) + F S ( 15
)
and thus F S is the space part of the curvature two-form. We take the Chern-Simon coupling constant to be k/4π = 1/β. D i is the charge carried by the particle, which in our case will be just its mass D i = µ i T 04 , while T i is the algebra element defined by

T i = h i D i h -1 i ( 16
)
Finally < •, • > denotes the invariant, symmetric, bilinear form on the gauge group algebra, which below will be taken to be a trace of product of appropriate matrices, normalized such that < 1 >= 1.

The one-dimensional delta function δ(∂M) reducing the expression to the boundary is defined by

M δ(∂M)( ) = ∂M ( ) (17)
In the action (15) the first line is the kinetic term, while the second is the constraint

k 2π F S δ(∂M) - n i=1 T i δ (3) (x -x (i) )dx 1 ∧ dx 2 ∧ dx 3 = 0 ( 18
)
Notice that because of topological invariance of the model the deformation of the trajectory of a particle x (i) → x i (x (i) ) does not change the value of the physical degrees of freedom of the model. The later are encoded in the group elements h i at the location of the particle and not in the position of the trajectory. In particular, by such deformation one can map the whole particle trajectory x (i) on the boundary ∂M of the manifold without changing h i .

However, the constraint (18) distinguishes the boundary, and in fact, mapping the particle trajectories to the boundary is the only way to satisfy it for a nonzero value of T i . Indeed, if some part of the particle trajectory does not belong to the boundary, at such points the first term in ( 18) is zero, while the second term is equal to a constant times T i . In such situation the constraint [START_REF] Fairbairn | On gravitational defects, particles and strings[END_REF] would force T i to be zero. As due to Bianchi identity the charge T i has to be conserved, T i = 0 along some part of trajectory means T i = 0 along the whole trajectory which in turn means that there is no particle. Thus, the only way to introduce a particle satisfying constraint [START_REF] Fairbairn | On gravitational defects, particles and strings[END_REF] is to map the whole particle trajectory on the boundary.

Mapping particle trajectories on the boundary is analogous to introducing the Dirac string singularity for magnetic monopole in electromagnetism. Through the Dirac string the magnetic flux can reach from the boundary to the point at which the monopole is located. In our model the situation is even simpler. In topological field theory the position of a monopole with respect to a manifold coordinates has no physical relevance. Therefore connecting a monopole to the boundary with a string is the same as placing the monopole on the boundary.

We have seen that the equations of motion of the model, namely the constraint [START_REF] Fairbairn | On gravitational defects, particles and strings[END_REF], force the particle trajectories to lie on the boundary ∂M of the manifold. So far our considerations were classical. But we should expect these results to hold also in quantum theory because the equations of motion which are constraints on initial data, such as [START_REF] Fairbairn | On gravitational defects, particles and strings[END_REF] hold in quantum theory exactly. To see this let us substitute the action (14) into the path integral over the connection A and perform integration over A 0 . Among the spatial components of connection A we will distinguish A S -the components lying within the boundary and A 3 -the component transverse to the boundary. As the action does not depend on A 3 we will not include integration over it. We obtain

W 0 ({g pi }, {g p f }) = DAe iS BF [A S ,A0,hi] = DA S e i R R dx 0 R S 2 n k 4π ∂0A S ∧A S - R R dx 0 P n p=1 Di , h -1 i ∂0hi x∈∂M δ k 2π F S,12 - n i=1 T i δ (2) (x -x (i) ) x ∈∂M δ n i=1 T i δ ('3) (x -x (i) ) (19)
The last factor in the path integral in [START_REF] Kowalski-Glikman | Can we see gravitational collapse in (quantum) gravity perturbation theory?[END_REF],

x ∈∂M δ n p=1 T i δ ('3) (x -x (i) )
forces the partition function to be zero whenever we have a particle away from the boundary. If all the particles are on the boundary this factor is (an infinite) constant which can be absorbed in the normalization of the partition function. Therefore the path integral [START_REF] Kowalski-Glikman | Can we see gravitational collapse in (quantum) gravity perturbation theory?[END_REF] has only contributions from the particles sitting on the boundary where the last factor can be ignored. The resulting path integral is a path integral for the action

S[A S , A 0 , h i ] = R dx 0 S 2 n k 4π ∂ 0 A S ∧ A S - R dx 0 n i=1 D i , h -1 i ∂ 0 h i (20) + R dx 0 S 2 n A 0 , k 2π F S - n i=1 T i δ (2) (x -x (i) )dx 1 ∧ dx 2
of Chern-Simons theory on the boundary coupled to point charges. We will study this action in the next section. The relation between quantum gravity in the bulk and Chern-Simons theory on the boundary was first studied in [START_REF] Smolin | Linking topological quantum field theory and nonperturbative quantum gravity[END_REF] and extended to include translational degrees of freedom in [START_REF] Starodubtsev | Topological excitations around the vacuum of quantum gravity. 1. The Symmetries of the vacuum[END_REF]. For the amplitudes considered in this paper the correspondence is precise. This picture is also supported by spinfoam studies [START_REF] Barrett | Observables in the Turaev-Viro and Crane-Yetter models[END_REF] where it was shown that the invariants of four dimensional Crane-Yetter model are equivalent to Turaev-Viro invariants on the three dimensional boundary.

From Chern-Simons to particle kinematics

Let us summarize what was achieved above. We showed that in the topological limit the action for gravity coupled to the particles is equivalent to the Chern-Simon action for the gauge group SO(4, 1), coupled to particles, carrying the charges (masses and spins) of the same group. The 3 dimensional manifold on which this theory is defined is assumed to be a product of a punctured 2sphere S 2 n , with each puncture corresponding to the particle, with real line R representing time. Thus in the zeroth order of perturbative expansion, i.e., in the "no local gravity limit of quantum gravity" our theory is described by a holographic quantum SO(4, 1) Chern-Simons theory with particles.

In what follows we recall the construction of Alekseev and Malkin [START_REF] Alekseev | Symplectic structure of the moduli space of flat connection on a Riemann surface[END_REF] (see also [START_REF] Meusburger | Phase space structure of Chern-Simons theory with a non-standard puncture[END_REF], from which we borrowed the notation). In the case when the topology of the boundary of spacetime is simple (no handles) and thus reflects only the presence of particles exhibited by punctures, the topological degrees of freedom of gravity can be absorbed by particle's. Thus effectively we obtain a (possibly deformed) theory of particles kinematics. As stressed above this theory could only depend on gauge group elements, as positions in spacetime do not play any role.

In what follows we will consider only the classical theory. The symplectic form on the space of gauge field and particles configurations can be easily found from the action [START_REF] Cattaneo | Four-dimensional Yang-Mills theory as a deformation of topological BF theory[END_REF] and reads

Ω = k 4π S 2 n δA S ∧ δA S + n i=1 δ D i , h -1 i δh i (21)
This symplectic form is subject to the constraint

k 4π F S = n i=1 T i δ (2) (x -x (i) )dx 1 ∧ dx 2 (22)
which makes the curvature zero everywhere except for the positions of the particles. It follows that the connection A S takes simple form at appropriately defined submanifolds of S 2 n . Alekseev and Malkin [START_REF] Alekseev | Symplectic structure of the moduli space of flat connection on a Riemann surface[END_REF] define these submanifolds as follows. Take a point p 0 away from the punctures and construct loops m i with this base point around each puncture (Figure 1). Along these loops we will later calculate holonomies. Now cut the surface along the loops and remove the so obtained discs. As a result we obtain n punctured discs Q i and the polygon Q 0 with no punctures inside, whose boundary will contain exactly n edges which we denote m i and vertices 2). Now the first term in ( 21) can be decomposed as follows

p i p 0 m1 → p 1 • • • p n-1 mn → p 0 (Figure
k 4π S 2 n δA S ∧ δA S = k 4π Q0 δA S ∧ δA S + k 4π n i=1 Qi δA S ∧ δA S ( 23
)
The virtue of the decomposition ( 23) is that on each region the form of the connection is quite simple. Consider the region Q 0 first. Since there are no punctures in this simply connected region the constraint [START_REF] Smolin | Linking topological quantum field theory and nonperturbative quantum gravity[END_REF] tells that connection A S is trivial there Then by direct calculation one can convince oneself that the first integral in [START_REF] Starodubtsev | Topological excitations around the vacuum of quantum gravity. 1. The Symmetries of the vacuum[END_REF] reduces to the boundary one

A S | Q0 = γ 0 d S γ -1 0 (24)
k 4π Q0 δA S ∧ δA S = k 4π ∂Q0 δγ -1 0 γ 0 , d(δγ -1 0 γ 0 ) (25)
The contribution from the discs Q i can be found by similar analysis. The only difference is that now we have to do with a region with puncture inside, carrying the charge, being an element of gauge algebra. Since the curvature must have the delta singularity at the puncture, the most general connections on the disc must be given by gauge transformations of the canonical ones. The latter are defined by

B i = 1 k Di dφ i , Di = k 2π µ i T 04 = k 2π D i
where φ i are angular coordinates, that along with the radial ones (defined such that the singularity corresponds to r i = 0) can be introduced on the disks globally. Thus

A S | Qi = 1 k γ i Di dφ i γ -1 i + γ i dγ -1 i (26)
It is easy to see that since

ddφ i = 2πδ(x -x i )dx ∧ dy this connection solves the constraint (22) if γ i (x i ) = h i . Plugging (26) to one of the Q i integrals in (23) we find Qi δA S ∧ δA S Figure 2: = -δ Di , h -1 i δh i + k 4π ∂Qi δγ -1 i γ i d δγ -1 i γ i - 1 2π ∂Qi δ Di δγ -1 i γ i dφ i (27)
Notice that the sum of the first terms in (27) cancels exactly the second term in [START_REF] Rovelli | Graviton propagator from background-independent quantum gravity[END_REF]. Thus the symplectic form is a sum of (25) and the last two terms in [START_REF] Meusburger | Generalised Chern-Simons actions for 3d gravity and kappa-Poincare symmetry[END_REF]. Its form can be simplified further by observing that the connections on the segments of ∂P 0 = i m i must be equal to the connection on appropriate ∂P i , in order to make connection continuous. We have therefore

A S | mi = γ 0 d S γ -1 0 mi = 1 k γ i Di dφ i γ -1 i + γ i dγ -1 i ∂Qi = A S | ∂Qi (28)
This equation can be solved to give

γ -1 0 mi = N i exp 1 k Di φ i γ -1 i mi , dN i = 0 (29)
Substituting this expression to [START_REF] Alekseev | Symplectic structure of the moduli space of flat connection on a Riemann surface[END_REF] and then, along with ( 27) to (23) will give us the final result. Before doing this final step, let us consider again the polygon P 0 . It has n vertices p 0 , p 1 , . . . , p n-1 connected by edges m i . Since the connection is trivial on the polygon, the parallel transport along an edge is given simply by a product of the gauge parameters at the beginning and the end of the edge

P T pi-1→pi = γ 0 (p i )γ -1 0 (p i-1 )
But because of the continuity of the connection the values of γ 0 at the vertices of the of the polygon P 0 are just the products of holonomies M i of the connection of the disks

γ 0 (p i ) ≡ K -1 i = M i M i-1 • • • M 1 , γ 0 (p 0 ) ≡ K -1 0 = 1 (30)
Next, since the connection γ i is single valued on each disk it is clear from [START_REF] Meusburger | Poisson structure and symmetry in the Chern-Simons formulation of (2+1)-dimensional gravity[END_REF] that each holonomy is of the form

M i = g i C -1 i g -1 i , C i = exp 2π k Di , g i = γ i (p i ) (31)
Now we are ready to present the central result of this section, the Alekseev-Malkin theorem [START_REF] Alekseev | Symplectic structure of the moduli space of flat connection on a Riemann surface[END_REF]: The symplectic form Ω (21) gets contributions only from vertices of the polygon Q 0 and reads

Ω = k 4π n i=1 C i g -1 i δg i C -1 i ∧ g -1 i δg i -δK i K -1 i ∧ δK i-1 K -1 i-1 (32)
The first term in this expression is a deformed kinetic term (a deformation of the standard free particle symplectic form δp µ ∧δq µ , and the second describes some "topological interaction" of the particles. Notice that as expected the final form of the symplectic structure does contain no trace of the connection; there is as many degrees of freedom as the number of particles, each described by a gauge group element. In this way the topological degrees of freedom of gravity has been "eaten" by the particle's ones.

We will discuss the theory of particles provided by this symplectic structure in the next section.

The deformed particle

Having obtained the symplectic form for single particle (32), let us try to understand in which sense it describes a deformed particle. To do that let us first compare it with with an analogous expression for free particle moving in 4 dimensional de Sitter spacetime in the limit of vanishing cosmological constant.

The Lagrangian for such particle is given by equation ( 7) and reads (following [START_REF] Freidel | Particles as Wilson lines of gravitational field[END_REF] we represent SO(4, 1) generators by gamma matrices, T IJ = γ IJ ≡ 1/2[γ I , γ J ], I, J = 0, . . . , 4 with γ 5 matrix denoted here γ 4 )

L =< DA h τ >= µγ 04 (h -1 A (0) τ h + h -1 ḣ) , (33)
where µ = m with being the cosmological length scale, related to the cosmological constant, cf. (4). Since de Sitter background is gauge equivalent to zero configuration we take A (0) τ = 0 Then (33) takes the form

L = µγ 04 h -1 ḣ (34)
Let us decompose the SO(4, 1) group element into translational and Lorentz SO(3, 1) parts (using Cartan decomposition)

h = T L, T = 1 + q a γ a4 + O 1 2 , (35)
Substituting this into (34) and keeping only the leading term after simple calculation we find (the L -1 L term is proportional to γ ab and cancels under the trace)

L = µ γ 04 L -1 γ a4 L qa , a = 0, . . . , 3 (36)
After some γ matrices algebra this expression gives

L = p a qa , p a = m γ 0 L -1 γ a L (37)
and it follows that p a defined above is to be identified with particle momentum. Notice that it follows from (37) that

p a γ a = mLγ 0 L -1 (38)
This equation just says that there is one to one correspondence between momentum of a particle and the Lorentz transformation that boosts the particle from the rest to its actual velocity. It follows that the components of momenta on the left hand side are restricted to be on shell, p 2 + m 2 = 0. As usual we can treat these components as independent adding to the lagrangian the on shell constraint p 2 + m 2 . Thus we conclude that the lagrangian (33) describes a free relativistic particle, as it should.

Returning to the starting point, eq. (33) it can be easily computed that the Lagrangian L =< Dh -1 ḣ > corresponds to the symplectic form

Ω f ree =< Dδh -1 ∧ δh >= - 1 2 [D, h -1 δh] ∧ h -1 δh
Comparing this with the first term in (32) we see that as the result of deformation instead of commutator with algebra element D = µγ 04 , in the deformed case we have to do with conjugation with group element

C = e 2πD/k = exp(2πµγ 04 /k) = cosh 2πµ k 1 + sinh 2πµ k γ 04
It can be checked that the first "free" term in (32) reduces to

Ω = k 4π cosh 2πµ k sinh 2πµ k [γ 04 , g -1 δg] ∧ g -1 δg (39)
Now we have to recall that the Chern-Simons coupling constant is related to the original coupling constant β: k/4π = 1/β ∼ 2 κ 2 , where κ = /G. Therefore the prefactor in Ω in (39) becomes, in the leading order in , 1 2 m . One can then easily check using definition of momenta (38) and the expansion (35) that (39) exactly reproduces (up to the sign and a prefactor depending on Immirzi parameter) the free particle symplectic structure. We conclude that the first term in (32) is just a sum of free particle actions.

Let us now turn to the second, interaction term, in the symplectic form (32). As an example consider the case of two particles. In this case the interaction term reads

- k 4π δK 2 K -1 2 ∧ δK 1 K -1 1 , K 1 = M -1 1 , K 2 = M -1 1 M -1 2 (40)
To calculate this let us first expand the holonomy in powers of , by making use of the Cartan decomposition of the group (35) and definition of momenta (38)

M -1 = gCg -1 = 1 + 1 κ 2 p a γ a γ 4 + O 1 2 (41) Thus - k 4π δK 2 K -1 2 ∧ δK 1 K -1 1 = - 1 κ 2 δp (2) a + δp (1) a ∧ δp (1) b γ a γ 4 γ b γ 4 = 1 κ 2 δp (2) a ∧ δp (1)a (42)
Equation ( 42) can be easily generalized to an arbitrary number of particles, and knowing the symplectic form one can readily reproduce the lagrangian for the n-particles system. It is a sum of the standard kinetic terms and the on shell constraints for each particle along with the interaction terms

L = n i=1   p (i) • d dτ q (i) + λ (i) p (i)2 + m (i)2 + 1 κ 2 p (i) • i-1 j=1 d dτ p (j)   ( 43
)
where τ is the affine parameter and λ(τ ) is the Lagrange multiplier enforcing the mass shell constraints, and • is the standard Minkowski product. Equation ( 43) is the final result of our paper. It shows that in the action of n-particles system, the actions of each particle is not deformed while the deformation arises in the form of the presence of the additional interaction terms. It is worth noticing that this deformation does not change equations of motion (because the momenta p (i) are constants of motion), while it certainly changes the Poisson brackets. In particular the brackets of positions will not vanish anymore. However, it is easy to see that by simple change of variables q (i) -→ q(i) = q (i) + 1

κ 2 i-1 j=1 d dτ p (j) (44)
the lagrangian (43) can be cast into the form of the one of the standard system of free relativistic particles.

Discussion

The results of the preceding section can be regarded surprising. Building on the experience with 3d gravity one would expect that the deformation of the final particle action is essentially guaranteed, and the recent investigations reported in [START_REF] Meusburger | Generalised Chern-Simons actions for 3d gravity and kappa-Poincare symmetry[END_REF] suggested that the deformation was to be of κ-Poincarè form [START_REF] Lukierski | New quantum Poincare algebra and k deformed field theory[END_REF]. Instead what we got was just a standard undeformed relativistic particles action. Let us therefore try to understand this result.

To get some more insight let us try to investigate how the results above change if we go beyond the the leading order in large cosmological scale expansion. Our starting point is again the formula (39)

Ω = k 4π sinh 4πµ k γ 04 g -1 δg ∧ g -1 δg (45)
Using the definition of momenta and Cartan decomposition as above one computes

Ω = k 4πm sinh 4πµ k p a γ a4 T -1 δT ∧ T -1 δT -δp a ∧ γ a4 T -1 δT (46)
with momenta restricted to be on-shell p 2 + m 2 = 0. It is now convenient to parametrize the translational part of the group as follows

T = Q 4 1 + Q a γ a γ 4 (47)
with (dimensionless) positions Q belonging to de Sitter space

Q 2 4 + Q 2 = 1, Q 2 ≡ -Q 2 0 + Q 2 (48)
Plugging this into the symplectic form above after straightforward calculations we find

Ω = k 4πm sinh 4πµ k δ Q 4 p a δQ a + 1 Q 4 p a Q a Q b δQ b (49)
It follows that the lagrangian is again up to the irrelevant prefactor just the one of a relativistic particle moving on de Sitter background: defining q a = Q a we have

L = 1 1 -q 2 / 2 p a qa + 1 3 1
1 -q 2 / 2 p a q a q b qb + λ p 2 + m 2 (50)

Notice that again we see no trace of any deformation and both the symplectic structure and the lagrangian are linear in momenta, which makes the position space commutative. However we were not able yet to calculate the second "interaction" terms of the symplectic structure beyond the leading term in large expansion, which may change the picture considerably. It is a general result of Alekseev and Malkin [START_REF] Alekseev | Symplectic structure of the moduli space of flat connection on a Riemann surface[END_REF] that such interactions can be removed by an appropriate symplectic transformation. In this transformation the Borel subroup of the original gauge group is known to play a role, and on the other hand it is directly related to κ-Poincarè algebra and κ-Minkowski space [START_REF] Freidel | Field theory on κ-Minkowski space revisited: Noether charges and breaking of Lorentz symmetry[END_REF]. One should also remember that in order to get the κ-Poincarè algebra as a contraction limit of → ∞, a nontrivial rescaling of momenta is required (see [START_REF] Lukierski | New quantum Poincare algebra and k deformed field theory[END_REF] and [START_REF] Amelino-Camelia | Quantum symmetry, the cosmological constant and Planck scale phenomenology[END_REF] for details.) Although such rescaling is well motivated mathematically, it is not clear how to justify it physically. We are going to address all these questions in the forthcoming paper.

Figure 1 :

 1 Figure 1:

Strictly speaking it is an element of the dual algebra so(4, 1) * . Here we describe it by a canonically conjugated element of the algebra.

Acknowledgment

We would like to thank Florian Girelli, Catherine Meusburger, and especially Bernd Schroers for many discussions. Thanks are also due to Lee Smolin for reading the preliminary version of the manuscript and his constant encouragement. For JKG this research was supported as a part of 2007-2010 research project N202 081 32/1844. AS is grateful to ENRAGE network for the support at Utrecht University where this work started.