text:
# Review of Vector Condensation at High Chemical Potential

## 1 Introduction

Relativistic vector condensation has been proposed and studied in different realms of theoretical physics. However the condensation mechanism and the nature of the relativistic vector mesons themselves is quite different. Linde Linde:1979pr ; Ferrer:jc ; Ferrer:ag , for example, proposed the condensation of the intermediate vector boson $`W`$ in the presence of a superdense fermionic matter while Ambjørn and Olesen Ambjorn:1989gb ; Kajantie:1998rz  investigated their condensation in presence of a high external magnetic field. Manton Manton:1979kb  and later on Hosotani Hosotani:1988bm  considered the extension of gauge theories in extra dimensions and suggested that when the extra dimensions are non simply connected the gauge fields might condense. Li in Li:2002iw  has also explored a simple effective Lagrangian and the effects of vector condensation when the vectors live in extra space dimensions111The effective Lagrangian and the condensation phenomenon in extra (non compact) space dimensions has also been studied/suggested in Sannino:2001fd ..

In Brown:kk ; Langfeld  it has also been suggested that the non gauge vectors fields such as the (quark) composite field $`\rho `$ in QCD may become light and possibly condense in a high quark matter density and/or in hot QCD. Harada and Yamawaki’s Harada:2000kb  dynamical computations within the framework of the hidden local gauge BKY  symmetry support this picture. Recently the Bose-Einstein relativistic condensation phenomenon has been suggested as physical mechanism underlying the electroweak symmetry breaking scenario Sannino:2003mt ; Sannino:2003ai

We consider another type of condensation. If vectors themselves carry some global charges we can introduce a non zero chemical potential associated to some of these charges. If the chemical potential is sufficiently high one can show that the gaps (i.e. the energy at zero momentum) of these vectors become light Sannino:2002wp ; Sannino:2001fd ; Lenaghan:2001sd  and eventually zero signaling an instability. If one applies our results to 2 color Quantum Chromo Dynamics (QCD) at non zero baryon chemical potential one predicts that the vectors made out of two quarks (in the 2 color theory the baryonic degrees of freedom are bosons) condense. Recently lattice studies for 2 color at high baryonic potential Alles:2002st  seem to support our predictions. This is the relativistic vectorial Bose-Einstein condensation phenomenon. A decrease in the gap of vectors is also suggested at high baryon chemical potential for two colors in Muroya:2002ry .

I review here the general structure of vector condensation Sannino:2002wp .

## 2 Vacuum Structure and Different Phases

We choose to consider the following general effective Lagrangian for a relativistic massive vector field in the adjoint of $`SU(2)`$ in $`3+1`$ dimensions and up to four vector fields, two derivatives and containing only intrinsic positive parity terms ARS ; DRS 222For simplicity and in view of the possible physical applications we take the vectors to belong to the adjoint representation of the $`SU(2)`$ group.:

$`L`$ $`=`$ $`{\displaystyle \frac{1}{4}}F_{\mu \nu }^aF^{a\mu \nu }+{\displaystyle \frac{m^2}{2}}A_\mu ^aA^{a\mu }+\delta ϵ^{abc}_\mu A_{a\nu }A_b^\mu A_c^\nu `$ (1)
$``$ $`{\displaystyle \frac{\lambda }{4}}\left(A_\mu ^aA^{a\mu }\right)^2+{\displaystyle \frac{\lambda ^{}}{4}}\left(A_\mu ^aA^{a\nu }\right)^2,`$
with $`F_{\mu \nu }^a=_\mu A_\nu ^a_\nu A_\mu ^a`$, $`a=1,2,3`$ and metric convention $`\eta ^{\mu \nu }=\mathrm{diag}(+,,,)`$. Here, $`\delta `$ is a real dimensionless coefficient, $`m^2`$ is the tree level mass term and $`\lambda `$ and $`\lambda ^{}`$ are positive dimensionless coefficients with $`\lambda \lambda ^{}`$ when $`\lambda ^{}0`$ or $`\lambda 0`$ when $`\lambda ^{}0`$ to insure positivity of the potential. The Lagrangian describes a self interacting $`SU(2)`$ Yang-Mills theory in the limit $`m^2=0`$, $`\lambda =\lambda ^{}>0`$ and $`\delta =\sqrt{\lambda }`$.

We set $`\delta =0`$ and the theory gains a new symmetry according to which we have always a total number of even vectors in any process Sannino:2002wp . The effect of a nonzero chemical potential associated to a given conserved charge - (say $`T^3=\tau ^3/2`$) - can be included by modifying the derivatives acting on the vector fields according to $`_\nu A_\rho _\nu A_\rho i[B_\nu ,A_\rho ]`$ with $`B_\nu =\mu \delta _{\nu 0}T^3V_\nu T^3`$ where $`V=(\mu ,\stackrel{}{0})`$. The chemical potential breaks explicitly the Lorentz transformation leaving invariant the rotational symmetry. Also the $`SU(2)`$ internal symmetry breaks to a $`U(1)`$ symmetry. If the $`\delta `$ term is absent we have an extra unbroken $`Z_2`$ symmetry which acts according to $`A_\mu ^3A_\mu ^3`$. These symmetries suggest introducing the following cylindric coordinates:

$`\varphi _\mu `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(A_\mu ^1+iA_\mu ^2),\varphi _\mu ^{}={\displaystyle \frac{1}{\sqrt{2}}}(A_\mu ^1iA_\mu ^2),\psi _\mu =A_\mu ^3,`$ (2)
on which the covariant derivative acts as follows:

$`D_\mu \varphi _\nu =\left(+iV\right)_\mu \varphi _\nu ,D_\mu \psi _\nu =_\mu \psi _\nu ,V_\nu =(\mu ,\mathrm{𝟎}).`$ (3)
The vacuum structure of the theory is explored via the variational ansatz Sannino:2002wp :

$`\psi ^\mu =0,\varphi ^\mu =\sigma \left(\begin{array}{c}0\\ 1\\ e^{i\alpha }\\ 0\end{array}\right).`$ (8)
Substituting the ansatz in the potential expression yields:

$`V=2\sigma ^4\left[(2\lambda \lambda ^{})\lambda ^{}\mathrm{cos}^2\alpha \right]+2\left(m^2\mu ^2\right)\sigma ^2.`$ (9)
The potential is positive for any value of $`\alpha `$ when $`\lambda >\lambda ^{}`$ if $`\lambda ^{}0`$ or $`\lambda >0`$ if $`\lambda ^{}<0`$. Due to our ansatz the ground state is independent of $`\delta `$. The unbroken phase occurs when $`\mu m`$ and the minimum is at $`\sigma =0`$. A possible broken phase is achieved when $`\mu >m`$ since in this case the quadratic term in $`\sigma `$ is negative. According to the value of $`\lambda ^{}`$ we distinguish three distinct phases:

### 2.1 The polar phase: $`\lambda ^{}>0`$

In this phase the minimum is for

$`\varphi ^\mu =\sigma \left(\begin{array}{c}0\\ 1\\ 1\\ 0\end{array}\right),\mathrm{with}\sigma ^2={\displaystyle \frac{1}{4}}{\displaystyle \frac{\mu ^2m^2}{\lambda \lambda ^{}}},`$ (14)
and we have the following pattern of symmetry breaking $`SO(3)\times U(1)SO(2)`$, where $`SO(3)`$ is the rotational group. We have three broken generators and 3 gapless excitations with linear dispersion relations. All of the physical states (with and without a gap) are either vectors (2-component) or scalars with respect to the unbroken $`SO(2)`$ group. The dispersion relations for the 3 gapless states can be found in Sannino:2002wp .

At $`\mu =m`$ the dispersion relations are no longer linear in the momentum. This is related to the fact that the specific part of the potential term has a partial conformal symmetry discussed first in Sannino:2001fd . Some states in the theory are curvatureless but the chemical potential present in the derivative term prevents these states to be gapless. There is a transfer of the conformal symmetry information from the potential term to the vanishing of the velocity of the gapless excitations related to the would be gapless states. This conversion is due to the linear time-derivative term induced by the presence of the chemical potential term Schafer:2001bq ; Sannino:2001fd ; Sannino:2002wp .

### 2.2 Enhanced symmetry and type II Goldstone bosons: $`\lambda ^{}=0`$

Here the potential has an enhanced $`SO(6)`$ in contrast to the $`SU(2)\times U(1)`$ for $`\lambda ^{}0`$ global symmetry which breaks to an $`SO(5)`$ with 5 broken generators. Expanding the potential around the vacuum we find 5 null curvatures Sannino:2002wp . However we have only three gapless states obtained diagonalizing the quadratic kinetic term and the potential term. Two states (a vector of $`SO(2)`$) become type II goldstone bosons while the scalar state remains type I Nielsen . This latter state is the goldstone boson related to the spontaneously broken $`U(1)`$ symmetry333According to the Nielsen and Chadha counting scheme in absence of Lorentz invariance if $`n_I`$ denotes the number of gapless excitations of type $`I`$ with linear dispersion relations (i.e. $`Ep`$) and $`n_{II}`$ the ones with quadratic dispersion relations (i.e. $`Ep^2`$) the Goldstone theorem generalizes to $`n_I+2n_{II}\mathrm{\#}\mathrm{broken}\mathrm{generators}`$.. Using the Nielsen Chadha theoremNielsen  the type II states are counted twice with respect to the number of broken generators while the linear just once recovering the number of generators broken by the vacuum. We can be more specific since we discovered Sannino:2001fd  that the velocity of the associated gapless states is proportional to the curvatures (evaluate on the minimum) of the would be goldstone bosons which is zero in the $`\lambda ^{}=0`$ limit. Again we have an efficient mechanism for communicating the information of the extra broken symmetries from the curvatures to the velocities of the already gapless excitations.

### 2.3 The apolar phase: $`\lambda ^{}<0`$

In this case the potential is minimized for:

$`\varphi ^\mu =\sigma \left(\begin{array}{c}0\\ 1\\ ı\\ 0\end{array}\right),\mathrm{with}\sigma ^2={\displaystyle \frac{1}{2}}{\displaystyle \frac{\mu ^2m^2}{2\lambda \lambda ^{}}},`$ (19)
with 3 broken generators. However the unbroken generator is a linear combination of the $`U(1)`$ and a $`SO(3)`$ generator and in Sannino:2002wp  it has been shown that only two gapless states emerges. One of the two states is a type I goldstone boson while the other is type II. The two goldstone bosons are one in the $`z`$ and the other in the $`xy`$ plane. Interestingly in this phase, due to the intrinsic complex nature of the vev, we have spontaneous $`CP`$ breaking. We summarize in Fig. 1 the phase structure in terms of the number of goldstone bosons and their type according to the values assumed by $`\lambda ^{}`$.

### 2.4 The case $`\lambda =\lambda ^{}`$: the gauge theory limit

Here the potential is:

$`V=2\sigma ^4\lambda \mathrm{sin}^2\alpha +2\left(m^2\mu ^2\right)\sigma ^2,`$ (20)
with two extrema when $`\mu >m`$, one for $`\alpha =0`$ and $`\sigma =0`$ which is an unstable point and the other for $`\alpha =\pm \pi /2`$ and $`\sigma ^2=\frac{(\mu ^2m^2)}{2\lambda }`$ corresponding to a saddle point (see the potential in Fig. 2).

At first the fact that we have no stable solutions seems unreasonable since we know that in literature we often encounter condensation of intermediate vector mesons such as the $`W`$ boson. However (except for extending the theory in higher space dimensions) in these cases one often introduces an external source. For example one adds to the theory a strong magnetic field (say in the direction $`z`$) which couples to the electromagnetically charged intermediate vector bosons $`W^+`$ and $`W^{}`$). In this case the potential is (see Ambjorn:1989gb ):

$`V=2\sigma ^4\lambda \mathrm{sin}^2\alpha +2\left(m^2eH\mathrm{sin}\alpha \right)\sigma ^2,`$ (21)
where $`e`$ is the electromagnetic charge and $`H`$ is the external electromagnetic source field. This potential has a true minimum for $`\alpha =\pi /2`$ and $`\sigma =\frac{eHm^2}{2\lambda }`$ whenever the external magnetic field satisfies the relation $`eH>m`$. We learn that the relativistic vector theory is unstable at large chemical potential whenever the non derivative vector self interactions are tuned to be identical. This is precisely the limit often used in literature when writing effective Lagrangians that in QCD describe the $`\rho `$ vector field. In principle we can still imagine to stabilize the potential in the gauge limit by adding some higher order operators. This might be the case if one introduces massive gauge bosons as in BKY . Due to the gauge limit $`\lambda ^{}=\lambda `$ is positive and assuming the higher order corrections to be small one predicts a polar phase within this model. Another solution to this instability is that the chemical potential actually does not rise above the mass of the vectors even if we increase the relative charge density. This phenomenon is similar to what happens in the case of an ideal bose gas at high chemical potential kapusta . If the strict gauge limit is taken (i.e. also the mass term is set to zero) we need to impose gauge neutrality and the analysis modifies kapusta .

Interestingly by studying the vector condensation phenomenon for strongly interacting theories on the lattice at high isospin chemical potential we can determine the best way of describing the ordinary vector self-interactions at zero chemical potential.

## 3 Physical Applications and Conclusions

We presented the phase structure of the relativistic massive vector condensation phenomenon due to a non zero chemical potential associated to some of the global conserved charges of the theory Sannino:2002wp . The phase structure is very rich. According to the value assumed by $`\lambda ^{}`$ we have three independent phases. The polar phase with $`\lambda ^{}`$ positive is characterized by a real vacuum expectation value and 3 goldstone bosons of type I. The apolar phase for $`\lambda ^{}`$ negative has a complex vector vacuum expectation value spontaneously breaking CP. In this phase we have one goldstone boson of type I and one of type II while still breaking 3 continuous symmetries. The third phase has an enhanced potential type symmetry and 3 goldstone bosons one of type I and two of type II.

We also discovered that if we force the self interaction couplings $`\lambda `$ and $`\lambda ^{}`$ to be identical, as predicted in a Yang-Mills massive theory, our ansatz for the vacuum does not lead to a stable minimum when increasing the chemical potential above the mass of the vectors. We suggest that lattice studies at high isospin chemical potential in the vector channel for QCD might be able to, indirectly, shed light on this sector of the theory at zero chemical potential. More generally the hope is that these studies might help understanding how to construct consistent theories of interacting massive higher spin fields not necessarily related to a gauge principle.

The present knowledge can be used for a number of physical phenomena of topical interest. For example in the framework of 2 color QCD at non zero baryon chemical potentialSplittorff:2002xn ; Dunne:2003ji  vector condensation has been predicted in Lenaghan:2001sd ; Sannino:2001fd . Recent lattice studies Alles:2002st  seem to support it. StudiesHands:2003dh  of the Gross-Neveu model in 2+1 dimensions with a baryon chemical potential might also shed light on the vector meson channel. The present analysis while reinforcing the scenario of vector condensation shows that we can have many different types of condensations with very distinct signatures. Other possible physical applications are discussed in Sannino:2001fd . The analysis has been extended to a general number of space dimensions Sannino:2001fd  and is useful for various scenarios related to the phenomenon of vector condensation Li:2002iw ; Moffat:2002nu .


format: text

text:
# Burgers Equation Revisited

## I Introduction

The problem of the randomly driven Burgers equation, or 1D pressureless turbulence, attracted a lot of attention in the literature in the last decade Yakhot ; Polyakov ; BMP ; GM ; Sinai ; Frisch ; KG ; Bec . Indeed, it is a tantalizing problem. On the one hand, it incorporates many features one would expect from a real 3D turbulent fluid. On the other hand, it is much simpler than the 3D turbulence of the incompressible fluid, usually studied in the framework of the Navier-Stokes equation. Its simplicity prompted some researchers to suggest that it will become the “Ising model” of turbulence. Yet many published papers later, we are still far from a complete solution to this problem. In fact, many of the aspects of the Burgers equation remain the subject of an ongoing debate. Several competing methods exist in the literature, each with its own range of applicability, which provide different, sometimes mutually exclusive answers to various aspects of the 1D turbulence.

The situation could be improved if a way existed to check various statements made about the Burgers equation. Such a check could be provided by numerical simulations. But the numerical simulations in turbulence are hard to conduct, they require large expenditure of resources, and the answers they provide are often too ambiguous to either confirm or disprove a particular theory. The situation in Burgers turbulence seems to be no different. Despite several influential papers on the numerical aspects of Burgers turbulence, those studies failed to rule out completely or confirm any of the competing theories currently on the market. The exception to this seems to be the paper of J. Bec Bec , where he confirmed the asymptotics of a certain probability distribution function (discussed below) predicted earlier in Ref. Sinai .

In this paper, I would like to point out that there exists an alternative way to study the Burgers equation numerically. The method involves minimizing the Hopf-Cole functional instead of trying to solve a complicated nonlinear differential equation. That turned out to be extremely easy to do numerically. A ten line program, which can be written in half an hour, produces, after several minutes of computer time, clean graphs of probability distribution functions. With the help of this program, I report confirmation of some of the existing theories. I also hope that the method, being so simple, can be used by other researchers to do a quick check if their proposed theory indeed does not contradict the numerical experiment.

The rest of this paper is organized as follows. In section II I formulate the problem of Burgers turbulence and summarize many of the known facts about it. In section III I show how Burgers equation must be discretized in order to facilitate doing numerics. In section IV I discuss the statistics of the velocity gradients of the Burgers equation using the methods of Lagrangian trajectories. I also discuss the “anomaly” hypothesis of Polyakov and explain its meaning within the Lagrangian approach. In section V I discuss the numerical simulations, show that they are compatible with the asymptotic results of Refs. Sinai ; Bec  and show that they are also compatible with the Polyakov hypothesis if the “anomaly” parameter $`\beta `$ is chosen to be equal to $`3/2`$. Finally in the last section VI I also show that the compatibility with the known results on the behavior of pinned charge density waves AR ; Fogler ; GC ; GC1  also forces the value $`\beta =3/2`$. I speculate that since the problem of pinned charge density waves and the asymptotics of Refs. Sinai ; Bec  relate to different properties of Burgers turbulence, there is a chance that the Polyakov’s theory with $`\beta =3/2`$ is indeed exact, at least as far as the calculation of probabilities of velocity gradients is concerned, unless a remarkable and unlikely coincidence is at play here.

## II Formulation of the Problem

The Burgers turbulence problem can be formulated in just a few lines. Take a pressureless fluid in 1D, introduce its velocity field $`u(x,t)`$ where $`x`$ is space and $`t`$ is time and write down the Navier-Stokes equation

$$_tu+u_xu\nu _x^2u=f(x,t).$$
(1)
Here $`\nu `$ is the viscosity of the fluid which is taken to be very small. $`f(x,t)`$ is a force acting on the fluid, which is taken to be random, white noise in time, and a smooth function in space

$$f(x,t)f(x,t^{})=\delta (tt^{})F(xx^{}),$$
(2)
where $`F(x)`$ is a function which is positive at $`x=0`$ and smoothly goes to zero when $`xL`$, where $`L`$ sets the force length scale. As far as boundary conditions are concerned, it is usually assumed that $`u`$ is either periodic in space (in which case $`F`$ has to be periodic as well), or the space is infinite and $`u0`$ at infinity. Note that the center of mass motion decouples in both cases, to give

$$\frac{d}{dt}𝑑xu=𝑑xf,$$
(3)
which is but the standard Brownian motion. For this reason one usually selects the force $`f`$ to be a derivative of another function $`h(x,t)`$, $`f=_xh`$. With this restriction, the average velocity vanishes $`u=0`$ and it makes sense to define the root mean square velocity as $`U_{\mathrm{rms}}=\sqrt{u^2}`$. It measures the average amplitude of the fluctuation of the velocity field. A basic assumption in the theory of turbulence is that this quantity does not depend on the viscosity $`\nu `$ as it is taken to zero. That allows to determine $`U_{\mathrm{rms}}`$ via simple hydrodynamic scaling to be $`U_{\mathrm{rms}}(LF(0))^{\frac{1}{3}}`$.

Once the equation (1) is written, one would like to calculate various averages and probabilities. A particular popular question focuses on the probability distribution function of the velocity at two different points. Define $`P(\sigma )`$ as the probability of observing a velocity gradient $`\sigma =_xu(x)`$ in a particular point in space. Similarly $`P_d(v,r)`$ is the probability of observing a difference in velocities $`v=u(x+r)u(x)`$ at two points separated by a distance $`r`$ from each other. These quantities received a lot of attention in the literature, possibly because they are related to the turbulent advection of particles suspended in the fluid, which is a very important topic in real 3D turbulence.

The following is a highly subjective historical overview of the progress made in understanding (1). As always, the choice of highlights is very personal one, and I apologize to the authors whose contributions are not mentioned here.

The first to address the behavior of $`P_d(v)`$ were V. Yakhot and A. Chekhlov Yakhot , who argued that it should not be symmetric under $`vv`$. They were also the first to study this function numerically.

Soon thereafter, A. Polyakov Polyakov  suggested an approach to calculation of these functions. He suggested an equation that $`P_d(v)`$ should satisfy. The drawback of his approach was that it was not just one equation but a family of equations parametrized by a parameter (referred to as “$`b`$-anomaly”) whose value remained undetermined. At about the same time, Bouchaud, Mezard and Parisi BMP  attacked the problem using replica trick and the Gaussian variational ansatz. The problem with their approach was that it was not well suited to computing quantities such as $`P(\sigma )`$ and $`P_d(v)`$.

In subsequent work GM  of the author with A. Migdal, it was shown that at least the asymptotics of $`P(\sigma )`$ and $`P_d(v)`$ at $`\sigma ,v+\mathrm{}`$ can be derived carefully in the saddle point approximation. It was shown that

$`P(\sigma )\mathrm{exp}\left({\displaystyle \frac{2\sigma ^3}{3G}}\right),\sigma G^{\frac{1}{3}}{\displaystyle \frac{U_{\mathrm{rms}}}{L}}`$ (4)
$`P_d(v)\mathrm{exp}\left({\displaystyle \frac{2v^3}{3r^3G}}\right),vU_{\mathrm{rms}}`$
where $`G=_x^2F(x)|_{x=0}`$. The drawback of the saddle point approach was that it allowed to determine the right hand side asymptotics only. But that asymptotics is expected to depend on whether the force $`f`$ of (1) is Gaussian. The asymptotics (4) is true only if the random force of the Burgers equation is Gaussian. But it does not have to be. For example, one could take the random force whose strength is always less than a certain value, $`|f|<f_{\mathrm{max}}`$. Then the asymptotics (4) would not be expected to hold. But if it is so, then one could legitimately ask if at least some part of the functions $`P(v)`$ and $`P_d`$ are universal, that is, they are independent on the details of the random force $`f`$, except via the force length scale $`L`$. On top of it, since any random Gaussian force leads to Eq. (4), one could ask if at least for the Gaussian force, the entire functions $`P(\sigma )`$ and $`P_d(v)`$ are universal (that is, independent of the force correlation function $`F(x)`$).

Some time later, an important paper by W. E, K. Khanin, A. Mazel and Y. Sinai Sinai  came out. It was argued in that paper that the asymptotics of $`P(\sigma )`$ at large negative $`\sigma `$ is

$$P(\sigma )\frac{1}{|\sigma |^{\frac{7}{2}}},\sigma \frac{U_{\mathrm{rms}}}{L}.$$
(5)
The arguments leading to Eq. (5) seemed to be well defined and valid independent of the details of the random force setup. Those arguments were later elaborated in a later paper by U. Frisch, J. Bec and B. Villone Frisch . In all of these papers, the asymptotics Eq. (4) was not disputed.

At the same time, these authors argued that the only universal feature of the probability distribution function $`P(\sigma )`$ is its “left tail” Eq. (5). The rest of the function was proposed to be nonuniversal, depending on the details of the random force.

These claims were later checked numerically by R. Kraichnan and T. Gotoh KG , but they were not able to unambiguously confirm or reject the power law (5). Some time later J. Bec Bec  did a very thorough numerical analysis of the “left tail” and concluded that the power law Eq. (5) is indeed correct.

Here I would like to report a numerical confirmation of the tail Eq. (5). However, I would also like to suggest that the left tail is not the only universal feature of the function $`P(\sigma )`$.

## III Discrete Burgers Equation

To study the Burgers equation numerically, its space and time have to be discretized. One way to do that would be to write derivatives as differences. That is a very crude scheme, however. Burgers equation allows a much more gentle way of discretizing its time, which is based on the Hopf-Cole transformation. In this section I am going to present the scheme as it is, and show how it produces the Burgers equation in the continuum limit.

Consider the energy function $`E_i(x)`$ where $`x`$ changes from $`0`$ to $`2\pi `$ and $`i`$ is the integer index (which will later become time). Now write down the following definition

$$E_i(x)=\underset{0y2\pi }{\mathrm{min}}\left[1\mathrm{cos}\left(xy\right)+E_{i1}(y)\right]+h_i(x).$$
(6)
Here $`h_i`$ is a random function, which I take as

$$h_i(x)=A_i\mathrm{cos}\left(x\varphi _i\right),$$
(7)
where $`A_i`$ is a random variable (for example, random Gaussian with or taking values uniformly distributed from $`0`$ to some maximum $`A_{\mathrm{max}}`$), and $`\varphi _i`$ is a random phase uniformly distributed from $`0`$ to $`2\pi `$. The sign $`\mathrm{min}`$ in Eq. (6) should be understood as a minimization over $`y`$ which varies from $`0`$ to ‘$`2\pi `$ of the expression in the square brackets. Eq. (6) plays the role of a definition of $`E_i`$ in terms of $`E_{i1}`$. As a starting point, one can take $`E_0(x)=0`$. Finally, the functions $`E_i(x)`$ will automatically be periodic in $`x`$.

The equation Eq. (6) is the discrete Burgers equation. To see that, let us take the continuum limit. In the notations of (6), the continuum limit is achieved when $`h_i1`$, $`E_i1`$. The minimization condition reads

$$\mathrm{sin}(xy)_yE_{i1}(y)=0.$$
(8)
Since $`E_i`$ is small, one can expand

$$x=y+_yE_{i1}.$$
(9)
Substituting it back to (6) and expanding $`E_{i1}(x+yx)`$ in powers of $`yx`$, one finds

$$E_i(x)E_{i1}(x)+\frac{1}{2}\left(_xE_{i1}(x)\right)^2=h_i(x).$$
(10)
The final step is replacing the discrete index $`i`$ with the continuous time $`t`$, to find

$$_tE+\frac{1}{2}\left(_xE\right)^2=h,$$
(11)
which is the so-called KPZ equation. Introducing the velocity $`u=_xE`$, one gets the Burgers equation

$$_tu+u_xu=_xh,$$
(12)
with $`_xh`$ identified as the force $`f`$. In this equation, the viscosity term $`\nu u_{xx}`$ is absent. But it is easy to see that find a global minimum in (6) is equivalent to the infinitesimal viscosity term in (12), following, for example, Feigelman Feigelman . The derivation below is somewhat technical and is given for completeness only, so it is possible to skip the equations Eqs. (13)-(15) without losing any essential information.

Introduce an alternative way of writing down (6)

$$\mathrm{exp}\left(E_i(x_i)\right)=\underset{\nu 0}{lim}\underset{j<i}{}dx_j\mathrm{exp}\left\{\frac{2}{\nu }\underset{ji}{}\left[1\mathrm{cos}\left(x_jx_{j1}\right)+h_j(x)\right]\right\}.$$
(13)
It is clear that the integral, in the limit of vanishing $`\nu `$ effectively minimizes the expression in the exponential, doing the work which is accomplished by the $`\mathrm{min}`$ sign in (6). Passing to the continuum limit, one finds the functional integral

$$\mathrm{exp}\left(E(y,T)\right)=\underset{\nu 0}{lim}𝒟x(t)\mathrm{exp}\left\{\frac{2}{\nu }_0^T𝑑t\left[\frac{1}{2}\dot{x}^2+h(x,t)\right]\right\},x(T)=y.$$
(14)
But the right hand side is none other than the propagator of the Schrödinger equation, which enables us to conclude that $`\mathrm{exp}\left(E\right)`$ satisfies the Schrödinger equation and as a consequence,

$$_tE+\frac{1}{2}\left(_xE\right)^2\nu _x^2E=h,_tu+u_xu\nu _x^2u=_xh,$$
(15)
which concludes the derivation. (In the last line of (15) I replaced the arguments of the function $`E(y,T)`$ by $`x`$ and $`t`$.)

All this suggests that the following steps must be taken to simulate the Burgers equation numerically. First, generate the random “force” $`h_i(x)`$ using (7). $`A_i`$ can be taken, for example, uniformly distributed on the interval from $`0`$ to $`A_{\mathrm{max}}`$. Then solve the relationship (6) recursively to find $`E_i(x)`$. The actual value of $`E_i`$ is not important for the numerics of the Burgers equation. What is important is the value of $`y`$ which minimizes the expression in the square brackets of (6) for each value of $`x`$. Given the relationship between $`x`$ and $`y`$, one finds for the Burgers velocity $`u_i=_xE_i(x)`$

$$u_i(x)=u_{i1}(y)+_xh_i(x).$$
(16)
This is the discrete analog of the obvious equation $`\frac{d}{dt}u(x(t),t)=f(x,t)`$ where $`\frac{d}{dt}x(t)=u(x(t),t)`$, the so-called Lagrangian trajectory (see below). Finally, for the velocity gradients $`\sigma _i(x)=_x^2E_i(x)`$ one finds

$$\sigma _i(x)=\mathrm{cos}\left(xy\right)\frac{\sigma _{i1}(y)}{\mathrm{cos}\left(xy\right)+\sigma _{i1}(y)}+_x^2h_i(x),$$
(17)
which is the discrete version of the Eq. (21) considered below. The relations Eq. (6), Eq. (16), and Eq. (17) are those one has to iterate to solve the Burgers equation numerically. To stay close to the continuum limit in time, one has to take $`A_{\mathrm{max}}1`$. And finally, to do the minimization in (6) in practice one has to discretize space as well. With 2000 spacial discretization point, reliable data can be accumulated over about 10000 steps in time, which takes about 10 minutes computer time on a Sun workstation. The result of one such calculation, with $`A_{\mathrm{max}}=0.01`$ is shown on Fig 1.

The random force chosen according to (7) fixes the function $`F`$ to be $`F(x)\mathrm{cos}(x)`$ and the characteristic length scale $`L`$ is then of the order of the size of the system. If one wants to use a shorter range force, one could take

$$h_i(x)=A_i\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{exp}\left\{\frac{\left(x2\pi n+\varphi _i\right)^2}{L^2}\right\},$$
(18)
or any other similar function, as long as it is periodic in space $`x`$.

## IV Statistics of the Velocity Gradients

The properties of the solution to Burgers equation has been discussed in many publication (see for example Frisch ). Here I list those which are relevant for the present discussion.

In the Burgers equation, the viscosity $`\nu `$ is supposed to be infinitesimally small. Yet it cannot be set to zero. The reason for that is, the solution to this equation develop shock waves, or sharp drops in the value of the velocity field. The width of this drops is of the order of $`1/\nu `$ and they become discontinuities in $`u`$ when $`\nu `$ is taken to zero. In the middle of the drop $`\nu _x^2u`$ has a finite limit as $`\nu 0`$.

To study the gradients of the velocity it is advantageous to introduce the notion of the lagrangian coordinate $`x(t)`$. It is the point which moves together with the fluid

$$\frac{dx(t)}{dt}=u(x(t)).$$
(19)
A discrete version of the lagrangian trajectory was given by Eq. (9). Once $`x(t)`$ reaches a shock wave, it becomes “trapped” in it. For the purposes of this paper, it is natural to assume that $`x(t)`$ ends as soon as it reaches the shock wave (this is precisely what Eq. (9) does). Introduce the notion of the gradient of the velocity on a trajectory defined by $`x(t)`$,

$$\sigma (t)=_xu(x(t),t).$$
(20)
As a consequence of the Burgers equation,

$$\dot{\sigma }+\sigma ^2=_xf(x(t),t).$$
(21)
Notice that the $`\nu `$ term can be neglected since, by construction, $`x(t)`$ avoids shock waves and terminates once it reaches them.

The equation Eq. (21) can be interpreted as a Langevin equation with the random force $`g(t)=_xf(x(t),t)`$. Before proceeding further, one needs to establish the statistics of $`g(t)`$. $`x(t)`$ is correlated with $`f`$ in a complicated way, and it is not a priori clear if those correlations can be neglected.

Fortunately, the discrete Burgers equation discussed in the previous section provides an answer to this question. $`g(t)`$ is none other than $`_x^2h_i(x)`$. But $`x`$ is related to $`y`$ via Eq. (9), and therefore, $`x`$ knows nothing about the random phase $`\varphi _i`$ included in the definition of $`h_i`$ in Eq. (7). Therefore, $`g(t)`$ is indeed just a white noise random force

$$g(t)g(t^{})=G\delta (tt^{})$$
(22)
where, as before, $`G=_x^2F(x)|_{x=0}(U_{\mathrm{rms}}/L)^3`$. Given a Langevin equation Eq. (21), a Fokker-Planck equation can be written down for the propability $`P(\sigma )`$

$$\frac{}{\sigma }\left(\frac{G}{2}\frac{}{\sigma }+\sigma ^2\right)P=\frac{P}{t}.$$
(23)
In the present context this equation was first written down in Bouchaud , but in that paper it was given an interpretation different from the one here. In fact, a slightly more general form of this equation will be useful below. It has the form

$$\frac{}{\sigma }\left(\frac{G}{2}\frac{}{\sigma }+\sigma ^2\right)P+\beta P=\frac{P}{t}$$
(24)
where $`\beta `$ is an arbitrary parameter. At this stage, $`\beta =0`$. However, later a version of Eq. (24) with nonzero $`\beta `$ will be useful. Notice that the large $`\sigma `$ asymptotics of $`P`$ obtained from the Eq. (24) is Eq. (4), irrespective of the value of $`\beta `$.

Eq. (21) is unstable. As soon as $`\sigma `$ becomes large negative, it quickly reaches infinity. This process is none other than the formation of a shock wave. A quick estimate of the probability $`P(\sigma )`$ in this regime is easy to do. Once $`\sigma `$ is large enough and negative, one can neglect $`g(t)`$ in Eq. (21) completely and solve the equation to find

$$\sigma (t)=\frac{1}{tt_0},t<t_0$$
(25)
The probability can be estimated to be

$$P(\sigma )𝑑t\delta (\sigma \sigma (t))\frac{1}{\sigma ^2},\sigma \frac{U_{\mathrm{rms}}}{L}.$$
(26)
This estimate was used in some publications to conclude that this should be the correct power law in the probability distribution function of the Burgers equation Eq. (5). However, to say so would be premature.

It is instructive to reproduce the power law Eq. (26) with the help of the Fokker-Planck equation. To do that, observe that one has to solve this equation with an absorbing boundary conditions at $`\sigma \mathrm{}`$, corresponding to the trajectories disappearing when shock waves are formed. The standard substitution

$$P=\mathrm{exp}\left(\frac{\sigma ^3}{3G}\right)\mathrm{\Psi }(\sigma ).$$
(27)
leads to the Schrödnger equation of a quantum mechanical particle moving in a potential $`U(\sigma )=\sigma ^4/(2G)(\beta +1)\sigma `$

$$\frac{G}{2}\frac{^2\mathrm{\Psi }}{\sigma ^2}+\left(\frac{\sigma ^4}{2G}\left(\beta +1\right)\sigma \right)\mathrm{\Psi }=\frac{\mathrm{\Psi }}{t}.$$
(28)
The solution to this equation reads

$$\mathrm{\Psi }(\sigma )=\underset{n=0}{\overset{\mathrm{}}{}}C_n\mathrm{exp}\left(E_nt\right)\mathrm{\Psi }_n(\sigma ),$$
(29)
where $`C_n`$ are arbitrary coefficients, $`E_n`$ are the eigenstates, and $`\mathrm{\Psi }_n`$ are the eigenfunctions of the Schrödinger equation satisfying

$$\frac{G}{2}\frac{^2\mathrm{\Psi }_n}{\sigma ^2}+\left(\frac{\sigma ^4}{2G}\left(\beta +1\right)\sigma \right)\mathrm{\Psi }_n=E_n\mathrm{\Psi }_n.$$
(30)
At large times $`t(E_1E_0)^1`$, only the ground state survives

$`P(\sigma ,t)\mathrm{exp}\left(E_0t\right)P_0(\sigma ),t\left(E_1E_0\right)^1,`$ (31)
$`P_0(\sigma )=\mathrm{exp}\left({\displaystyle \frac{\sigma ^3}{3G}}\right)\mathrm{\Psi }_0(\sigma ).`$ (32)
A detailed analysis of Eqs. (28)-(31) at arbitrary $`\beta `$ can be found in the paper by S. Boldyrev Boldyrev . In the case of interest, $`\beta =0`$. Ref. Boldyrev  shows that in this case $`E_0>0`$, and it confirms the asymptotics Eq. (26) for $`P_0(\sigma )`$. The asymptotics at large positive $`\sigma `$ coincides with Eq. (4). The fact that $`E_0>0`$ is not surprising. It means that the total probability $`𝑑\sigma P(\sigma )`$ dicreases with time. This is due to elimination of trajectories which ended on shock waves.

However, this cannot be the correct probability of observing a velocity gradient $`\sigma `$ at a given point in space. Indeed, introduce the density of trajectories $`\rho (x,t)`$. It satisfies the continuity equation

$$_t\rho +_x\left(u\rho \right)=0.$$
(33)
In particular, the density around the given Lagrangian trajectory $`x(t)`$ is given by $`\rho _L(t)=\rho (x(t),t)`$ and it satisfies

$$\frac{d}{dt}\rho _L+\sigma \rho _L=0$$
(34)
It is clear that the probability that a particular trajectory happens to go through an observation point where a velocity gradient is measured is inversely proportional to the density of trajectories at this point. Thus a more physically relevant probability distribution function is defined as

$$P(\sigma )=\frac{1}{\rho _L}\delta (\sigma \sigma (t)).$$
(35)
This function satisfies a Fokker-Planck equation Eq. (24), but with $`\beta =1`$. In this case, it is possible to show that $`E_0<0`$ and

$$P(\sigma )\frac{1}{|\sigma |^3},\sigma \frac{U_{\mathrm{rms}}}{L}.$$
(36)
This asymptotics was proposed in the paper KG . But it cannot be the correct asymptotics either. Indeed, while the above analysis correctly removes the trajectories which end at the points where a new shock wave is formed, it fails to take into account those which fall into “mature” shock waves, the ones which formed in the past. That is why $`E_0<0`$ and the probability grows with time. That comes from overcounting the trajectories.

Unfortunately, the Eq. (21) does not seem to contain enough information to determine if there is a shock wave nearby. If that is indeed the case, then the full Burgers equation Eq. (1) has to be solved. That is the core of the problem of Burgers turbulence, and the reason the solution to it seems so close and yet so hard to get. The full Burgers equation Eq. (1) leads to a full fledged 2D quantum field theory (via Martin-Siggia-Rose formalism, for example). That theory is presumably very hard, if not impossible, to solve. The “reduced” equation Eq. (21), on the other hand, leads to a quantum mechanics encoded in Eq. (28). The solution to quantum mechanics problems are clearly within reach. The question is, therefore, if it is possible to extract all the relevant information about the Burgers equation from the Lagrangian equation Eq. (21). At this stage of the present analysis, it does not seem possible.

To move further, I will make the following leap of faith. The equation Eq. (24) was studied in Ref. Boldyrev  for arbitrary $`\beta `$. It was found that while the large positive $`\sigma `$ asymptotics of $`P`$ is always given by Eq. (4), the other asymptotics is given by

$$P(\sigma )\frac{1}{\sigma ^{\beta +2}},\sigma \frac{U_{\mathrm{rms}}}{L}.$$
(37)
On the other hand, it is known from the analysis of Ref. Sinai  and Frisch  that the correct power of this asymptotics is $`7/2`$. To match this asymptotics, I choose $`\beta =3/2`$. Then I propose to solve the equation Eq. (24) and compare what it gives with the numerics. The results of the comparison is shown in next section.

The equation Eq. (24) with an arbitrary $`\beta `$ was proposed in Polyakov  (but for $`P_d(v)`$ as opposed to $`P(\sigma )`$) where $`\beta `$ was referred to as $`b`$-anomaly. But whether the present discussion has anything to do with the methods of Ref. Polyakov  is unclear. To derive Eq. (24) with $`\beta =3/2`$, one needs to show that the rate at which the trajectories fall into mature shock waves is proportional to the extra $`\frac{1}{2}\sigma P`$ term generated in the Fokker-Planck equation Eq. (24). How to do that is not known to me. Next, the Eq. (24) cannot be exact. After all, it clearly breaks down even if the force is not Gaussian. However, it may be correct in some mean field sense, for moderate values of $`\sigma `$. Finally, I’d like to note that for this value of $`\beta `$, $`E_0<0`$ and the probability still grows with time. That must be because this description gives the relative probability of observing one value of $`\sigma `$ over another, not the absolute probability, and the overall normalization must be enforced by $`PP/\mathrm{\Pi }`$ where $`\mathrm{\Pi }=𝑑\sigma P(\sigma )`$.

## V Numerical Simulations

With the setup descibed in the preceeding sections, it is very easy to calculate $`P(\sigma )`$ numerically. Unlike some of the previous work on this subject, I did not try to compute $`P(\sigma )`$ from the data. Instead, I computed the integrated probability density

$$N(\sigma )=_{\mathrm{}}^\sigma 𝑑\mu P(\mu ).$$
(38)
Obviously its left tail is expected to go as $`N(\sigma )1/|\sigma |^{\frac{5}{2}}`$. To plot it numerically, it is enough to collect the data into an array, to sort it in the increasing order, and plot the position of an array entry as a function of its value. Fig 2. depicts $`N(\sigma )`$ in the case when the force was chosen to be as in Eq. (7), $`A_{\mathrm{max}}=0.01`$, and $`U_{\mathrm{rms}}0.02`$.

Fig 3. shows the left tail in the log-log format. For comparison, a straight line with the slope $`5/2`$ is plotted.

The agreement is striking. One could object that the power law in this graph does not extend over a sufficiently long range of $`\sigma `$. To extend it further, one needs to do numerics at weaker force for longer periods of time (to keep close to the continuum limit even at large velocity fluctuations). It was not the task of the present work to do large scale numerical simulations, and I believe the agreement I have here is good enough (especially in view of the work reported in Ref. Bec ).

Fig 4. shows the same probability function, but plotted together with the (numerical) solution to the equation Eq. (24) with $`\beta `$ chosen to be $`\beta =3/2`$.

To arrive at this picture, one parameter - the units of measure of $`\sigma `$ \- had to be adjusted to make the graphs overlap as much as possible. The agreement is not perfect. Yet the curves are close enough to suggest that perhaps it is not coincidential either. One interesting number is the probability $`p`$ that the velocity gradient is negative. It is especially interesting, because it is independent of the adjustments of the units of $`\sigma `$. Numerically $`p0.21\pm 0.01`$. The solution to the Eq. (24) gives about $`0.18`$. According to some claims in the literature, this number is not universal. However, the value of $`p`$ is reproducible whether I do the numerics with the long ranged force Eq. (7) or shorter range force Eq. (18) with various values of $`L`$. So numerically $`p`$ appears to be universal, or perhaps approximately universal.

To conclude, the numerical data does not confirm that the Eq. (24) with the parameter $`\beta =3/2`$ is the exact solution to the problem. But it may be an approximate, perhaps in some mean field sense, solution to the problem, which reproduces the actual probability density rather close. I believe this question deserves further investigations.

## VI Pinned Charge Density Waves

In the papers of the author with J.T. Chalker GC ; GC1  a remarkable correspondence was established between the randomly driven Burgers equation and the physics of pinned charge density waves. It would go beyond this paper to discuss charge density waves, but mathematically the correspondence can be formulated in the following way. Consider the functional

$$[x(t)]=_0^T𝑑t\left[\frac{1}{2}\dot{x}^2+h(x,t)\right],$$
(39)
where $`h(x,t)`$ is the same random function as the one in Eq. (13). Let us find a function $`x_0(t)`$ which, when substituted for $`x(t)`$ in Eq. (39), gives an absolute minimum for the functional $``$. Obviously it satisfies the minimization equation

$$\frac{d^2x_0}{dt^2}+\frac{h(x,t)}{x}|_{x=x_0(t)}=0.$$
(40)
Now consider the energy cost of $`x(t)`$ deviating from $`x_0(t)`$. Writing $`x(t)x_0(t)+\psi (t)`$ we find the equation for the normal modes of oscillations about the absolute minimum, given by

$$\left[\frac{d^2}{dt^2}+\frac{^2h(x,t)}{x^2}|_{x=x_0(t)}\right]\psi =\omega ^2\psi .$$
(41)
With the problem thus set up, one needs to calculate the average number of modes with frequencies less than a given frequency $`\omega `$. Such a function is denoted as $`N(\omega )`$.

This problem, which was first formulated in the context of charge density waves in Ref. FL , received a substantial amount of attention in the literature (see Refs. Feigelman ; Gia ; AR ; Fogler ; GC ; GC1 . After initial disagreements, a consensus was built which gave for the function $`N(\omega )`$ the value

$$N(\omega )\omega ^5,$$
(42)
for sufficiently small $`\omega `$. This was first proposed in Ref. AR  and then the derivation was improved in subsequent publications until Ref. GC  derived Eq. (42) in a systematic way.

One of the results of Refs. GC ; GC1  relates the calculation of $`N(\omega )`$ to the following specific question in Burgers turbulence. Consider a Burgers fluid Eq. (1) driven by a random force $`f(x,t)=_xh(x,t)`$. Consider a Lagrangian trajectory $`x_0(t)`$ which moves with a fluid and never gets absorbed by a shock wave for times $`t<T`$. Consider a velocity gradient $`\sigma (t)`$ defined in Eq. (20). Let us find time intervals $`t_1<t<t_2`$ such that $`\sigma (t)<0`$ within those intervals. Now let us calculate the probability that

$$_{t_1}^{t_2}𝑑t\sigma (t)<\mathrm{log}(\omega ),$$
(43)
for some small value of $`\omega `$. Then the function $`N(\omega )`$ of the charge density wave problem coincides with this probability.

Notice that knowing the tails of the probability distribution $`P(\sigma )`$ would not help calculating $`N(\omega )`$. This is because the typical functions $`\sigma (t)`$ which satisfy Eq. (43) are not those for which $`\sigma `$ is very large negative, but rather moderate negative $`\sigma (t)`$ which however persist over long time intervals. Therefore, the results of W. E at al Sinai  are useless if we were to calculate $`N(\omega )`$ with the help of Eq. (43).

Let us however use the methods of this paper, in particular Eq. (28), to calculate $`N(\omega )`$. To do so, we need to write down a Feynman path integral formalism which corresponds to Eq. (28) and then use the Lagrange multiplier method described in Ref. GC1 , section VI C, subsection 1, corrected however for the presence of a parameter $`\beta `$ in Eq. (28). Then we find

$$N(\omega )\omega ^{2\beta +2}.$$
(44)
Now we know the exact answer to the problem, Eq. (42). If the methods described in this paper were to reproduce the exact answer, we have to choose $`\beta =3/2`$. However, this is the same value of $`\beta `$ required to reproduce the asymptotics of W. E et al!

Two possible explanations of this are possible. First is, perhaps the Eq. (28) with $`\beta =3/2`$ indeed reproduces the right asymptotics of W. E et al and for some reason also reproduces the function $`N(\omega )`$ correctly, but this equation is still generally not correct and the probability distribution derived with its help and shown on Fig. 4 is not correct either. Second is, the Eq. (28) is indeed correct and gives the correct solution to the Burgers problem.

The first explanation seems far fetched, as a remarkable coincidence must be at play to give rise to it. The second explanation sounds much more likely. And yet, no derivation of Eq. (28) is known at this point.

## VII Conclusion

In this paper, I demonstrate that it is easy to obtain reliable numerical data on the behavior of the 1D Burgers equation. I show that the data confirms the $`7/2`$ tail of the probability distribution suggested in the literature. The data also suggests that the entire probability distribution function could be universal, determined by the solution to the equation Eq. (24) with the “anomaly” parameter $`\beta =3/2`$. I also discuss that it is possible to reproduce the known solution to the pinned charge density wave problem by applying the techniques discussed here and choosing $`\beta =3/2`$, which gives extra weight to these methods.

## VIII Acknowledgements

I am grateful to J.T. Chalker for many discussions. This paper is based on the results obtained in part in the course of work on paper GC . The work presented here was supported in part by EPSRC through its Advanced Research Fellowship programme.


text:
# Cage diffusion in liquid mercury

## I Introduction

The collective microscopic dynamics of dense hard-sphere fluids can be understood on the basis of cage diffusion. In this process, the particles with hard sphere diameter $`\sigma _{hs}`$, while all moving diffusively, find themselves locked up in a cage formed by their nearest neighbors. As a consequence, the coherent (collective) dynamic structure factor, $`S_c(q,E)`$ of the system becomes, as a function of wave number $`q`$ and energy $`E`$, very narrow in energy for wave numbers near $`q\sigma _{hs}2\pi `$, i.e., near the peak in the static structure factor $`S(q)`$. This narrowing is referred to as de Gennes narrowing of $`S_c(q,E)`$. It is believed that cage diffusion plays an important part in the dynamics of real fluids, such as noble gas fluidspeter ; cage , concentrated colloidal suspensionspusey ; schepper  and liquid metalsliviaprl ; livia .
Over the past decennia, the microscopic dynamics of liquid metals have been the subject of intensive study by neutron scattering and X-ray techniquesru ; na ; na2 ; na3 ; cs ; li2 ; li3 ; k ; k2 ; hg ; ge , and molecular dynamics simulationslivia ; li ; na4 . Liquid metals pose a particularly interesting problem because of the long range of the interatomic potential. The (extended) hydrodynamics modes have been investigated in these systems, and the influence of the mode coupling mechanism on these modesna2 ; na3 , leading to long-time tails in the correlation functions. In here, we present results on the fast (in time) decay mechanism in liquid mercury pertinent to cage diffusion, as first investigated by Bove et al.liviaprl .
Recent neutron scattering experiments by Bove et al.liviaprl  on liquid mercury at room temperature showed that the incoherent part of the scattering function (i.e., the dynamics of individual Hg atoms) was characterized by two time scales. The authors attributed the presence of more than one time scale to the process of cage diffusion. These two time scales were also present in Molecular Dynamics (MD) simulations carried out by the same grouplivia . However, their neutron scattering setup did not allow for a detailed investigation of the incoherent dynamic response, as it was setup primarily to investigate $`q`$-dependence of the sound modes outside the hydrodynamic regime.
We report inelastic neutron scattering experiments on Hg at room temperature (number density 0.0408 Å-3) at a higher energy resolution, allowing us to investigate whether cage diffusion is indeed present in Mercury, and determine the origin of the two time scales found previouslylivia . Our results show that cage diffusion is the dominant dynamic process in Hg. We find that the main contribution to the incoherent cross section is determined by a self-diffusion process. The second time scale is manifest in the neutron scattering response as a broad quasi-elastic line with an energy halfwidth of 2 meV. The mode associated with this second time scale disappears from the spectra for $`q>`$ 1 Å-1, however, its intensity at the smallest $`q`$-values is an order of magnitude larger than what could be expected from cage diffusion. We show that this mode cannot be viewed as an extension of the hydrodynamic Rayleigh mode to finite $`q`$-values. Instead, we propose that the intensity of the cage diffusion mode in the neutron scattering spectra at small wave numbers has been amplified through a valence fluctuation mechanism. Other than the presence of this broad quasi-elastic contribution, we find that the response of liquid mercury can be described in terms similar to a range of liquids investigated previously, albeit that Hg has a very large anomalous sound dispersion reflecting the presence of the electron gasliviaprl .

## II Experiments and Data Analysis

The inelastic neutron scattering experiments were performed on natural mercury, using the time-of-flight spectrometer MARI at the Rutherford Appleton Laboratories. Mercury scatters neutrons both coherently and incoherently, yielding information on both the collective and individual dynamics of the Hg atoms, respectively. Unfortunately, natural-enrichment mercury also contains neutron absorbing isotopes, thus rendering the data correction procedure more cumbersome than usual.
In order to minimize neutron absorption, we used a thin disk container of 0.08 cm nominal internal thickness, and 11.2 cm radius. Niobium disks of 0.12 cm thickness were used for the walls of the container. Two sets of measurements were taken on MARI with different incident neutron energies, 20 meV and 40 meV, with the sample cell oriented perpendicular to the beam, and at an angle of 45o, respectively. For each incident energy, we collected sample spectra, empty container spectra on a separate container of identical dimensions, vanadium spectra on an vanadium disk of 0.084 cm thickness in the same orientation as the sample, and we measured the neutron transmission ratio for each of these measurements. The latter allowed us to determine the actual thickness of the mercury in the beam by comparing it to the empty container transmission. We found the thickness to be 0.0816 $`\pm `$ 0.0004 cm, and the orientation for the 40 meV measurements to be 44o.8 $`\pm `$ 0o.6. In our correction procedure, we use these measured values for calculating the attenuation of the sample. All measurements were taken under ambient conditions.
The neutron scattering experiments were carried out in 1995 and were originally setup to measure the sound dispersion in Hg. Unfortunately, we underestimated the propagation velocity of the extended sound modes, and we failed to observe a clear signal associated with the anomalous sound dispersionliviaprl  as it mostly lay outside of our kinematic scattering window. However, we could determine the halfwidths of the coherent and incoherent part of the scattering function for wave numbers $`q`$ with 0.3 $`<q<`$ 7 Å-1, the results for which we report in this paper.
For the data reduction procedure, we first subtracted the time-independent background, obtained from an empty spectrometer run, and corrected the data for the detector efficiencies, as obtained from a vanadium reference sample. Then, we multiploied the empty container spectra by multiplying them with the energy and angle dependent sample attenuation factor for Hgsears . These sample attenuation factors are typically of the order of 5-20 for our highly-absorbing Hg sample. Next, the attenuated container spectra were subtracted from the sample plus container spectra, and the results were divided by the calculatedsears  attenuation factor of Hg. Finally, the results were corrected for multiple scattering, a 5-10% effect, by carrying out a numerical integration over the sample volume using the measuredubaldo  static structure factor as input and assuming a Lorentzian lineshape for the dynamic structure factor with halfwidth given by the mercury self-diffusion coefficient (0.1436 Å2/ps at 283 Klobo ).
As a check on our data reduction procedure, we show in Fig.  the measured 20 meV data, integrated over energy($`|E|`$ $`<`$ 4 meV) as a function of scattering angle after subtraction of container scattering, but before the sample attenuation correction. In this figure, the dip at 90 degrees is due to the large absorption cross section of Hg and the flat plate design of the sample cell, so that the sample effectively shields itself at these scattering angles. The width of the dip represents Cd shielding on the side of the container. We compare the energy integrated results to the calculated total (single plus double) scattering obtained using the published static structure factorubaldo  and diffusion constantslobo ; handbook , multiplied with the calculated attenuation factor. It is clear from this figure that the agreement is good at least to within 10%, allowing for an accurate determination of the halfwidths of the spectra. For our analysis, we only use the spectra with $`\varphi `$ $`<`$ 800 for the 20 meV dataset, and $`\varphi `$ $`<`$ 1200 for the 40 meV dataset. We show the fully corrected data for both incident neutron energies as a contour plot in Fig. . Close inspection of this figure shows the presence of a spurion at $`E`$ = 5 meV energy transfer for the 20 meV data set, and at $`E`$ = 10 meV for the 40 meV data set. This spurion corresponds to neutrons that are scattered elastically, but have traveled an additional 60 cm. We believe that this corresponds to neutrons that were scattered by the sample, and subsequently scattered by the flange holding the sample in position. In the 1995 MARI setup, this flange was insufficiently shielded. For our data analysis, we ignored the region around the spurion. In addition, the corrected data show small remnants due to the Bragg peaks of the Nb container (such as at $`q`$ = 4.7Å-1 in Fig. ). Affected $`q`$-values were also ignored in our data analysis.
In order to determine the coherent and incoherent halfwidths from the neutron scattering spectra, we fit the data to a number of Lorentzian lines in a procedure similar to the one reported by Bove et al.liviaprl . We find that a two Lorentzian line description gives good agreement ($`\chi ^2`$ 1) with our measurements $`q<`$ 1.9 Å-1. We therefore fit these spectra to

$$S_{tot}(q,E)=\frac{sE}{1e^{E/k_BT}}[\frac{\sigma _i}{\sigma _c}\frac{\mathrm{\Gamma }_i(q)}{E^2+\mathrm{\Gamma }_i^2(q)}+\frac{\chi (q)\mathrm{\Gamma }_c(q)}{E^2+\mathrm{\Gamma }_c^2(q)}].$$
(1)
In here, $`k_B`$ is Boltzmann’s constant, $`\sigma _i`$ and $`\sigma _c`$ are the incoherent and coherent cross sections for Hg ($`\sigma _i`$/$`\sigma _c`$= 0.324), $`\mathrm{\Gamma }_i(q)`$ is the incoherent halfwidth at half maximum (HWHM), $`\mathrm{\Gamma }_c(q)`$ is the coherent halfwidth, and $`\chi (q)`$ is the static susceptibility. At room temperature, $`\chi (q)`$ is almost identical to the static structure factor $`S(q)`$. The prefactor $`s`$ takes the minor difference between $`\chi (q)`$ and $`S(q)`$ into account, as well as any remaining scale factor due to the attenuation correction. The fit function was convoluted with the experimental resolution function obtained from the vanadium measurements, which is well described by an Ikeda-Carpenter form. The full width at half maximum in energy for neutrons scattered elastically for the 20 meV (40 meV) measurement is 0.4 meV (0.8 meV), independent of scattering angle (cf. Fig. ).
For $`q>`$ 1.9 Å-1, we find that we obtain a better fit if we describe the coherent part of the scattering by multiple Lorentzian lines. We find that two Lorentzian lines (damped harmonic oscillatorwouter ) give a good description at most $`q`$-values, but that around the first peak of the static structure factor the quality of the fit can be improved further by fitting to three Lorentzian lines (visco-elastic modellovesey ). The (coherent part of the) damped harmonic oscillator model has one free parameter, the visco-elastic model has two, and both models satisfy the f-sum rule.
Note that we only use the fit function to reliably determine the halfwidths of the coherent and incoherent parts of the spectra, corrected for the spectrometer resolution function. We do not address the individual fit parameters as they come out of the visco-elastic or harmonic oscillator models; we find that the fit parameters in the visco-elastic model are highly correlated due to the presence of both coherent and incoherent scattering. However, Eq.  and modifications thereof give an overall good description of the data, and the derived $`\mathrm{\Gamma }`$’s should be an accurate representation of the halfwidths.
Examples of the fits are shown in Figs.  and . At low $`q`$, the spectra are clearly a combination of a narrow contribution, and a broad contribution. We verified that this broad contribution was not an artifact of our correction procedure by noting that it was completely absent in the vanadium reference spectra (see Fig. ), nor could it be attributed to container scattering. At low $`q`$, the incoherent line width is so small that the incoherent scattering is identical to the (asymmetric) MARI resolution function. It is only possible to determine $`\mathrm{\Gamma }_i(q)`$ for $`q>1.5`$Å-1 where the incoherent halfwidth becomes comparable to the width of the energy resolution function. We fitted $`\mathrm{\Gamma }_i(q)`$ for $`1.5<q<`$ 2.1 Å-1, and for $`q`$ 3.3 Å-1. In these $`q`$ regions the coherent and incoherent line widths are significantly different from eachother and allow for a reliable determiniation of the incoherent line width. We find that the incoherent line width is determined by the constant of self-diffusion $`D_s`$, with $`D_s`$ = 0.143 $`\pm `$ 0.007 Å2/ps. This value is slightly smaller than the expected value at room temperature (0.1436 at T= 283 K and 0.159 at T= 298 Klobo ; handbook ). We did not have a temperature sensor on the sample while it was positioned inside the evacuated MARI spectrometer, nor did we log the temperature of the ISIS experimental hall during our October experiment. In all, our value seems to be reasonable, albeit on the low side of what is expectedlobo ; handbook .
In order to determine the coherent line widths, we fixed the halfwidth $`\mathrm{\Gamma }_i(q)`$ of the incoherent part at $`D_sq^2`$ with $`D_s`$ = 0.143 Å2/ps, corrected for the expected decrease with $`q`$ as calculated from the Enskog theoryenskog  (see Fig. \[5b\]). This decrease can be calculated without ambiguity, and represents the slow approach to ideal gas behavior at very large $`q`$. Having fixed the incoherent line widths, we then fitted the coherent part of the scattering to one, two or three Lorentzians. When we fitted to more than one Lorentzian, we fixed $`\chi (q)`$ to the value measured for $`S(q)`$ by neutron diffraction experimentsubaldo . From the fit results, we then determined the coherent line width, the results for which we show in Fig. , together with the incoherent line width and the static structure factorubaldo  for comparison.

## III Discussion

Inspecting Fig. \[5a\], one observes that for $`q>`$ 2 Å-1 the coherent halfwidth oscillates around the incoherent halfwidth, given by the coefficient for self-diffusion. We emphasize this oscillatory behavior at the higher momentum transfers by dividing the halfwidth by $`q^2`$ (see Fig. \[5b\]). From this figure we conclude that the oscillations of the coherent halfwidth are in phase with the oscillations of the static structure factor (see Fig. \[5c\]). In other words, whenever the static structure factor reaches a maximum (minimum), the coherent halfwidth reaches a minimum (maximum). This is what one expectscage  from cage diffusion: it should take longer to diffuse out of the cage formed by its nearest neighbors for wave numbers corresponding to a peak in the structure factor.
In Fig.  we also plot the predictions of the Enskog theoryenskog  for hard spheres. Comparing the properties of a liquid to those of an equivalent hard-sphere fluid greatly facilitates the interpretation of experimental results, and it has been applied successfully to a host of liquids and gaseshost . The one problem in this comparison lies in determining an equivalent hard-sphere diameter $`\sigma _{hs}`$, something which cannot be done unambiguously. For instance, comparing the height of the first maximum of $`S(q)`$ to that of the equivalent hard-sphere fluid yields $`\sigma _{hs}`$= 2.73 Å, but comparing the position of this first maximum to that of the hard-sphere fluid gives $`\sigma _{hs}`$= 2.89 Å, while comparing the density of the liquid at the melting point to the density at which a hard-sphere fluid meltsvanloef  results in $`\sigma _{hs}`$= 2.81 Å. Somewhat arbitrarily, we have opted for the latter method. This choice has the advantage that the calculated incoherent halfwidth from the Enskog theoryenskog  coincides with our experimentally determined self-diffusion coefficient. The results for the coherent and incoherent halfwidths for the hard-sphere fluid are shown in Fig. \[5a\] and \[5b\]. The comparison with the static structure factor is shown in Fig. \[5c\].
The agreement in the phase of the oscillations between mercury and the hard-sphere fluid is rather good (dashed-dotted curves in Fig. ), but the amplitude of the oscillations is somewhat off, in particular at intermediate $`q`$-values. The situation improves when we use the measured $`S(q)`$ instead of the calculated hard-sphere $`S(q)`$ (solid curves in Fig. ). In all, from the qualitative agreement between the liquid mercury data and a hard-sphere fluid that is known to exhibit cage diffusioncage , combined with the observation that the oscillations of the coherent halfwidth are in phase with those of the static structure factor, we conclude that cage diffusion is the dominant mechanism determining the dynamic response of liquid mercury at room temperature.
Next, we focus on the low-$`q`$ region ($`q<1`$Å-1) where Bove et al. have performed their neutron scattering experiments and MD computer simulationsliviaprl ; livia . Bove et al. observed in their neutron scattering experiments the existence of a mode of characteristic half width at halh maximum (HWHM) of 2 meV, and an additional contribution to the single-particle dynamics (other than simple diffusion) in their MD simulations with a characteristic time scale of 0.3 ps. Since a time scale of $``$0.3 ps would show up as an energy width of $``$ 2 meV in neutron scattering experiments, the authors suggested that the neutron scattering mode and the MD mode were the same, and identified this mode as being related to cage diffusion. We agree with the identification of the MD mode as being related to cage diffusion, but we argue that the 2 meV mode in the neutron scattering experiments cannot be fully explained by a straightforward cage diffusion mechanism.
In Fig.  we show our experimental data for $`q`$= 0.45 Å-1. Data at $`q`$-values close to this value show a similar pattern: a resolution limited central peak, and a broader contribution centered around zero energy transfer with a characteristic HWHM of 2 meV for all $`q`$ with $`q<`$ 1 Å-1 (see Fig. \[5a\]). In addition, the absolute values of the intensities for the broad contribution agree well with those measured by Bove et al.liviaprl . In our fit, we had tentatively attributed this broad mode at small $`q`$-values to the coherent part of the scattering (see Eq. ), however, we argue in the following that this is unlikely. There are two obvious candidates for the broad mode: the aforementioned cage diffusion mode seen in MD simulations, and the extension of the hydrodynamic Rayleigh mode (heat diffusion) to larger momentum transfers. We discuss both candidates.
First, it is expected that the neutron scattering spectra will show the equivalent of the hydrodynamic Rayleigh mode, a coherent process associated with heat diffusion. This mode is commonly referred to as the extended heat mode, hence our initial assignment of this broad mode to the coherent part of the scattering. In hydrodynamics, the expected intensity associated with the Rayleigh mode is $`S(0)`$($`\gamma `$-1)/$`\gamma `$, with the static structure factor $`S(0)`$ (= 0.005) at $`q`$ = 0 given by the compressibility, while $`\gamma `$ is the ratio of the specific heats (= 1.14 at T=283 K). Thus, the expected intensity of this heat mode (0.0006) compared to the expected intensity of the incoherent scattering ($`\sigma _i`$/$`\sigma _c`$= 0.324) is completely negligible in the hydrodynamic region. With increasing $`q`$, the heat-mode intensity will increase predominantly due to the increase in $`S(q)`$ and to a lesser extent due to a change in $`\gamma (q)`$. The halfwidth associated with the extended heat mode is shown as the dashed-dotted line in Fig. \[5a\] calculated using the Enskog theoryenskog . It is clear from this figure that for $`q<`$ 1 Å-1, that the measured halfwidth reaches a constant value, instead of decreasing $`q^2`$ as expected for the extended heat mode. Moreover, the intensity of the 2 meV mode actually increases with decreasing $`q`$, an observation which also directly follows from Fig. . Thus, the constant halfwidth of the 2 meV mode, its intensity ($``$ 0.07 at $`q`$ = 0.45Å-1) being far greater than what could be expected of the extended heat mode, combined with its increase in intensity with decreasing $`q`$ instead of the expected decrease, all lead us to conclude that the 2 meV mode is not an extended heat mode for $`q<`$-1.
Second, if a particle is locked up in a cage, it will rattle around on short time scales before it escapes. Typical time scales correspond to how long it will take a particle to cover the distance to its neighbor (this distance is typically of the order of $`\sigma _{hs}`$/10) with an average thermal velocity $`\sqrt{k_BT/m}`$, $`m`$ being the mass of a Hg atom. For mercury, one can expect characteristic times very similar to those found in the MD simulationslivia . However, the influence of this process to the decay of the self-correlation function is very small at small $`q`$-values. The reason for this is that a particle just rattles around in its cage, never moving very far from its original position, so when it is probed with a long wavelength (small $`q`$) density fluctuation (such as a neutron scattering with small momentum transfer), only a very small effect is to be expected. The effect will be of the order 1-cos($`q\mathrm{\Delta }r`$), with $`\mathrm{\Delta }r`$ the distance between the surfaces of neighboring hard-spheres. For $`q`$= 0.3 Å-1 and $`\mathrm{\Delta }r`$= 0.3 Å, the expected effect compared to the standard self-diffusion mechanism will be less than 0.5%. In the MD simulationslivia , the observed effect was 0.4$`\%`$. Based on these numbers, it seems entirely justifiable to ascribe the additional mode observed in MD simulations to cage diffusion. Conversely, the expected intensity in neutron scattering experiments should then also be 0.4% of the intensity associated with the standard self-diffusion mode. However, when we compare the intensity of the broad mode to the resolution limited self-diffusion mode, we find a ratio in excess of 20 % for the lowest $`q`$-values. This observation is at odds with the expected intensity ratio of $`<`$ 0.5%, and we therefore rule out a straighforward cage diffusion mechanism as the dominant constituent of the broad mode at the lower $`q`$-values. We do not dispute the existence of a mode related to cage diffusion in the neutron scattering spectra, we just calculate that its predicted intensity is too small to account for the observed intensity.
We can only speculate as to the origin of the 2 meV mode in the neutron scattering spectra. First, this mode (see Fig. ) is not an artifact of our data correction procedure as our time-of-flight results are quantitatively similar to the neutron scattering experimentsliviaprl  performed on a triple-axis spectrometer. Second, this mode corresponds to excess intensity compared to what is expected for coherent plus incoherent scattering by mercury (see Fig. ). One possible explanation could be paramagnetic scattering by Hg ions. Whereas liquid mercury is diamagnetic, it is possible that valence fluctuations occur on very short time scales, as these time scales are not probed in a magnetic susceptibility experiment. If only a small fraction of the Hg-ions are not Hg++, then a significant paramagnetic signal can be expected: 1.8 barns per magnetic ($`S`$= 1/2) atom compared to 6.5 barns for the incoherent scattering cross section of an Hg ion. Of course, the magnetic form factor for a paramagnetic liquid would drop off precipitously with increasing $`q`$, and the magnetic scattering intensity would rapidly diminish with increasing $`q`$.
Supposing that the agreement between the neutron scattering line width of 2 meV ($``$3 ps-1) and the MD time scale of $``$1/3 pslivia  is more than a coincidence, we propose the following mechanism. When Hg ions aproach each other and move out again, as is the case with cage diffusion, then the overlap between the orbitals will change significantly. On close approach, the Coulomb repulsion could force one or more d-electrons into the sea of conduction electrons, whereas if a Hg ion happens to be relatively far removed from surrounding Hg ions, a non-bound electron could occupy the 4s orbital. Either of these processes would lead to a magnetic neutron scattering cross section, and these valence fluctuations would naturally occur on the same time scales as those of the cage diffusion process. This mechanism would therefore lead to an enhanced cross section for the cage diffusion process. Extrapolating the observed intensity of the 2 meV mode to $`q`$= 0, we find that the cross section for the 2 meV mode can be as high as $``$1.5 barns per Hg atom. This would correspond to the Hg ions spending significant part of their time with a valency different from 2+. This proposed scenario can easily be verified by future neutron scattering experiments: an increased temperature would lead to an even shorter time scale for the valence fluctuations, and therefore to an increased linewidth in energy.
Finally, we point out that somewhat similar effects have been observed in other liquid metals. For example, in liquid lithium it was observedli  that the measured static structure factor was larger than the calculated structure factor for $`q<`$ 1 Å-1, in liquid cesium a larger than the literature value of $`\sigma _i`$ had to be usedcs  in order to describe the neutron scattering data, and in liquid potassium it was foundk2  that the quasi-elastic coherent line width at $`q`$ = 0.4 Å-1 had decreased by only $``$ 30 % from its maximum value at $`q`$ 1.4 Å-1. We are currently investigating whether these observations can be explained by a scenario similar to the one outlined in the preceding paragraphs.
In conclusion, we have shown that cage diffusion plays a dominant role in the collective dynamics of liquid mercury. This is reflected in the wave number dependence of the characteristic width of the coherent dynamic structure factor. The expected intensity in the incoherent dynamic structure factor from the signal attributable to cage diffusion is too small to account for the observed intensity at small wave numbers, however, it is possible that valence fluctuations on the time scale of cage diffusion greatly enhance the neutron scattering cross section for this process at the smaller wave numbers.

###### Acknowledgements.

We are deeply indebted to the late Dr. Peter Verkerk, without whom this project never would have materialized. Also, we would like to thank Dr. L.E. Bove for very helpful discussions. The experiments at ISIS were supported by the Netherlands organization for Scientific Research (NWO), and part of this work is supported by grants from the Japan Society for the Promotion of Science (JSPS).

Corresponding author.


text:
# Mercury’s spin-orbit model and signature of its dynamical parameters

## 1 Introduction

Pettengill and Dyce discovered in 1965 the 3:2 spin-orbit resonance state of Mercury (the rotational and orbital periods are 56.646 and 87.969 days respectively). This particular resonance results from a functional dependance on a tidal friction adding to a non-zero eccentricity with a permanent asymmetry in the equatorial plane of the planet. Besides, in a tidally evolved system, the spin pole is expected to be trapped in a Cassini state: the orbital and rotational parameters are indeed matched in such a way that the spin pole, the orbit pole, and the solar system invariable pole remain coplanar while the spin and orbital poles on average precess at the same rate (Colombo 1965; Peale 1969, 1973).

The upcoming missions, MESSENGER (Solomon et al. 2001) and BepiColombo (Milani et al. 2001) with onboard instrumentation capable of measuring the rotational parameters of Mercury stimulate the objective to attempt to an accurate theory of the Mercury’s rotation. The current method for obtaining constraints on the state and structure of the hermean nucleus is based on the assertions introduced by Peale in 1976. According to this author, the determination of the four parameters, namely $`C_{20}`$, $`C_{22}`$, $`\eta `$ and $`\phi `$, should be sufficient to constrain size and state of the Mercury’s core ($`C_{20}`$ and $`C_{22}`$ are spherical harmonics of the second degree, $`\eta `$ is the obliquity of Mercury, and $`\phi `$ is the libration in longitude of 88 days).

In the first section of this paper, we present our approach which takes into account the coupled spin-orbit motion such as it derives from the integration of the entire solar system. In the second section, we present the dynamical behavior of the Mercury’s spin-orbit motion. In the last section, following our previous analysis of the different families of hermean rotational librations (see our previous paper Rambaux and Bois 2003), we make in evidence some accurate signatures due to parameters related to the dynamical figure of Mercury (C/MR2, $`C_{20}`$, and $`C_{22}`$) in the librations as well as in the obliquity of the planet. We find 0.45 arcseconds (as) and 2.4 milliarcseconds (mas) respectively on $`\phi `$ and $`\eta `$ for a variation of 1% on the C/MR2 value. These determinations provide new constraints on the internal structure of Mercury to be discussed with the expected accuracy forecasted in the BepiColombo mission (respectively 3.2 and 3.7 as according to Milani et al. 2001).

## 2 The Spin-Orbit N-Body Relativistic model

For obtaining the real motion of Mercury, we have used our BJV model of solar system integration including the coupled spin-orbit motion of the Moon. This model, expanded in a relativistic framework, had been previously built in accordance with the requirements of the Lunar Laser Ranging observational accuracy (Bois 2000; Bois & Vokrouhlický 1995). We extended the BJV model by generalizing the spin-orbit couplings to the terrestrial planets (Mercury, Venus, Earth, and Mars). The model is at present called SONYR (acronym of Spin-Orbit N-BodY Relativistic model). As a consequence, the SONYR model gives an accurate simultaneous integration of the spin-orbit motion of Mercury. The integration of the solar system, including the Mercury’s spin-orbit motion, uses a global reference system given by the solar system barycenter. The model is solved by modular numerical integration and controlled in function of the different physical contributions and parameters taken into account. It permits the analysis of the different families of librations and the identification of their causes (such as the planetary interactions) or their interdependences with the parameters involved in the dynamical figure of the planet. A complete description of the SONYR model is given in our previous paper (Rambaux and Bois 2003). The spin-orbit motion of Mercury is characterized by two proper frequencies (namely $`\mathrm{\Phi }`$ = 15.847 years and $`\mathrm{\Psi }`$ = 1066 years) and its 3:2 resonance presents a second synchronism, which can be understood as a spin-orbit secular resonance ($`\mathrm{\Pi }`$ = 278 898 years).

Using the SONYR model and its method of analysis (modular and controlled numerical integration, differentiation method), the librations are accurately surrounded and identified. Our previous paper Rambaux and Bois (2003) presents a detailed analysis of the whole spin-orbit mechanism of Mercury as well as its main librations.

## 3 The Mercury’s spin-orbit motion

Using our SONYR model, the present work is devoted for studiing the impact of the dynamical figure of Mercury on its spin-orbit motion. We have then analyzed the impact of the three dynamical parameters C/MR2, $`C_{20}`$, and $`C_{22}`$ on the rotational motion of Mercury.

Figure 1 presents the dynamical behavior of the Mercury’s rotation expressed by the Euler angles $`\psi ,\theta `$, and $`\phi `$ (three first panels) as well as its obliquity $`\eta `$ plotted in the fourth panel over 500 days. The Euler angles ($`\psi ,\theta ,\phi `$) related to the 3-1-3 angular sequence describe the evolution of the body-fixed axes $`Oxyz`$ with respect to the axes of the local reference frame $`OXYZ`$. Let us recall the definition used for these angles : $`\psi `$ is the precession angle of the polar axis $`Oz`$ around the reference axis $`OZ`$, $`\theta `$ is the nutation angle representing the inclination of $`Oz`$ with respect to $`OZ`$, and $`\phi `$ is the rotation around $`Oz`$ and conventionally understood as the rotation of the greatest energy (it is generally called the proper rotation). The axis of inertia around which is applied the proper rotation is called the axis of figure and defines the North pole of the rotation (Bois 1992). Let us remark that in this Figure, we have removed the mean rotation of $`58.646`$ days in the $`\phi `$ angle in order to better distinguish the librations. The resulting amplitude of libration appearing in $`\phi `$ is 20 as, which is in good agreement with Balogh and Giamperi (2002).

The integration of the SONYR model including the simultaneous spin-orbit motion of Mercury gives also the dynamical behavior of the hermean obliquity by the way of the following relation :

$$\mathrm{cos}\eta =\mathrm{cos}i\mathrm{cos}\theta +\mathrm{sin}i\mathrm{sin}\theta \mathrm{cos}(\mathrm{\Omega }\psi )$$
(1)
where $`i`$ and $`\mathrm{\Omega }`$ are respectively the inclination and ascending node of the Mercury’s orbital plane relative to the ecliptic plane (the obliquity is the angle between the polar axis and the orbital plane). In our computations presented in this paper, the used initial conditions are listed in Table 1. For $`\eta `$, $`\dot{\psi }`$, and $`\dot{\theta }`$, we use the values determined in our previous work (Rambaux and Bois 2003).

The dynamical structure of the Mercury’s rotation given in Figure 1 results from the direct solar torque whose coordinates include the planetary interactions. The effect of the gravitational torque due to Venus is five orders of magnitude smaller than the one due to the Sun. The $`P_\phi `$ rotation period of 58.646 days appears in the $`\psi `$ and $`\theta `$ angles, as well as in $`\eta `$. Signature of the $`P_\lambda `$ orbital period of 87.969 days is clearly visible in the $`\phi `$ angle (this angle is called in literature the libration in longitude of 88 days). A third period appears in the $`\psi `$, $`\theta `$, and $`\eta `$ angles, namely 176.1 days. This period $`\stackrel{~}{P}`$ results from the $`3`$:$`2`$ spin-orbit resonance ($`\stackrel{~}{P}=2P_\lambda =3P_\phi `$).

## 4 The impact of the dynamical parameters

The gravity field of Mercury is globally unknown. The tracking data from the three fly-by of Mariner 10 in 1974-75 have been re-analyzed by Anderson et al. (1987) in order to give an accurate estimation of the normalized coefficients $`C_{20}`$ and $`C_{22}`$. The resulting nominal values are C/MR2=0.34, $`C_{20}=(6.0\pm 2.0)\mathrm{\hspace{0.33em}10}^5`$, and $`C_{22}=(1.0\pm 0.5)\mathrm{\hspace{0.33em}10}^5`$. As a consequence, the dynamical coefficients $`\alpha =(CB)/A`$, $`\beta =(CA)/B`$, and $`\gamma =(BA)/C`$ (where $`A,B`$, and $`C`$ are the three principal moments of inertia of Mercury) are infered from $`C_{20}`$ and $`C_{22}`$ by the following formulae :

$`\alpha `$ $`=`$ $`{\displaystyle \frac{C_{20}2C_{22}}{C_{20}2C_{22}+C/MR^2}}`$
$`\beta `$ $`=`$ $`{\displaystyle \frac{C_{20}+2C_{22}}{C_{20}+2C_{22}+C/MR^2}}`$ (2)
$`\gamma `$ $`=`$ $`{\displaystyle \frac{4C_{22}}{C/MR^2}}`$
In these conditions, we may assume the uncertainties on the Mercury’s dynamical parameter values as dynamical variations on the Mercury’s rotation. We have computed the impact of the variations of the greatest principal moment of inertia, C/MR2, as well as $`C_{20}`$ and $`C_{22}`$, both on the obliquity and the libration in longitude. Figures 2 and 3 present the results plotted over 500 days.

These Figures show three panels corresponding to signatures of the variations of C/MR2 (a), $`C_{20}`$ (b), and $`C_{22}`$ (c) on $`\eta `$ and $`\phi `$. In each panel, 10 integrations with a step of 0.1% over the given parameter are plotted over 500 days; the maximum of amplitude is obtained for 1%. In each plot, the maximum amplitude is reached with the characteristic period of $`176.1`$ days for the obliquity $`\eta `$, and 87.969 days for the libration in longitude $`\phi `$. The results related to variations of 1% are listed in Table 2.

In the Figure 2, a variation of 1 % on the C/MR2 parameter induces a peak to peak amplitude of 2.4 mas in the obliquity, whereas with the $`C_{20}`$ and $`C_{22}`$ coefficients, the amplitudes are 1.8 and 0.7 mas respectively. The characteristic period of these signatures is $`\stackrel{~}{P}=176.1`$ days, which expresses the spin-orbit coupling of Mercury. The dynamical behavior of $`\eta `$ proves to be the right tracer for analyzing the relationship between the rotation of Mercury and its dynamical figure.

In the Figure 3, the impact of variations of the $`C_{20}`$ value on the libration in longitude is 30 times smaller than the one due to variations of C/MR2 and $`C_{22}`$. This is make clear by the fact that the dynamical parameter $`\gamma `$ depends only on $`C_{22}`$ and C/MR2. As a consequence, $`C_{20}`$ has only an indirect action on the $`\phi `$ angle. In addition, the dynamical behaviors in Figure 3, panels (a) and (c), are symmetrical; this fact reflects structure of the third equation in the formulae (2). The relation between $`C_{20}`$ and $`\phi `$ proves to be a good tracer of the $`\stackrel{~}{P}`$ spin-orbit period while the relations between C/MR2 and $`\phi `$, or $`C_{22}`$ and $`\phi `$, express the $`P_\lambda `$ period.

## 5 Conclusion

One important objective of the upcoming space missions, MESSENGER and BepiColombo, is to constrain the nature of the Mercury’s core. The theoretical framework follows the procedure proposed by Peale (1976) and Peale et al. (2002). The onboard instrumentations are planed for measuring the Mercury’s rotational parameters with a great accuracy. However, our work shows that for obtaining an excepted accuracy up to 1% on the C/MR2 coefficient, the measurements have to reach at least 0.45 as on the Mercury’s rotation and 2.4 mas on the obliquity. To our knowledge, these high accuracies are not yet achieved within the simulations.

The Mercury’s rotation strongly depends on the C/MR2 (0.45 as) and $`C_{22}`$ (0.45 as) coefficients. Besides, the obliquity of Mercury hierarchically depends on C/MR2 (2.4 mas), $`C_{20}`$ (1.8 mas), and $`C_{22}`$ (0.7 mas). The dynamical behavior of Mercury’s obliquity is finally the tracer for analyzing the non-linear relations between the rotation of Mercury and parameters of its dynamical figure.

###### Acknowledgements.

The authors thank Anne Lemaître and Sandrine D’Hoedt for useful discussions.