text:
# UT-KOMABA 01-06 Variational Approach to Dynamics of Quantum Fields 11footnote 1Lectures delivered at the NATO Advanced Study Institute on QCD Perspectives on Hot and Dense Matter, Cargèse, Aug. 6 - 18, 2001.

## 1 Introduction

In the summer of 1995, Dominique Vautherin came to Kyoto and stayed three months at the Yukawa Institute of Kyoto University where he delivered a series of very informal, pedagogical lectures on the application of variational methods to quantum field theories.  This initiated our long lasting enjoyable collaboration on the subject and Yasuhiko Tsue joined the force later. Throughout this collaboration Vautherin was always enthusiastic to explore new physics problems and always came up with new innovative ideas, much of which I suspect have origins in his expertise in nuclear many-body theory. We are very sorry that this fruitful collaboration with Vautherin came to end so soon by his untimely death and that we are no longer able to be inspired by his deep physics insights, charmed by his elegance, and most of all cheered up by his very warm presence. I like to dedicate these lectures to the memory of Dominique Vautherin from whom I learned most of the material presented below.

First, I like to say a few words about physics motivations of our works. In recent years, several authors have constructed non-trivial time-dependent solutions of classical field equations of effective meson fields.  Such solutions are relevant in considering the fate of defects which might be produced in dynamical order-disorder phase transitions in ultrarelativistic nucleus-nucleus collisions or in the evolution of very early universe. We like to study the effect of quantum fluctuations which has been ignored in these classical analyses. We may consider the classical solution as a collective mode of underlying microscopic degrees of freedom. It is well-known in nuclear many-body theory that the internal microscopic motions of nucleons are influenced by the motion of the nuclear mean field such as rotation or vibration. One may expect similar effects would arise in the field theory.

Our method is based on the Schrödinger wave functional representation of the quantum field theories in which the time-evolution of the wave function is explicitly considered. Another advantage of using the Schrödinger picture is that one can introduce the variational method to construct approximate but non-perturbative solutions of the problem as in similar problems in quantum mechanics. The extension of the variational principle for time-dependent wave function has also been developed.  The method can also be extended to deal with statistical ensembles described by Gaussian form of density matrices.

## 2 A Simple Exercise in Quantum Mechanics

We start with a very elementary example in quantum mechanics which all students learn in an introductory course of quantum mechanics: the harmonic oscillator problem. Although our interests lies in the time-dependent solution of quantum field theories, it would be instructive to first study the time-dependent solution of this simple exactly soluble problem in order to gain some physical insights into our more difficult problems since quantum field theories are nothing but an assemble of infinite number of coupled harmonic oscillators.

Let us consider a simple one-dimensional harmonic oscillator whose Hamiltonian is given by

$$H_0=\frac{1}{2}\left(p^2+\omega ^2x^2\right),$$
(1)
where we set the mass of the particle $`m=1`$ for simplicity. We look for solutions of the time-dependent Schrödinger equation:

$$i\frac{}{t}\mathrm{\Psi }(t)=H\mathrm{\Psi }(t)$$

with the quantization condition: $`[x,p]=i`$.

There are two ways to solve this problem. The easier one is the algebraic method in which one introduces ”creation” and ”annihilation” operators:

$$a^{}=\frac{1}{\sqrt{2\omega }}\left(ip+\omega x\right),a=\frac{1}{\sqrt{2\omega }}\left(ip+\omega x\right)$$

which satisfy a usual commutation relation: $`[a,a^{}]=1`$ and diagonalize the Hamiltonian: $`H=\omega \left(a^{}a+\frac{1}{2}\right)`$. The ground state of the Hamiltonian is given by the condition: $`a|0=0`$. This method is usually transcribed to the quantization of fields and one obtains particle creation and annihilation operators as basic building blocks in describing physical processes. In this approach, one usually does not refer to the wave function of the system explicitly but instead concentrates on amplitudes of particular process for given initial and final states specified by the particle number and other quantum numbers.

Alternatively, one can solve the problem by the analytic method in which one expresses the Schrödinger equation in coordinate representation with the differential operator $`p=i/x`$ and solve the resultant second order differential equation. The eigenfunction of this differential equation are Hermite’s polynomials with the Gaussian ground state wave function:

$$\mathrm{\Psi }_0(x)=x|0=\left(\frac{\omega }{\pi }\right)^{1/4}e^{\frac{1}{2}\omega x^2}$$
(2)
It is straightforward to write down the ground state wave functional of the quantized free field $`\phi (x)`$. For one component scalar field theory, it becomes just a product of the Gaussian wave functions of normal modes each specified by the momentum with the oscillator frequency $`\omega _k=\sqrt{m^2+k^2}`$:

$$\mathrm{\Psi }_0[\phi (x)]=𝒩\mathrm{exp}\left[\frac{1}{2}\underset{k}{}\omega _k\phi _k^{}\phi _k\right].$$

with a proper normalization condition.

Time-dependent variational wave function:

Now let us consider time-dependent solutions of the harmonic oscillator. We first modify the ground state wave function by adding extra complex phase factor $`e^{i\left(p_0i\omega x_0\right)x}`$. One then obtains the Gaussian wave function with its the center shifted

$$\mathrm{\Psi }(x,t;x_0,p_0)=\mathrm{exp}\left[\frac{1}{2}\omega \left(xx_0(t)\right)^2+ip_0(t)x\right]$$
(3)
Here $`p_0(t),x_0(t)`$ are time-dependent parameters which are to be determined by imposing that the above function is a solution of the time-dependent Schrödinger equation. Here we use the variational method to derive the equations of motion of $`p_0(t),x_0(t)`$.

The Schrödinger equation can be obtained by imposing a stationary condition $`\delta S=0`$ for the action:

$$S[\mathrm{\Psi }]=𝑑t\mathrm{\Psi }(t)|Hi\frac{}{t}|\mathrm{\Psi }(t)$$

with respect to variation of the wave function $`\mathrm{\Psi }(t)|`$. For the variational wave function in the form of 3, the integrand can be computed easily with $`\mathrm{\Psi }|p^2|\mathrm{\Psi }=p_0^2`$, $`\mathrm{\Psi }|x^2|\mathrm{\Psi }=x_0^2`$, $`\mathrm{\Psi }|i\frac{}{t}|\mathrm{\Psi }=\dot{p}_0x_0`$:

$$\mathrm{\Psi }|Hi\frac{}{t}|\mathrm{\Psi }=p_0^2+\omega ^2x_0^2+\dot{p}_0x_0.$$
(4)
Taking stationary conditions of the action with respect to two time-dependent parameters, $`\delta S/\delta p_0=\delta S/\delta x_0=0`$, one finds

$$\dot{p}_0=\omega ^2x_0,\dot{x}_0=p_0$$
(5)
which are just the classical equation of motion of the harmonic oscillator. The center-shifted Gaussian wave function may be obtained, beside the phase factor $`e^{ip_0x}`$, by operating a translation operator on the ground state:

$$|0;x_0,p_0\mathrm{exp}[ix_0p]|0=\mathrm{exp}\left[i\sqrt{\frac{\omega }{2}}x_0(a^{}a)\right]|0.$$

It thus describes a phase-coherent mixture of infinite number of excited states $`|n`$. The corresponding states in quantum field theory are called coherent states. One thus sees that classical field equations arises if one describes the wave functional of the quantized fields in terms of coherent states.

Keeping Gaussian form of variational wave function one can go one step further by introducing additional phase factor quadratic in $`x`$. The wave function now takes the form of

$$\mathrm{\Psi }_{\mathrm{sq}.}(x,t)=\left(\frac{\mu }{\pi }\right)^{1/4}\mathrm{exp}\left[\frac{1}{2}(\mu +i\sigma )(xx_0)^2ip_0x\right]$$
(6)
where we have introduced two new real parameters, $`\mu ,\sigma `$, characterizing the width of the wave function, which we shall consider as time-dependent variational parameters. The time-dependent parameter $`\mu (t)`$ is related to the quantum fluctuation of the position of the particle around its mean $`x_0=x`$ by $`(xx_0)^2=1/\mu (t)`$. The another parameter $`\sigma `$ is related to the rate of change of $`\mu `$ as we shall see below. This modification of the wave function thus describes the breathing motion (squeezing and stretching) of the wave function centered at $`x_0(t)`$. In the quantum field theory, these generalized coherent states are called squeezed states.

With the new variational wave function the action integrand becomes

$$\mathrm{\Psi }(t)|Hi\frac{}{t}|\mathrm{\Psi }(t)_{\mathrm{sq}.}=\frac{1}{2}\left(p_0^2+\omega ^2x_0^2\right)+\frac{1}{4}\left(\mu +\frac{\sigma ^2}{\mu }+\frac{\omega ^2}{\mu }\right)+\dot{p}_0x_0\frac{\dot{\sigma }}{4\mu }$$

Taking the stationary conditions with respect to the variations of $`\mu (t),\sigma (t)`$, in addition to classical parameters, $`x_0(t),p_0(t)`$, one finds

$$\dot{\mu }=2\sigma \mu ,\dot{\sigma }=\sigma ^2+\omega ^2\mu ^2,$$
(7)
while the equations of motion of $`x_0(t),p_0(t)`$ are unchanged. We see that the imaginary part $`\sigma `$ of the Gaussian width parameter plays a role similar to the velocity of the motion of the center of the Gaussian $`x_0`$. For the harmonic oscillator Hamiltonian, coherent states and squeezed states are exact solutions of the Schrödinger equations. The classical motion and the quantum fluctuation decouple in this exactly soluble problem. This is not the case, however, when the potential is not quadratic in $`x`$.

Anharmonicity:

To illustrate the effect of non-harmonic part of the potential, we add a term quartic in $`x`$ in our Hamiltonian:

$$H_I=\frac{\lambda }{4!}x^4$$

In the scalar field theory, this corresponds to adding a $`\phi ^4`$ self-interaction term. This anharmonic term generates new terms in the integrand of the action:

$$\mathrm{\Psi }(t)|H_I|\mathrm{\Psi }(t)_{\mathrm{sq}.}=\frac{\lambda }{4!}\left(\frac{3}{4\mu ^2}+3\frac{x_0^2}{\mu }+x_0^4\right)$$

which cause coupling between the classical motion of the mean of particle position and the quantum fluctuation around it:

$`\dot{x}_0=p_0,`$ $`\dot{p}_0=\omega ^2x_0{\displaystyle \frac{\lambda }{6}}x_0^3{\displaystyle \frac{\lambda }{4\mu }}x_0,`$
$`\dot{\mu }=2\sigma \mu ,`$ $`\dot{\sigma }=\sigma ^2+\omega ^2\mu ^2+{\displaystyle \frac{\lambda }{4\mu }}+{\displaystyle \frac{\lambda }{2}}x_0^2`$
Let us see how the time independent solution (the ground state) is modified by the anharmonic term. The conditions $`\dot{x}_0=\dot{\mu }_0=0`$ demand that $`p_0=\sigma =0`$, while the remaining two conditions $`\dot{p}_0=\dot{\sigma }_0=0`$ determine the values of $`x_0`$ and $`\mu `$. There are two types of solutions depending on the sign of $`\omega ^2`$. When $`\omega ^2>0`$ we have a ”normal” solution centered at the origin $`x_0=0`$ but the width is modified slightly as determined by:

$$\mu ^2=\omega ^2+\frac{\lambda }{4\mu }$$

This equation may be called the gap equation for the reason to be discussed in the next section. If $`\omega ^2<0`$, the potential has double minima at $`x_{\mathrm{min}.}=\pm \sqrt{6/\lambda }\omega `$ and the $`x_0=0`$ point becomes local maximum of the potential. In this case, there are two ”symmetry breaking” solutions centered at the two solutions of

$$\omega ^2+\frac{\lambda }{4\mu }+\frac{\lambda }{6}x_0^2=0$$

where the value of the width parameter $`\mu `$ is determined self-consistently with the modified gap equation:

$$\mu ^2=2\omega ^2\frac{\lambda }{2\mu }.$$

We note that the position of the center of our Gaussian variational wave function shift slightly toward the origin from the position of the minimum of the potential.

## 3 Scalar Field Theory: $`\phi ^4`$ model

Having been warmed up by a much simpler problem in quantum mechanics, it is now our task to transcribe the result to a problem in quantum field theories. We consider first a prototype scalar field theory with $`\phi ^4`$ self-interaction.

The Hamiltonian density of the theory is given by

$$(𝐱)=\frac{1}{2}\pi ^2(𝐱)+\frac{1}{2}\left(\phi (𝐱)\right)^2+\frac{m_0^2}{2}\phi ^2(𝐱)+\frac{\lambda }{24}\phi ^4(𝐱),$$
(8)
where $`\pi (𝐱)`$ is an operator conjugate to the field $`\phi (x)`$ and is expressed by a functional derivative $`i\delta /\delta \phi (𝐱)`$. 222A care must be taken to give precise mathematical meanings for these expressions, but we do not go into such problems here.  We write a Gaussian time-dependent variational wave functional for this Hamiltonian formally as

$$\mathrm{\Psi }\left[\phi (𝐱)\right]=𝒩\mathrm{exp}\left(i\overline{\pi }|\phi \overline{\phi }\phi \overline{\phi }|\frac{1}{4G}+i\mathrm{\Sigma }|\phi \overline{\phi }\right),$$
(9)
where $`G`$, $`\mathrm{\Sigma }`$, $`\overline{\phi }`$, $`\overline{\pi }`$ define respectively the real and imaginary part of the kernel of the Gaussian width and its average position and momentum. We have used the short hand notation $`\overline{\pi }|\phi =\overline{\pi }(𝐱,t)\phi (𝐱)𝑑𝐱.`$ Although it looks a little horrible, it is just a straightforward generalization of (6). The correspondences between the previous quantum mechanical example and the present case are summarized in the Table 1.

The equations of motion are found to be

$$\begin{array}{ccc}\dot{\overline{\phi }}\hfill & =\hfill & \overline{\pi },\hfill \\ & & \\ \dot{\overline{\pi }}\hfill & =\hfill & \left(\mathrm{\Delta }+m_0^2+\frac{\lambda }{6}\overline{\phi }^2(𝐱)+\frac{\lambda }{2}G(𝐱,𝐱)\right)\overline{\phi },\hfill \\ & & \\ \dot{G}\hfill & =\hfill & 2(G\mathrm{\Sigma }+\mathrm{\Sigma }G),\hfill \\ & & \\ \dot{\mathrm{\Sigma }}\hfill & =\hfill & \frac{1}{8}G^22\mathrm{\Sigma }^2\frac{1}{2}\left(\mathrm{\Delta }+m_0^2+\frac{\lambda }{2}\overline{\phi }^2+\frac{\lambda }{2}G(𝐱,𝐱)\right).\hfill \end{array}$$
(10)
where it is understood that $`\overline{\phi }`$ and $`\overline{\pi }`$ denote vectors with components $`\overline{\phi }(𝐱)`$ and $`\overline{\pi }(𝐱)`$, $`G`$ and $`\mathrm{\Sigma }`$ denote matrices with matrix elements $`G(𝐱,𝐲)`$ and $`\mathrm{\Sigma }(𝐱,𝐲)`$, and the matrix product $`G\mathrm{\Sigma }`$ is given by $`𝑑𝐳G(𝐱,𝐳)\mathrm{\Sigma }(𝐳,𝐲)`$.

For the vacuum, we should have time-independent solution so that $`\overline{\pi }=\mathrm{\Sigma }=0`$, $`\overline{\phi }(𝐱)=\phi _0`$, and

$$G(𝐱,𝐲,t)=G_0(𝐱𝐲)=𝑑𝐤\frac{e^{i(𝐱𝐲)𝐤}}{2\sqrt{\mu ^2+k^2}}.$$
(11)
where the effective mass $`\mu `$ is determined self-consistently by the non-linear integral equation

$$\mu ^2=m_0^2+\frac{\lambda }{2}G_0(0)+\frac{\lambda }{2}\phi _0^2.$$
(12)
which is usually called the gap equation because it determines the mass gap self-consistently. In the symmetric phase the expectation value of the field $`\phi _0`$ vanishes while in the symmetry broken phase it must be such that

$$m_0^2+\frac{\lambda }{2}G_0(0)+\frac{\lambda }{6}\phi _0^2=0.$$
(13)
This last equation implies that $`\mu ^2`$=$`\lambda \phi _0^2/3`$. The equations of motion (10) determine the time-evolution of the variational wave functional (9).

We comment on two well-known problems which are absent in quantum mechanics but are characteristic in quantum field theory: divergences and covariance. In the above expression the momentum integral in $`G_0(0)`$ is quadratically divergent. As well-known, this divergence originates from couplings of infinite degrees of freedom (all momentum modes) in the interaction of local fields. In the renormalizable field theories, these divergences (or cut-off dependence) may be absorbed into an appropriate redefinition of finite number of physical parameters, such as the mass and the coupling constant. This renormalization procedure works at least order by order in the power series expansion in terms of the coupling constant. In our non-perturbative calculation scheme with the variational method, we may absorb the cut-off dependence into the mass and the coupling constant. We encounter, however, well-known triviality problem: renormalized coupling constant $`\lambda _R`$ becomes zero as one sends the cut-off $`\mathrm{\Lambda }`$ to infinity keeping the original coupling $`\lambda `$ positive for the stability of the ground state. Despite this well-known pathology, the model still makes sense as an effective theory with the finite cut-off.

The another problem is the lack of manifest covariance in our formulation: the time coordinate $`t`$ have been treated differently from spatial coordinates throughout calculations. This is an old problem which was originally solved in QED by Tomonaga, Schwinger, Feynman and Dyson who developed the covariant perturbation theory based on the interaction representation. Vautherin invented a new ingenious trick to derive manifest covariant form of the equations of motion which I shall now describe.

Mean field equations in the Hartree-Bogoliubov form:

A key ingredient of his method is the reduced density matrix defined by

$$(𝐱,𝐲;t)=\left(\begin{array}{cc}i\widehat{\phi }(𝐱)\widehat{\pi }(𝐲)1/2& \widehat{\phi }(𝐱)\widehat{\phi }(𝐲)\\ \widehat{\pi }(𝐱)\widehat{\pi }(𝐲)& i\widehat{\pi }(𝐱)\widehat{\phi }(𝐲)1/2\end{array}\right),$$
(14)
where $`\widehat{\phi }=\phi \overline{\phi }`$ , $`\widehat{\pi }=\pi \overline{\pi }`$, and expectation values are calculated with the Gaussian functional $`\mathrm{\Psi }(t)`$ and is given by

$$=\left(\begin{array}{cc}2iG\mathrm{\Sigma }& G\\ \frac{1}{4}G^1+4\mathrm{\Sigma }G\mathrm{\Sigma }& 2i\mathrm{\Sigma }G\end{array}\right).$$
(15)
Using (10), one can show that the equation of motion of the reduced density matrix can be cast into the Liouville-von Neumann form

$$i\dot{}=[,],$$
(16)
where the generalized Hamiltonian $``$ is given by

$$=\left(\begin{array}{cc}0\hfill & 1\hfill \\ \mathrm{\Gamma }\hfill & 0\hfill \end{array}\right).$$
(17)
with

$$\mathrm{\Gamma }=\mathrm{\Delta }+m_0^2+\frac{\lambda }{2}\overline{\phi }^2+\frac{\lambda }{2}G(𝐱,𝐱).$$
(18)
This form of equations is known in many-body theory as the time-dependent Hartree-Bogoliubov equations.

The reduced density matrix satisfies $`^2=\frac{1}{4}𝐈`$ so that it has two eigenvalues of $`\pm 1/2`$. We write the $`n`$-th eigenvector of $``$ with eigenvalue $`1/2`$ as $`(u_n(𝐱),v_n(𝐱))`$:

$$\left(\begin{array}{c}u_n\hfill \\ v_n\hfill \end{array}\right)=\frac{1}{2}\left(\begin{array}{c}u_n\hfill \\ v_n\hfill \end{array}\right),$$
(19)
then one can show that vectors $`(u_n^{}(𝐱),v_n^{}(𝐱))`$ give eigenvectors for eigenvalue $`1/2`$. The $`u`$ and $`v`$ components of eigenvectors are called mode functions. With these eigenvectors the reduced density matrix has a spectral decomposition as

$$=\frac{1}{2}\underset{n>0}{}\left[\left(\begin{array}{c}u_n\hfill \\ v_n\hfill \end{array}\right)(v_n^{},u_n^{})+\left(\begin{array}{c}u_n^{}\hfill \\ v_n^{}\hfill \end{array}\right)(v_n,u_n)\right].$$
(20)
The Liouville-von Neumann equation (16) can be rewritten in terms of the mode functions as

$$i_t\left(\begin{array}{c}u_n\hfill \\ v_n\hfill \end{array}\right)=\left(\begin{array}{cc}0\hfill & 1\hfill \\ \mathrm{\Gamma }\hfill & 0\hfill \end{array}\right)\left(\begin{array}{c}u_n\hfill \\ v_n\hfill \end{array}\right).$$
(21)
Eliminating $`v_n`$ and inserting (18) we obtain a modified Klein-Gordon type equation for the mode functions

$$\left(\mathrm{}+m_0^2+\frac{\lambda }{2}\overline{\phi }^2+\frac{\lambda }{2}G(𝐱,𝐱)\right)u_n=0,$$
(22)
where the spectral representation (20) implies

$$G(𝐱,𝐱)=𝐱|G(t)|𝐱=\frac{1}{2}\underset{n}{}|u_n(𝐱,t)|^2.$$

To write the equation of motion fully covariant way, we introduce a Feynman propagator in terms of the mode functions

$$\begin{array}{ccc}x|S|y\hfill & =\hfill & \theta (x_0y_0)_{n>0}u_n^{}(𝐱,x_0)u_n(𝐲,y_0)\hfill \\ & & \\ & & +\theta (y_0x_0)_{n<0}u_n^{}(𝐱,x_0)u_n(𝐲,y_0).\hfill \end{array}$$
(23)
so that

$$𝐱|G(x_0)|𝐱=x|S|x$$
(24)
Then we finally arrive at a fully covariant, self-consistent equations of motion for the mode functions which reproduces the same equations obtained earlier by the functional integral method.

This way of writing the equations of motion also paves a way to generalize the calculation at finite temperatures.

## 4 Statistical Ensembles

Foregoing discussions are limited to evolution of a single coherent (Gaussian) state. In realistic physical situations, we are more interested in the evolution of the statistical ensemble which is described by the density matrix:

$$\widehat{\rho }(t)=\underset{n}{}|\mathrm{\Psi }_n(t)p_n(t)\mathrm{\Psi }_n(t)|$$
(25)
where $`p_n(t)`$ is a probability distribution specifying the ensemble hence It should satisfy $`_np_n(t)=1`$. The expectation value of an observable $`O`$ is given by $`O=\mathrm{Tr}\widehat{\rho }O`$ and the statistical entropy $`S`$ is given by

$$S=\mathrm{Tr}\widehat{\rho }\mathrm{ln}\widehat{\rho }=\underset{n}{}p_n\mathrm{ln}p_n$$
(26)
Hence if $`p_n`$ is time-independent then the entropy is conserved. In equilibrium at temperature $`T=1/\beta `$, the density matrix is given by the canonical ensemble: $`p_n^{eq.}(t)=e^{\beta E_n}/Z`$ where $`Z=_ne^{\beta E_n}`$ which maximizes the entropy $`S(t)`$ under the condition of fixed expectation value of the energy $`E=H`$. The density matrix with time-independent $`p_n`$ obeys the Liouville-von Neumann equation: $`i_t\widehat{\rho }=[H,\widehat{\rho }]`$.

The variational method has been extended for the time-dependent density matrix by Eboli, Jackiw and Pi.  Without going into detail, we illustrate the essence of their method in terms of simple quantum mechanical example with harmonic oscillator Hamiltonian (1). We first introduce the coordinate representation of the density matrix by

$$\rho (x,y;t)=x|\widehat{\rho }(t)|y=\underset{n}{}\mathrm{\Psi }_n(x,t)p_n(t)\mathrm{\Psi }_n^{}(y,t)$$
(27)
and observe that with single Gaussian variational wave function the density matrix is just a product of two Gaussian: $`\rho (x,y;t)e^{\omega (x^2+y^2)/2}`$. One can show that in the other extreme limit of thermal equilibrium, the density matrix $`\rho (x,y;t)`$ can also be expressed by a mixed Gaussian form. This is so because the equilibrium density matrix $`\widehat{\rho }_{eq.}(\beta )=e^{\beta H}`$ obeys the imaginary time Schrödinger equation (the Bloch equation): $`_\tau \widehat{\rho }_{eq.}(\tau )=H\widehat{\rho }_{eq.}`$. Using $`\rho _{eq.}(x,y)=x|\widehat{U}(i\tau )|y`$ and well-known path-integral expression of the matrix elements of the unitary evolution operator $`\widehat{U}(t)=e^{itH}`$, one finds

$$\rho _{eq.}(x,y)=\left(\frac{\omega }{2\pi \mathrm{sinh}\omega \beta }\right)^{1/2}\mathrm{exp}\left[\frac{\omega }{2\mathrm{sinh}\omega \beta }\left\{(x^2+y^2)\mathrm{cosh}\omega \beta 2xy\right\}\right]$$
(28)
which is again Gaussian with an extra term containing a product of two coordinates $`xy`$.

For more general mixed states we introduce a generalized Gaussian density matrix

$$\rho (x,y;t)=𝒩\mathrm{exp}\left[ip_0(\widehat{x}\widehat{y})\frac{\mu }{2}\left(\widehat{x}^2+\widehat{y}^22\xi \widehat{x}\widehat{y}\right)+i\frac{\sigma }{2}(\widehat{x}^2\widehat{y}^2)\right]$$
(29)
where we have used the time-dependent shifted coordinates: $`\widehat{x}=xx_0`$ and $`\widehat{y}=yx_0`$. Equations of motion of the time-dependent parameters of the generalized Gaussian density matrix (29) can be derived from the Liouville equation of the density matrix and we obtain a set of equation similar to (30) with a small modification:

$$\dot{\mu }=2\sigma \mu ,\dot{\sigma }=\sigma ^2+\omega ^2(1\xi ^2)\mu ^2.$$
(30)
The parameter $`\xi `$ is called the mixing parameter which measures the degree of mixture of different pure states in the ensemble; it remains constant for an adiabatic evolution of the system. In equilibrium, $`p_0=x_0=\sigma =0`$ and other two parameters are given by the specific functions of temperature:

$$\mu _{eq.}=\omega \mathrm{coth}\omega \beta ,\xi _{eq.}=\mathrm{cosh}^1\omega \beta $$
(31)
as indicated by the formula (28). Extension of the Gaussian density matrix in quantum mechanics to that in quantum field theories is straightforward as indicated in the last row of the Table 1.

The reduced density matrices we have introduced in the previous section for a pure Gaussian state can be extended for a mixed state immediately by the replacements:

$$\widehat{\phi }(𝐱)\widehat{\pi }(𝐲)=\mathrm{Tr}(\widehat{\rho }(t)\widehat{\phi }(𝐱)\widehat{\pi }(𝐲)),\widehat{\phi }(𝐱)\widehat{\phi }(𝐲)=\mathrm{Tr}(\widehat{\rho }(t)\widehat{\phi }(𝐱)\widehat{\phi }(𝐲)),\mathrm{etc}.$$
(32)
and one can derive equations of motion of the reduced density matrix similar to the pure state case. In the case of equilibrium distribution, this amounts to introduce a factor containing the occupation number in the sum over the mode functions.

## 5 Rotating chiral condensate in the $`O(N)`$ sigma model

We briefly mention about an application of the above method to the $`O(N)`$ sigma model which is composed of $`N`$-components coupled scalar fields $`\phi _n`$ with continuous $`O(N)`$ symmetry. We expect that this global symmetry of the model is broken spontaneously at low temperatures, characterized by non vanishing expectation value of one component of fields, say $`\overline{\phi }_0=\mathrm{Tr}\widehat{\rho }\widehat{\phi }0`$, and the system exhibits an order-disorder phase transition to a state with $`\overline{\phi }_0=0`$. This is what is expected in QCD where the chiral symmetry is broken in vacuum and is expected to be restored at finite temperature. This chiral phase transition has been studied by an effective theory with pion and sigma fields with $`O(4)`$ global symmetry.

We have applied our method to describe a special kind of time-dependent condensate which rotates in a subspace of the internal symmetry space:

$$\left(\begin{array}{c}\overline{\phi }_1(x)\\ \overline{\phi }_2(x)\end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}(qx)& \mathrm{sin}(qx)\\ \mathrm{sin}(qx)& \mathrm{cos}(qx)\end{array}\right)\left(\begin{array}{c}\phi _0\\ 0\end{array}\right)=\mathrm{exp}[i(qx)\tau _2]\left(\begin{array}{c}\phi _0\\ 0\end{array}\right)$$
(33)
where $`\tau _2=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right)`$ and $`q_\mu =(\omega ,𝐪)`$ is a four vector specifying the direction of rotation: for a spatially uniform rotation in time it is a pure time-like vector; while a static condensate with oscillation in space, it is a pure space-like vector. The matrix $`U(x)=\mathrm{exp}[i(qx)\tau _2]`$ can be considered as a ”gauge transformation” to the local rest frame of the rotating condensate. In the rotating frame, the four derivative which appear in the equation of motion of the mode functions is ”gauge transformed” to $`U_\mu U^{}=_\mu iq_\mu \tau _2`$ and the effect of the rotation is seen in this frame as an appearance of the apparent ”centrifugal force”. Indeed, our result of phase diagram for such dynamical condensates shows that the amplitude of the chiral condensate with uniform time-like rotation increases due to the centrifugal force. This effect has been observed in the classical solutions of Anselm and Ryskin ; our quantum generalization of their solution shows that this effect is amplified by the coupling of quantum fluctuations to rotations, which the static condensate with spatial oscillations are more suppressed by the quantum fluctuations. A phase diagram of rotating condensate was obtained in  and the damping of the rotating condensate due to the symmetry breaking perturbation was computed by the method of the response function in .

## 6 Outlook

Vautherin started to work on the variational approach to quantum field theories many years ago with Arthur Kerman. They developed many important ideas in their unpublished works and tried to solve QCD non-perturbatively with their method with a hope to gain new insights in the quark confinement problem. The Gaussian Ansatz for the variational wave functional however has difficulty of breaking the local gauge invariance and the projection to color singlet state destroys a nice feature of the Gaussian wave functional. Vautherin continued to work on the problem with his students and brought a new insight into the problem again introducing a technique developed in nuclear many-body theory in his last paper. His efforts in this direction may be carried over to study the dynamical evolution of the quark-gluon plasma in ultrarelativistic nucleus-nucleus collisions. Vautherin was also interested in the recent experimental breakthrough of creating weakly interacting Bose-Einstein condensates in well-controlled laboratory environments. Our method can be also applied to such problem to investigate the effect of quantum fluctuations which are usually ignored in theoretical descriptions.

I am much indebted to Yasuhiko Tsue as well as to Dominique Vautherin for our works quoted above. I thank them for sharing the joy of physics in our collaboration.


format: text

text:
# Discerning the neutron density distribution of 208Pb from nucleon elastic scattering

## I Introduction

Interest in the matter distributions of 208Pb, and its neutron density profile particularly, is quite topical . There is a proposal to measure its neutron root-mean-square (rms) radius at the Jefferson Laboratory  from an analysis of parity-violating electron scattering data. In contrast to proton rms radii that are known to within an accuracy $`0.02`$ fm , neutron rms radii are less certain.

Recently, the neutron rms radius in 208Pb was assessed in terms of modern Skyrme-Hartree-Fock (SHF) models . With the Friedman-Pandharipande (FP) neutron equation of state  as a constraint, the neutron rms radius in 208Pb was expected to be $`0.16\pm 0.02`$ fm larger than the proton value. Previous estimates of this neutron skin, $`S=\sqrt{r_n^2}\sqrt{r_p^2}`$, ranged from 0.1 fm to 0.3 fm ; the lower values favored in general by SHF models, while relativistic mean field models predict values closer to 0.3 fm . Knowledge of the skin thickness in 208Pb then is a good constraint upon such models of structure .

The planned parity-violating electron scattering experiment  will only provide information about the neutron rms radius itself. We seek further information and address the question of whether analyses of nucleon scattering data establish a measure of the neutron density distribution in 208Pb. Hadron scattering data have been analyzed previously to deduce the neutron skin thickness in 208Pb. The ratio of $`\pi ^+`$ and $`\pi ^{}`$ reaction cross sections gave $`S=0.0\pm 0.1`$ fm . Analyses of elastic proton scattering data at 0.8 GeV gave $`0.14\pm 0.04`$ fm . However, a review of the analysis of proton scattering from 40Ca  gave values in the range of $``$0.4 fm to $``$0.2 fm for $`S`$ in that nucleus which are systematically smaller than all theoretical models which give $``$0.05 fm , and suggests that there is a systematic problem in the analysis using phenomenological models of high energy proton scattering data which would affect the extracted $`S`$ values at the level of $`0.2`$ fm. A more recent analysis of 650 MeV proton scattering data  gave $`S=0.20\pm 0.04`$ fm for 208Pb while that reaction model gave a result for 40Ca which is consistent with theoretical predictions. The excitation of the isovector giant dipole resonance in 208Pb by inelastic alpha scattering  was used to deduce $`S=0.19\pm 0.09`$ fm. These previous analyses, based on phenomenological models, produce a range of results for 208Pb and suggest that there may be small but systematic errors in those strong interaction models for hadron scattering limiting the accuracy in the extraction of $`S`$ to $`0.2`$ fm.

Our approach is based on coordinate space nonlocal optical potentials generated by a full folding of realistic effective nucleon-nucleon ($`NN`$) interactions with ground state density matrices (termed densities hereafter) of 208Pb. This allows us to distinguish between various model structures, including those of the SHF type proposed by Brown , even if two may have the same rms radii. As a calibration of the use of SHF models we consider the elastic scattering from 40Ca as well. Further, as the effective $`NN`$ interaction is dominated by the isoscalar $`{}_{}{}^{3}S_{1}^{}`$ channel , proton scattering predominantly will probe the neutron density and vice-versa and so we consider both proton and neutron elastic scattering at a given energy seeking as complete a map as possible of the nuclear matter distributions.

## II Nuclear models and the microscopic optical potential

The model with which predictions of nucleon-nucleus ($`NA`$) scattering observables are made has been given in detail in a recent review . Use of the complex, nonlocal, $`NA`$ optical potentials defined by that model prescription, without localization of the exchange amplitudes, gave predictions of differential cross sections and spin observables that are in good agreement with data from many nuclei (3He to 238U) and for a wide range of energies (40 to 300 MeV). Crucial to that success was the use of effective $`NN`$ interactions built upon $`NN`$ $`g`$ matrices; solutions of Brueckner-Bethe-Goldstone equations for realistic starting (free) $`NN`$ interactions. The $`NA`$ optical potentials result from folding those effective interactions with the densities of the target nucleus. That folding includes the antisymmetrization of the projectile-nucleus wave functions and therefore exchange (knock-out) amplitudes are treated explicitly. Consequently the $`NA`$ potentials inherently are non-local. The optical potentials that result from that process we term for brevity as $`g`$-folding potentials.

Recently, this approach was applied successfully to make predictions of the integral observables of nucleon elastic scattering . Thus, the use of the $`g`$-folding optical potentials give good predictions to both angular-dependent and integral observables; a result not guaranteed with the more common phenomenological approaches. Of import, however, is that the level of agreement with data in the $`g`$-folding model depends on the quality of the underlying model of structure. We seek to use that dependence as a sensitive evaluation of the densities considered.

Our theoretical density distributions are based upon the Skyrme Hartree-Fock model for 208Pb with a spherical closed-shell configuration. It was shown in  that there is a combination of parameters in the Skyrme Hamiltonian (dominated by the $`x_3`$ parameter) which has a strong influence on the neutron skin thickness but which are not well determined by the data on binding energies and charge radii. This combination can be correlated with the pressure in the neutron equation of state at normal nuclear density . It can also be related to the surface symmetry energy . The Skyrme interaction SKX  was obtained with a constraint on the neutron equation of state provided by the FP model  which constrains the neutron skin in 208Pb to be $`S=0.16\pm 0.02`$ fm. The results obtained with the SKX model will be denoted herein by SHF1. If one does not allow for any constraint from a model for the neutron equation of state then a much larger range of $`S`$ is allowed. In particular, we use a model which is constrained by the same nuclear properties used for SKX but with a value of the parameters which gives $`S=0.25`$ fm for 208Pb. This we denote as SHF2. The densities for 208Pb obtained from these models, as well as the others considered herein, are shown in Fig. 1.

An initial test of these interactions is provided by the elastic electron scattering data which yields information on the charge density which, in turn, gives information on the proton density. SKX appears to give an excellent reproduction of the charge-density distribution (Fig. 11 of ). However, there is some model dependence in the extraction of the charge density from electron scattering data. It is better to compare to a representation which is more closely associated with the actual data – this is the plane-wave transform of the charge density shown in Fig. 2. The experimental form factor is obtained from the charge density distribution given in . Experiment is compared with the SHF1 and SHF2 model results showing a disagreement with data which systematically increases with momentum transfer. The main feature of the distribution which affects the high-$`q`$ behavior is the surface thickness. As discussed in , the SKX interaction appears to have a surface thickness which is a little sharper than that determined from experiment. We have thus looked at other Skyrme interactions in terms of the data in Fig. 2 and find the older SKM\* interaction  is much better than others in this regard. SKM\* appears to achieve this improvement by a decrease in the power of the density, $`\rho ^\alpha `$, associated with the density-dependent part of the interaction from its value $`\alpha =1/2`$ for SKX to $`\alpha =1/6`$ for SKM\*. Coincidentally, SKM\* predicts a neutron skin of 0.17 fm which is essentially the same as that predicted by SKX (SHF1 model). This is obtained mainly because the $`x_3`$ parameter was set to zero by default, since it is not well determined by nuclear data.

The value of $`\alpha `$ is also associated with the incompressibility coefficient $`K`$ for infinite nuclear matter which ranges from $`K=270`$ MeV for SKX to $`K=217`$ MeV for SKM\*. It was found in  that $`\alpha =1/2`$ gave the best overall fit with the data set considered and that when $`\alpha `$ is decreased (as in the SKXm interaction, which has $`\alpha =1/3`$) the overall $`\chi ^2`$ increased mainly because of an increase in the contribution from the single-particle energies. It was also argued in  that an improvement in the surface properties may require an additional parameter in the Skyrme Hamiltonian associated with the next order $`d`$-wave term in the expansion in terms of the range of the $`NN`$ interaction. Thus the present models are not perfect but they are good enough for a discussion of the effects of the neutron skin and the surface thickness on the proton and neutron scattering data. We compare results obtained with three Skyrme interactions SHF1 (SKX), SHF2, and SKM\* which will enable us to explore the effect of neutron skin (a comparison of SHF1 and SHF2) with a fixed surface thickness, and surface thickness (a comparison of SHF1 and SKM\*) with a fixed neutron skin.

We compare also with results obtained with the simple harmonic-oscillator radial wave functions for 208Pb which were used in . Two sets of oscillator parameters were used. For HO1 we use $`\mathrm{}\omega =6.70`$ MeV for both protons and neutrons which is chosen to give the rms charge radius of 5.50 fm (a proton rms radius 5.45 fm). HO1 has a neutron rms radius of 5.84 fm. For HO2 the oscillator parameter for neutrons was changed to 7.25 MeV to decrease the rms neutron radius to 5.61 fm so that the neutron skin $`S=0.16`$ fm is close to that obtained with SHF1 and SKM\*.

The rms radii from all models of the ground state of 208Pb considered are listed in Table I. All five give essentially the same radius for the protons but they vary considerably in the radius for the neutrons. The difference between the neutron and proton rms radii given in the last column emphasizes that spread. Note that the HO2 model was chosen to give the same rms radii as those of the SHF1 model. The neutron radius obtained from the SKM\* model is similar to those of the SHF1 and HO2 models. However, as is evident in Fig. 1, each model gives distinctive density distributions. The normalization used is such that their volume integrals equate to the proton and neutron numbers of 82 and 126 respectively. In Fig. 1, the proton density $`\rho _p(r)`$ of both the SHF1 and SHF2 models are displayed by the solid curve. The density obtained from the SKM\* model is given by the double-dot-dashed line and exhibits a larger diffuseness compared to the densities of the other Skyrme models. Likewise, both the HO1 and HO2 models have proton densities as given by the dot-dashed curve. These quite distinct shapes nevertheless give the same proton rms radius. However, they differ in the longitudinal electron scattering form factor as shown in Fig. 2.

The five model neutron densities $`\rho _n(r)`$ are also shown in Fig. 1. The SHF1 and SHF2 neutron densities are displayed by the solid and dashed lines respectively. They are similar with the SHF2 density having a slightly more diffuse surface region; a property that results in the larger neutron rms radius. The density from the SKM\* model is given by the double-dot-dashed line and, as for the proton density, exhibits a larger diffuseness compared to the densities of the other Skyrme models. The neutron densities of the HO1 and HO2 models are displayed by the dotted and dot-dashed lines respectively. As with their proton densities, the neutron densities of both of these models are enhanced in the nuclear interior over the SHF values. But these HO densities also have increased neutron probability at the surface. Recall that the HO2 model was set to have the same neutron rms radius as the SHF1 prescription.

The SHF (SKX) densities for 40Ca are displayed in Fig. 3. Also displayed therein are the densities from the oscillator model used by Karataglidis and Chadwick , for which $`\mathrm{}\omega =10.25`$ MeV. In the surface region, the two models predict essentially the same densities. As those densities differ markedly only well within the nuclear volume, we expect influence on scattering primarily at high momentum transfer scattering results at energies for which absorption through the nucleus is not too large. We anticipate differences in cross sections for scattering at energies $`200`$ MeV and at momentum transfer values $`1`$ fm-1.

## III Results

We have analyzed both proton and (where available) neutron elastic scattering data from 40Ca and 208Pb at 40, 65, and 200 MeV. The choices of energies were predicated in part on the availability of data and of the momentum transfer values at which those data have been measured. In addition our choice was influenced by previous applications of the $`g`$-folding potentials with those energies being quite successful . Furthermore, the effective interactions defined for each energy are quite different due to the energy dependence of medium effects in the basic $`g`$ matrices so that the set of $`NA`$ scattering results we obtain provide a consistency check on the various models of structure used.

The differential cross sections for 40 MeV proton and neutron elastic scattering from 40Ca are presented in Fig. 4. Therein, the proton scattering data of Camis et al.  and the neutron scattering data of de Vito et al.  are compared to the results of the calculations made using SHF (solid line) and HO (dashed line) models. The SHF and HO model results equally well describe the data, although they underpredict the proton scattering in the regions of the minima. That disagreement may be due to problems in specifying the effective interaction at low energies . However, further comment should await the consideration of the comparisons of 40 MeV proton scattering from 208Pb. Whatever any deficiency at this energy there may be though does not affect results we find for, and conclusions we may draw from, scattering at 65 and 200 MeV.

The 65 MeV elastic scattering scattering cross sections for 40Ca are displayed in Fig. 5. Therein, the agreement between the model results and the proton scattering data of Sakaguchi et al.  is now very good; much better than that found in Fig. 4. Also, there is a slight preference for the SHF result at larger angles. In the case of the neutron scattering results, both model results agree but underpredict some of the data of Hjort et al. . However, those data are somewhat problematic as concluded from a recent analysis  in which several sets of data were compared with theory at that energy. A new measurement of this cross section is required to resolve any such problem.

It is with 200 MeV scattering from 40Ca, displayed in Fig. 6, that we observe significant differences between the predictions of scattering made using the SHF and HO models. For proton scattering that is displayed in Fig. 6(a), the SHF model result agrees well with the data of Hutcheon et al.  (circles) and Seifert et al.  (squares) and is a significant improvement on the result found using the HO model. This variation concurs with our expectation that the differences between the inner radial densities of the two models of structure would influence the results of scattering calculations at higher energies. Taking the results at all three energies, we believe that the SHF model is the better description of 40Ca.

We now turn to 208Pb and consider first the integral observables from nucleon scattering as a test of sensitivity to the matter distributions of 208Pb. Those integral observables at 40, 65, and 200 MeV are given in Tables II and III and we note that a study of these quantities with two of these structure models has been made recently  for energies 10 to 300 MeV. The total reaction cross sections for both proton and neutron scattering from 208Pb are listed in Table II. Comparing the results of the calculations at both 40 and 65 MeV with the available proton data indicates a preference for both SHF, the SKM\*, and HO2 models of the ground state density. However, comparing the model results of the neutron total reaction cross section with the available evaluated data at those energies  gives a preference for the Skyrme models only. The predicted proton and neutron total reaction cross sections at 200 MeV vary sufficiently that their measurements would be desirable.

The results of our calculations for the total neutron cross sections are given for the three energies in Table III and are compared with the data of Finlay et al. . For 40 MeV the Skyrme models are preferred, although all results overestimate the measured value. At 65 MeV the Skyrme model results agree well with the data, the SKM\* model result doing best of all, and do better than both HO predictions. This is not the case at 200 MeV, where all results predict the measured value reasonably well. While these total reaction and total cross section results together indicate a preference for the Skyrme models, we need additional evidence. We consider then the angular distributions of each scattering.

The differential cross sections for the scattering of 40, 65 and 200 MeV nucleons from 208Pb are presented in Figs. 7, 8, and 9 respectively. At 40 MeV, proton elastic scattering is shown in Fig. 7 and evidently there is little if any differentiation between the SHF1, SHF2, and SKM\* model results. They both compare well with the data of Blumberg et al. . Note, however, that both SHF and the SKM\* calculated results agree much better with the data than do those found using the HO models of structure. The quality of reproduction of the data in this case is in stark contrast to what we found at 40 MeV with 40Ca. If that is to remain a problem with specification of the effective interaction then it seems to be a nucleus-dependent effect. In the case of neutron scattering, the results of all model calculations agree quite well with the data of de Vito et al. . The results for 208Pb are more distinctive than those for 40Ca at this and other energies, but so are the distinctions between the model proton densities for both nuclei. However, these results indicate that while the nucleon densities in 208Pb are better described by the SHF and SKM\* models nucleon scattering at this energy is largely sensitive to the surface properties only. Only in the surface are the proton densities still sufficiently similar for neutron scattering results of all models to be as alike as they are.

Recall that the integral observables given in Table II make preference to the SHF and SKM\* models of the density. Similar preference is indicated by the differential cross sections for 65 MeV scattering that are displayed in Fig. 8. Considering the neutron scattering first, all models give similar results in good agreement with the data of Ibaraki et al. , although above $`40^{}`$ the SHF and SKM\* model results clearly do better. For proton scattering, the SHF and SKM\* models are both in good agreement with each other and with the data and differ only slightly from those of the HO models. As for 40 MeV scattering, the integral observables at 65 MeV concur with these findings that favor the Skyrme models of the density.

The differential cross sections for 200 MeV nucleon scattering from 208Pb are shown in Fig. 9. In this case only 200 MeV proton scattering data exist  for comparison with our predictions. Nevertheless, that comparison confirms the findings from analyses of the lower energy data, namely that the SHF and SKM\* models are the better descriptions of the densities of 208Pb. In addition, however, some discrimination between the three Skyrme model results is evident at this energy. Of the three Skyrme models, the SKM\* model result agrees best of all with the data while the SHF1 does worst. Given that the rms radii from these two models are very similar, the marked difference indicates a sensitivity to the diffuseness in the density. At this energy, unlike those at lower energies, the neutron scattering results show marked differences between the SHF and HO calculations. Again, we expect that this is due to the scattering of the higher energy probe being more influenced by the bulk nuclear medium properties of the densities. Consequently, a measurement of the 200 MeV neutron elastic scattering, angular distribution and associated integral observables, is certainly desirable.

## IV Conclusions

The distinctions between the predictions of nucleon scattering from 208Pb at these three energies, and found with the set of five model structures used, suffice to select that which most likely prescribes the actual matter densities of the nucleus. A similar conclusion is reached for 40Ca for the models considered herein. Specifically we contend that the use of $`g`$-folding potential model calculations can differentiate between different model structures of the neutron density so as to pin down the neutron rms radius far better than has been possible in the past. As well the process gives a good appraisal of the actual density distribution. For 208Pb in particular, our analyses indicate that the SKM\* model gives the best representation of the density. Together with analyses of the longitudinal elastic electron scattering form factor it suggests a neutron skin thickness for 208Pb of 0.17 fm; a value consistent with expectations of the SHF1 model nucleus, which is constrained by the FP neutron equation of state. The only difference between the two models is that of a larger diffuseness for the SKM\* model, accounting for both the agreement between the results and data for both nucleon and electron scattering. This would suggest a need to extend the SKX models to predict a larger diffuseness, for example, by the addition of a $`d`$-wave term in the Skyrme Hamiltonian.

One should also note that while a measurement of the skin thickness as proposed in the experiment for the Jefferson Laboratory  is important, that quantity is a volume property of the nucleon distributions. Other information is required to specify a more complete picture of the neutron density. From our studies it seems that simultaneous analyses of angular and integral observables are relevant. Given that the HO2, SHF1, and SKM\* models of structure we have used predict essentially the same skin thickness for 208Pb but give significantly different predictions when used to generate $`NA`$ $`g`$ folding optical potentials for the nucleon scattering, as well as electron scattering form factors, analyses of complementary nucleon and electron scattering data permit one to discern such finer details of densities.

This work was supported in part by a grant from the Australian Research Council, U.S. National Science Foundation grants no. PHY-9605207 and PHY-0070911, and DOE Contract no. W-7405-ENG-36. Two of us (S.K. and K.A.) would like to thank the National Superconduncting Cyclotron Laboratory of Michigan State University for their kind hospitality during part of this work, while another one of us (K.A.) would like to thank the Theoretical Division of the Los Alamos National Laboratory for its similar kind hospitality.


text:
# Varying speed of light cosmology from a stringy short distance cutoff

## I  Introduction

Varying speed of light (VSL) cosmology was conceived as an alternative solution to the cosmological horizon problem Moffat (1993a, b); Albrecht and Magueijo (1999); Magueijo (2000). In contrast with inflation, where the comoving horizon scale shrinks as a consequence of the accelerating growth of the scale factor, VSL keeps the scale factor dynamics of the standard hot big bang unchanged but assumes that the speed of light, $`c(t)`$, is a decreasing function of of cosmic time $`t`$. The horizon problem can be solved provided that $`c(t)`$ drops sufficiently steeply for some time. In the original formulation, $`c(t)`$ was treated as a free parameter of the theory.

In Ref. Alexander and Magueijo (2001), it was noted that VSL can effectively be obtained from the thermodynamical properties of theories with nonlinear dispersion relations which, in turn, were motivated in Alexander and Magueijo (2001) by noncommutative geometry. While there are indications from string/M-theory that noncommutativity of spatial coordinate operators may be relevant near the string scale, we are lacking a simple, concrete realization of noncommutative geometry that can be used in the cosmological context. Hence one is forced to reduce its effects to a modified dispersion relation.

A very similar yet more specific approach is taken in this work. The starting point is an observation that was made in numerous studies of quantum gravity and string theory: under rather general assumptions, one can show that there exists a minimal length scale that can be probed by experiments whereas the resolution of shorter distances is prohibited by quantum gravitational effects (e.g., Gross and Mende (1988); Amati et al. (1989); Garay (1995); Witten (1996)). Examined from the point of view of the low energy effective theory, it appears that quantum gravity modifies the Heisenberg uncertainty relation to

$$\mathrm{\Delta }x\mathrm{\Delta }p\frac{\mathrm{}}{2}\left(1+\beta (\mathrm{\Delta }p)^2+\mathrm{}\right)$$
(1)
which gives rise to a minimal short distance uncertainty $`\mathrm{\Delta }x_{\mathrm{min}}=\mathrm{}\sqrt{\beta }`$ even for infinite $`\mathrm{\Delta }p`$. For practical purposes, $`\mathrm{\Delta }x_{\mathrm{min}}`$ can be identified with either the string scale or the Planck length, where the former could be a few orders of magnitudes larger than the latter. The correction may of course originate from complicated dynamics of the underlying fundamental theory. On the other hand, it can be modeled by equivalent corrections to the canonical commutation relation:

$$[𝐱,𝐩]=i\mathrm{}(1+\beta 𝐩^2+\mathrm{})$$
(2)
as discussed in Kempf (1994). Indeed, the uncertainty relation above belongs to class of only very few types of short-distance structures of space-time that are admitted under very general assumptions Kempf (2000). Hilbert space representations of Eq. 2 were employed, for instance, for regularizing field theory Kempf and Mangano (1997) and, more recently, for analyzing the impact of short distance uncertainty on the predictions of inflation Kempf (2001); Kempf and Niemeyer (2001).

The finite spatial localization of particles also changes their thermodynamic behavior at very high density. This has already been observed in Lubo (2000); Kalyana Rama (2001) where some cosmological consequences have been pointed out. However, a complete analysis of this theory in the context of VSL cosmology has not been done and shall be provided in this work. The most relevant features of the minimal length uncertainty theory are summarized in Sec. II, and are used to derive the expressions for the energy density and pressure of a radiation fluid in Sec. III. In Sec. IV, the cosmological implications are discussed and a further outlook is presented in Sec. V.

## II  Quantum theory with minimal short distance uncertainty

An explicit example for the choice of commutation relations in three spatial dimensions that maintains translation and rotation invariance is given in Kempf and Mangano (1997):

$`[𝐱_i,𝐩_j]`$ $`=`$ $`i\mathrm{}\left({\displaystyle \frac{\beta 𝐩^2}{(1+2\beta 𝐩^2)^{1/2}1}}\delta _{ij}+\beta 𝐩_i𝐩_j\right)`$ (3)
$`[𝐱_i,𝐱_j]`$ $`=`$ $`0`$ (4)
$`[𝐩_i,𝐩_j]`$ $`=`$ $`0`$ (5)
Demanding translation and rotation invariance, this choice is unique to first order in $`\beta `$ 111Note that in Lubo (2000), a different choice for the commutation relations was made which is noncommutative in the spatial coordinates and breaks translation invariance.. It does, however, violate Lorentz invariance and thereby specify a preferred frame which, in the present context, coincides with the cosmological rest frame.

Owing to the finite short distance uncertainty, the usual Hilbert space representation in terms of positions eigenstates is no longer available (in contrast with the momentum representation $`|p`$, which still is). An alternative representation is found by introducing the translators $`𝐓_i`$ with eigenstates $`|\kappa `$ obeying $`𝐓_i|\kappa =\kappa _i/i\mathrm{}|\kappa `$, $`i=1\mathrm{}3`$, and $`\kappa ^2<2/\beta `$ Kempf and Mangano (1997). The projection onto momentum space is given by

$$\kappa |p=\delta \left(p_i\frac{\kappa _i}{1\beta \kappa ^2/2}\right).$$
(6)
In terms of $`\psi (\kappa )=\kappa |\psi `$, the operator representations and the scalar product become:

$`𝐱_i\psi (\kappa )`$ $`=`$ $`i\mathrm{}_{\kappa _i}\psi (\kappa )`$ (7)
$`𝐩_i\psi (\kappa )`$ $`=`$ $`{\displaystyle \frac{\kappa _i}{1\beta \kappa ^2/2}}\psi (\kappa )`$ (8)
$`\psi _1|\psi _2`$ $`=`$ $`{\displaystyle \underset{\kappa ^2<2/\beta }{}}\mathrm{d}^3\kappa \psi _1^{}(\kappa )\psi _2(\kappa ).`$ (9)
The same representation was employed in Kempf (2001); Kempf and Niemeyer (2001).

As demonstrated in Kempf and Mangano (1997), it is now possible to construct the closest analogue to the position representation by replacing the usual position eigenstates with states of maximal localization around a given position, $`|x^{\mathrm{ml}}`$. Their $`\kappa `$-space representation can be constructed by applying the translation operator $`e^{x𝐓}`$ to the maximally localized field around the origin, and therefore it varies as $`e^{ix\kappa /\mathrm{}}`$. Using (6), one finds the quasi-position representation of the plane wave with momentum $`p`$:

$$p|x^{\mathrm{ml}}=N(p^2)\mathrm{exp}\left(\frac{ix}{\mathrm{}}\frac{p\sqrt{1+2\beta p^2}1}{\beta p^2}\right),$$
(10)
where $`N(p^2)`$ is a normalization coefficient. It is now clear that $`\kappa /\mathrm{}`$ acts as a substitute of the usual wavenumber $`k`$, with the crucial difference that it only allows a minimum wavelength $`\lambda _{\mathrm{min}}=\pi \mathrm{}(2\beta )^{1/2}`$ for $`p\mathrm{}`$, as required at the onset. This will become important for the computation of the thermodynamic variables in the next section.

## III  Equation of state of a radiation fluid

In order to derive the energy density and pressure of an ideal gas of photons (with straightforward generalization to other ultrarelativistic particles), the grand canonical formalism can be used (e.g. Reichl (1980)). The grand potential is defined as

$$\mathrm{\Omega }=k_\mathrm{B}T\mathrm{ln}𝒵,$$
(11)
where $`𝒵`$ is the grand partition function,

$$𝒵=\mathrm{Tr}e^{(𝐇\mu 𝐍)/k_\mathrm{B}T},$$
(12)
$`𝐇`$ is the Hamiltonian, $`𝐍`$ is the number operator, $`k_\mathrm{B}`$ is Boltzmann’s constant, and $`T`$ is the temperature. Evaluating the trace in the number representation for bosons, one finds as usual

$$\mathrm{\Omega }=k_\mathrm{B}T\underset{l}{}\mathrm{ln}\left(1e^{(E^l\mu )/k_\mathrm{B}T}\right).$$
(13)
The summation is over all allowed momentum states $`p^l`$ and $`E^l=c|p^l|`$ is the energy of a relativistic particle in this state.

At this point, one usually transforms the sum into an integral over the three-dimensional wavenumber $`k`$ by partitioning space into boxes of volume $`V=L^3`$ and demanding $`L`$-periodicity of plane waves with momentum $`p=\mathrm{}k`$ in the position representation. In the present framework, a position representation is unavailable. Instead, the continuum limit is now most conveniently taken in $`\kappa `$-space: requiring the periodicity of $`e^{ix\kappa /\mathrm{}}`$ yields $`\kappa =2\pi \mathrm{}l/L`$, where $`lM`$ labels the discrete quantum state and $`M=\{l^n|l^2<(L/\lambda _{\mathrm{min}})^2\}`$. For sufficiently large $`V`$, the energy levels are nearly continuous and Eq. (13) can be written as

$$\mathrm{\Omega }=\frac{gVk_\mathrm{B}T}{(2\pi \mathrm{})^3}\underset{\kappa ^2<2/\beta }{}\mathrm{d}^3\kappa \mathrm{ln}\left(1e^{[E(\kappa )\mu ]/k_\mathrm{B}T}\right)$$
(14)
with $`E(\kappa )=c|p(\kappa )|`$ and $`g=2`$ for photons. After transforming the integral to energy space and writing it in terms of the dimensionless variables $`ϵ=E/E_\beta `$, $`\stackrel{~}{\mu }=\mu /E_\beta `$, and $`\tau =k_\mathrm{B}T/E_\beta `$, where $`E_\beta =c\beta ^{1/2}`$ is the characteristic energy scale where quantum gravitational effects become relevant, one obtains

$$\mathrm{\Omega }=\frac{4\pi gVE_\beta ^4\tau }{(2\pi \mathrm{}c)^3}\underset{0}{\overset{\mathrm{}}{}}dϵJ(ϵ)ϵ^2\mathrm{ln}\left(1e^{[ϵ\stackrel{~}{\mu }]/\tau }\right).$$
(15)
Here, $`J`$ is the Jacobian determinant of the transformation $`\kappa _ip_i=\kappa _i(1\beta \kappa ^2/2)^1`$, $`p_i`$, given in terms of $`ϵ`$ by

$$J(ϵ)=\frac{2}{ϵ^6}\left(2+ϵ^2\frac{2+3ϵ^2}{\sqrt{1+2ϵ^2}}\right).$$
(16)
$`J(ϵ)`$ essentially contains all of the information about the modification of the phase space volume that results from the short distance cutoff.

One can now derive the thermodynamic variables from the grand potential in the standard way. Using $`\rho (\tau )=V^1_{\tau ^1}\tau ^1\mathrm{\Omega }`$ for the energy density, $`P(\tau )=_V\mathrm{\Omega }`$ for the pressure, $`N(\tau )=V^1_\mu \mathrm{\Omega }`$ for the number density and setting $`\stackrel{~}{\mu }=0`$ for ultrarelativistic particles, one finds:

$`\rho (\tau )`$ $`=`$ $`{\displaystyle \frac{4\pi gE_\beta ^4}{(2\pi \mathrm{}c)^3}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dϵJ(ϵ){\displaystyle \frac{ϵ^3}{e^{ϵ/\tau }1}}`$ (17)
$`P(\tau )`$ $`=`$ $`{\displaystyle \frac{4\pi gE_\beta ^4\tau }{(2\pi \mathrm{}c)^3}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dϵJ(ϵ)ϵ^2\mathrm{ln}\left(1e^{ϵ/\tau }\right)`$ (18)
$`N(\tau )`$ $`=`$ $`{\displaystyle \frac{4\pi gE_\beta ^3}{(2\pi \mathrm{}c)^3}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dϵJ(ϵ){\displaystyle \frac{ϵ^2}{e^{ϵ/\tau }1}}`$ (19)
It is useful to study the asymptotic behavior of the integrals for very small and very large $`\tau `$ (e.g., Hinch (1991)). In the case of $`\rho (\tau )`$, rescaling the integrand by substituting $`u=ϵ/\tau `$ shows that the dominant contribution comes from $`u=O(1)`$ for $`\tau 0`$. One can therefore write $`J=1+O(u^2\tau ^2)`$ and find the usual result $`\rho (\tau )\tau ^4`$ for $`\tau 0`$, as expected. For large $`\tau `$, we first note that $`J(ϵ)`$ drops as $`ϵ^4`$ for large $`ϵ`$. Hence, for $`ϵ=O(1)`$ the integrand is $`O(\tau )`$ and so is the integral, whereas for $`ϵ=O(\tau )`$, the integrand is $`O(\tau ^1)`$ and the integral is $`O(1)`$. The dominant contribution therefore comes from $`ϵ=O(1)`$ so that the exponential can be approximated by $`1+ϵ/\tau +O(ϵ^2\tau ^2)`$, yielding

$$\rho _\tau \mathrm{}(\tau )\frac{4\pi gE_\beta ^4}{(2\pi \mathrm{}c)^3}\tau \underset{0}{\overset{\mathrm{}}{}}dϵJ(ϵ)ϵ^2,$$
(20)
where the integral evaluates to $`2^{3/2}/3`$.

In the case of $`P(\tau )`$, partial integration yields

$$P(\tau )=\frac{4\pi gE_\beta ^4}{3(2\pi \mathrm{}c)^3}\underset{0}{\overset{\mathrm{}}{}}dϵG(ϵ)\frac{ϵ^3}{e^{ϵ/\tau }1},$$
(21)
where $`G(ϵ)`$ is defined as

$`G(ϵ)`$ $`=`$ $`{\displaystyle \frac{3}{ϵ^3}}{\displaystyle dϵJ(ϵ)ϵ^2}`$ (22)
$`=`$ $`\sqrt{1+2ϵ^2}J(ϵ).`$
For small $`\tau `$, $`G(ϵ)=1O(ϵ^2)`$ behaves like $`J(ϵ)`$, and one obtains the expected result $`P=\rho /3`$. On the other hand, $`Gϵ^3`$ drops more slowly for $`\tau \mathrm{}`$ than $`J(ϵ)`$. It is only for this reason that the ratio of pressure and energy density, and hence the dynamics of the scale factor, differs from its usual behavior at high temperatures, as will be demonstrated in Sec. IV. Indeed, for $`ϵ=O(\tau )`$, the integrand is now $`O(1)`$ and the integral is $`O(\tau )`$, contributing equally strongly to the integral as $`ϵ=O(1)`$. In this case, a good approximation is given by

$`P_\tau \mathrm{}(\tau )`$ $``$ $`{\displaystyle \underset{O(1)}{\overset{O(\tau )}{}}}dϵ{\displaystyle \frac{\tau }{ϵ}}`$ (23)
$`=`$ $`\tau (A+B\mathrm{log}\tau ).`$
The cosmological evolution of $`\rho `$ is determined, via the energy conservation equation, by the equation of state parameter $`w=P/\rho `$. For small $`\tau `$, one recovers the usual equation of state, $`w=1/3`$, whereas $`\tau \mathrm{}`$ gives $`w_\tau \mathrm{}(\rho )=C+D\mathrm{log}\rho _\tau \mathrm{}`$. The same form of $`w`$ was discovered and analyzed in Alexander and Magueijo (2001), hence one can expect a similar result for the cosmological evolution. As shown in Sec. IV, this is indeed true. The result for $`w(\rho )`$ was confirmed numerically by inverting $`\rho (\tau )`$ in order to evaluate $`P(\tau )`$. A good high-density fit is provided by $`C1.2`$ and $`D2.3`$.

Finally, using the same arguments as for $`\rho (\tau )`$, the expression for $`N(\tau )`$ yields the standard result for $`\tau 0`$ while $`\tau \mathrm{}`$ gives

$$N_\tau \mathrm{}(\tau )\frac{4\pi gE_\beta ^3}{(2\pi \mathrm{}c)^3}\tau \underset{0}{\overset{\mathrm{}}{}}dϵJ(ϵ)ϵ=\frac{2\pi gE_\beta ^3}{(2\pi \mathrm{}c)^3}\tau .$$
(24)
All of the preliminaries are now in place to examine the cosmological impact of the modified equation of state, which will proceed in close analogy with Ref. Alexander and Magueijo (2001).

## IV  Cosmological implications

The evolution of $`\rho (t)`$ and the scale factor $`a(t)`$ in the very early universe are determined by the Friedmann and energy conservation equations with negligible curvature and cosmological constant terms, given by

$`\left({\displaystyle \frac{\dot{a}}{a}}\right)^2`$ $`=`$ $`{\displaystyle \frac{\rho }{3M_{\mathrm{Pl}}^2}}`$ (25)
$`{\displaystyle \frac{\dot{\rho }}{\rho }}`$ $`=`$ $`{\displaystyle \frac{3\dot{a}}{a}}\left(1+w(\rho )\right),`$ (26)
where a dot corresponds to a derivative with respect to cosmic time and $`M_{\mathrm{Pl}}=(8\pi G)^{1/2}`$ is the reduced Planck mass. These equations can only be solved numerically (starting with the known low-temperature solution and integrating backward), using the numerical evaluation of $`w(\rho )`$. However, if $`w(\rho )`$ varies sufficiently slowly with $`\rho `$, one can approximate the solution for the scale factor by the usual one for constant $`w`$, with $`w`$ replaced by $`w(\rho )`$ and an offset for the origin of the time coordinate $`t_0`$:

$$a(t)(tt_0)^{2/3(1+w(\rho ))}.$$
(27)
Fig. 1 demonstrates that this is a good approximation at all times except for the brief period where $`w`$ transitions from $`1/3`$ to the logarithmic high-temperature behavior.

Eq. (27) captures the relevant scale factor dynamics in our model: at very high density, the deceleration $`\ddot{a}/aH^2`$ is very large, giving rise to a smaller than usual Hubble rate $`H=\dot{a}/a`$, whereas the standard scale factor evolution is recovered at densities far below the string scale.

The high-density evolution of $`\rho `$ as a function of the scale factor can be obtained directly by solving Eq. (26) with $`w(\rho )w_\tau \mathrm{}`$, yielding

$$\mathrm{log}\rho _\tau \mathrm{}(a)=Ea^{3D}\frac{C+1}{D}.$$
(28)
This was also shown in Alexander and Magueijo (2001).

The question whether the cosmological horizon problem can be solved in this model hinges crucially on the choice for the effective propagation speed of information, $`c_{\mathrm{eff}}`$. The VSL mechanism works if the comoving horizon scale,

$$r_\mathrm{h}=\frac{c_{\mathrm{eff}}}{aH}=\frac{c_{\mathrm{eff}}}{\dot{a}},$$
(29)
declines for some period of time before turning around to its current growing behavior.

A case can be made that the most suitable measure for $`c_{\mathrm{eff}}`$ is given by the speed of sound, $`c_\mathrm{s}`$, of the radiation fluid rather than by the group or phase velocity of individual photons. It is given by

$$c_\mathrm{s}^2=\left(\frac{P}{\rho }\right)_{S=const}=w+\rho \frac{w}{\rho }.$$
(30)
In the high-density limit, $`c_\mathrm{s}`$ scales like $`(\mathrm{log}\rho )^{1/2}a^{3/2}`$, whereas $`\dot{a}a\rho ^{1/2}a\mathrm{exp}(1/2a^3)`$ which clearly wins as $`a0`$. Consequently, $`r_\mathrm{h}`$ is always a growing function of time in this case and the horizon problem remains.

The group and phase velocities are given by

$`c_\mathrm{g}(ϵ)={\displaystyle \frac{E}{|\kappa |}}`$ $`=`$ $`{\displaystyle \frac{c}{2}}\left(1+2ϵ^2+\sqrt{1+2ϵ^2}\right)`$ (31)
$`c_{\mathrm{ph}}(ϵ)={\displaystyle \frac{E}{|\kappa |}}`$ $`=`$ $`{\displaystyle \frac{c}{2}}\left(1+\sqrt{1+2ϵ^2}\right).`$ (32)
As pointed out in Alexander and Magueijo (2001), the evaluation of $`ϵ`$ in these expressions is somewhat ambiguous. The bulk of all photons will populate a region around the peak in the modified Planck distribution, which saturates at $`\tau 1`$. After saturation, the group and phase velocities of these photons will cease growing as a function of increasing density. However, the authors of Alexander and Magueijo (2001) argue that the “fast tail” of the Planck distribution may provide the required causal contact, and they propose to evaluate $`c_\mathrm{g}`$ and $`c_{\mathrm{ph}}`$ at $`ϵ\tau \rho _\tau \mathrm{}`$. In this case, $`c_\mathrm{g}\rho ^2`$ and $`c_{\mathrm{ph}}\rho `$ at high densities, so that $`r_\mathrm{h}`$ indeed declines as a function of time (consistent with Alexander and Magueijo (2001)). Consequently, this choice of effective velocity solves the horizon problem, provided that the model assumptions made in Sec. II remain valid well into the regime $`\tau 1`$.

As a final possibility, one can consider a thermal average of $`c_\mathrm{g}`$ (or $`c_{\mathrm{ph}}`$), defined for instance as the particle number weighed value:

$`c_{\mathrm{g},\mathrm{ph}}(\tau )`$ $`=`$ $`N^1(\tau ){\displaystyle \frac{4\pi gE_\beta ^3}{(2\pi \mathrm{}c)^3}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dϵ{\displaystyle \frac{K_{\mathrm{g},\mathrm{ph}}(ϵ)ϵ^2}{e^{ϵ/\tau }1}}`$ (33)
$`K_{\mathrm{g},\mathrm{ph}}(ϵ)`$ $`=`$ $`c_{\mathrm{g},\mathrm{ph}}J(ϵ).`$ (34)
The asymptotic properties of $`N(\tau )`$ were discussed in Sec. III. Noting that $`K_{\mathrm{g},\mathrm{ph}}=1O(ϵ^2)`$ for small $`ϵ`$ one finds $`c_{\mathrm{g},\mathrm{ph}}(\tau )c`$ as required. On the other hand, $`K_\mathrm{g}ϵ^2`$ for large $`ϵ`$, whereas $`K_{\mathrm{ph}}ϵ^3`$. The situation for the integral over $`K_\mathrm{g}`$ at $`\tau \mathrm{}`$ is therefore analogous to that of $`P_\tau \mathrm{}`$, yielding the same logarithmic growth. The normalization by $`N(\tau )`$ gives rise to the same behavior as found earlier for $`w_\tau \mathrm{}`$, i.e. $`c_\mathrm{g}_\tau \mathrm{}\mathrm{log}\rho `$. In contrast, the integral over $`K_{\mathrm{ph}}`$ at $`\tau \mathrm{}`$ is analogous to that of $`\rho _\tau \mathrm{}`$, hence the normalized result for $`c_{\mathrm{ph}}_\tau \mathrm{}`$ asymptotes to a constant independent of $`\tau `$. Explicitly, one finds $`c_{\mathrm{ph}}_\tau \mathrm{}=2c`$.

To summarize, the thermally averaged group or phase velocities grow at best logarithmically. This is insufficient to solve the horizon problem by means of the VSL mechanism, as argued above for $`c_\mathrm{s}`$.

Fig. 2 shows the size of the comoving horizon as a function of time, computed numerically using the “fast tail” group velocity and the thermally averaged one.

## V  Summary and Outlook

The analysis presented in this work is based on a rather simple premise: The spatial separation that can be probed by high-energy particles has a finite minimum value. Without referring to any specific quantum gravitational effect in particular, it is a very general model for the physics of relativistic particles at densities that occurred in the very early universe. It is very interesting that such a simple modification may provide an alternative solution to the horizon problem without resorting to the inflationary paradigm. Of course, one must now think of possibilities to resolve other quandaries that inflation has proved to be so helpful with, such as the flatness and relics problems and, perhaps most importantly, the generation of scale invariant fluctuations. Some ideas in this direction have already been put forward in the context of VSL Albrecht and Magueijo (1999); Alexander and Magueijo (2001), and they can be readily generalized to the model in this work.

As always, many open questions remain. One of them, which may be answered by a more detailed analysis of transport processes in the early universe, is the correct choice of the effective communication speed. On a more fundamental level, one may question the validity of the minimal distance uncertainty principle as early as may be required to solve the horizon problem. In other words, it is possible (and even likely) that additional, quantum gravitational degrees of freedom will be excited in the relevant regime of $`\rho `$ that will drastically change the equation of state. Unfortunately, the answer to this question might remain elusive until the fundamental theory is known in full detail.

###### Acknowledgements.

I am happy to thank Achim Kempf and Joao Magueijo for valuable comments.


text:
# Ehrenfest times for classically chaotic systems

## Abstract

We describe the quantum mechanical spreading of a Gaussian wave packet by means of the semiclassical WKB approximation of Berry and Balazs. We find that the time scale $`\tau `$ on which this approximation breaks down in a chaotic system is larger than the Ehrenfest times considered previously. In one dimension $`\tau =\frac{7}{6}\lambda ^1\mathrm{ln}(A/\mathrm{})`$, with $`\lambda `$ the Lyapunov exponent and $`A`$ a typical classical action.

today

According to Ehrenfest’s theorem , the propagation of a quantum mechanical wave packet is described for short times by classical equations of motion. The time scale at which this correspondence between quantum and classical dynamics breaks down is called the Ehrenfest time. If the classical dynamics is chaotic with Lyapunov exponent $`\lambda `$, then the Ehrenfest time $`\tau `$ is of order $`\lambda ^1\mathrm{ln}(A/\mathrm{})`$ (with $`A`$ a typical classical action of the dynamical system) . There is actually more than a single Ehrenfest time, corresponding to different types of semiclassical approximations. Although they differ only by a numerical coefficient, $`\tau _i=c_i\lambda ^1\mathrm{ln}(A/\mathrm{})`$, the structure of the wave function changes qualitatively from one time scale to the next.

Up to a time $`\tau _1`$, with $`c_1=1/6`$, the initial coherent state will retain its Gaussian form with vanishing error in the limit $`\mathrm{}0`$ . For longer times up to $`\tau _2`$, with $`c_2=1/2`$, the uncertainty in position and momentum of the particle remains small but the phase space structure of the wave packet deviates strongly from a Gaussian. For times greater than $`\tau _2`$ the wave function no longer has the form of a wave packet (this is the “mixing regime” of Refs. ), but up to a time $`\tau _3`$ it can still be described semiclassically by the time-dependent WKB approximation of Berry and Balazs . As we will show in this paper, the WKB representation implies $`c_3=7/6`$ for a single degree of freedom (with simple generalizations for higher dimensions). This is larger than the value $`c_3=2/3`$ obtained by Bouzouina and Robert  from a different semiclassical approximation.

Let us start with the Gaussian one-dimensional wave packet

$$\mathrm{\Psi }(x)=\left(\frac{\alpha }{\pi \mathrm{}}\right)^{1/4}\mathrm{exp}\left(i\frac{p_0x}{\mathrm{}}+(i\beta \alpha )\frac{(xx_0)^2}{2\mathrm{}}\right).$$
(1)
Initially $`\beta (t=0)=0`$ and $`\alpha (t=0)=p_F/L`$, where $`p_F`$ and $`L`$ are the typical classical momentum and length. The typical classical action is $`A=p_FL`$. The parameters $`x_0(t),p_0(t)`$ follow the classical trajectory for $`\mathrm{}A`$. We will measure momentum and coordinate in units of $`p_F`$ and $`L`$, respectively, so that $`\alpha (0)=1`$ and $`A=1`$. For chaotic dynamics with Lyapunov exponent $`\lambda `$ one has $`\alpha (t)\mathrm{exp}(2\lambda t)`$, hence $`\alpha 1`$ for $`t1/\lambda `$.

To describe the time evolution in phase space we consider the Wigner function

$`W(x,p)={\displaystyle \mathrm{\Psi }(x+\frac{y}{2})\mathrm{\Psi }^{}(x\frac{y}{2})e^{ipy/\mathrm{}}\frac{dy}{2\pi \mathrm{}}}`$ (2)
$`={\displaystyle \frac{1}{\pi \mathrm{}}}\mathrm{exp}\left({\displaystyle \frac{\alpha (xx_0)^2}{\mathrm{}}}{\displaystyle \frac{[pp_0\beta (xx_0)]^2}{\alpha \mathrm{}}}\right).`$ (3)
The wave packet is centered at $`x_0(t),p_0(t)`$ and for $`\alpha (t)1`$ becomes highly elongated and tilted with slope $`\mathrm{\Delta }p/\mathrm{\Delta }x\beta `$. It has length $`l_{}=\sqrt{\mathrm{}(1+\beta ^2)/\alpha }`$ and width $`l_{}=\sqrt{\mathrm{}\alpha /(1+\beta ^2)}`$, so that the area in phase space is conserved exactly, $`l_{}l_{}=\mathrm{}`$. The Gaussian quantum wave packet satisfies the classical Liouville theorem.

The Gaussian form (1) takes into account the elongation of the wave packet, but not the curvature that develops in time and results in a bending of the packet. To describe the curvature we add an imaginary cubic term in the exponent in Eq. (1),

$$\mathrm{\Psi }(x)=\left(\frac{\alpha }{\pi \mathrm{}}\right)^{1/4}\mathrm{exp}\left(i\frac{p_0x}{\mathrm{}}+\frac{(i\beta \alpha )x^2}{2\mathrm{}}+i\frac{\gamma x^3}{3\mathrm{}}\right).$$
(4)
(For simplicity we have put $`x_0=0`$.) The cubic term leads to an appreciable phase shift over a length $`l_{}(\mathrm{}/\alpha )^{1/2}`$ when $`(\gamma /\mathrm{})(\mathrm{}/\alpha )^{3/2}1`$, hence when $`\alpha (t)\mathrm{}^{1/3}\gamma ^{2/3}`$.

For $`\alpha \mathrm{}^{1/3}\gamma ^{2/3}`$ the Wigner function takes again a simple form, in terms of the Airy function Ai:

$$W(x,p)=\frac{\alpha ^{1/2}\mathrm{exp}(\alpha x^2/\mathrm{})}{\pi \mathrm{}^{1/2}(\gamma \mathrm{}^2/4)^{1/3}}\text{Ai}\left(\frac{p_0+\beta x+\gamma x^2p}{(\gamma \mathrm{}^2/4)^{1/3}}\right).$$
(5)
One can check that $`W(x,p)\delta (x)\delta (pp_0)`$ when $`\mathrm{}0`$ (at fixed $`\alpha `$), by means of the identity $`lim_{\epsilon 0}\text{Ai}(z/\epsilon )/\epsilon =\sqrt{\pi }\delta (z)`$. At finite $`\mathrm{}`$ the wave packet is extended along the curved line $`p=p_0+\beta x+\gamma x^2`$. Since $`p,p_0,x`$ are of order unity, the two parameters $`\beta `$ and $`\gamma `$ are of order unity as well (in contrast to $`\alpha `$, which is $`1`$). The transverse width is of order

$$l_{}\gamma ^{1/3}\mathrm{}^{2/3}(1+\beta ^2)^{1/2}.$$
(6)
The length of the packet remains at $`l_{}\sqrt{\mathrm{}(1+\beta ^2)/\alpha }`$. Since now $`l_{}l_{}\mathrm{}`$, the Liouville theorem no longer holds.

To obtain the Ehrenfest time we parameterize time as

$$t=\frac{c}{\lambda }\mathrm{ln}\frac{1}{\mathrm{}}.$$
(7)
The classical limit for a chaotic system means $`\mathrm{}0`$, $`t\mathrm{}`$ at fixed $`c`$. Different coefficients $`c`$ follow from different semiclassical approximations. If we use the Gaussian wave packet (1), without the cubic term to account for the curvature, then we need $`\alpha (t)\mathrm{}^{1/3}\gamma ^{2/3}`$. Since $`\alpha e^{2\lambda t}\mathrm{}^{2c}`$ we need $`c<1/6`$. The upper limit of $`c`$ gives the first Ehrenfest time $`\tau _1=\frac{1}{6}\lambda ^1\mathrm{ln}(1/\mathrm{})`$.

The classical limit can be reached for longer times if we use the wave packet (4), including the cubic term. The dimensions of the packet for $`t>\tau _1`$ scale with $`\mathrm{}`$ as

$$l_{}\mathrm{}^{2/3},l_{}\mathrm{}^{1/2c}.$$
(8)
For $`c<1/2`$ the length of the packet approaches zero in the classical limit. This upper limit of $`c`$ gives the second Ehrenfest time $`\tau _2=\frac{1}{2}\lambda ^1\mathrm{ln}(1/\mathrm{})`$.

For $`t>\tau _2`$ the length of the wave packet exceeds the size of the system and is no longer small compared to the radius of curvature. For these large times we may adopt the semiclassical WKB approximation of Berry and Balazs . Consider a curve in phase space $`p(x)`$ and a phase space distribution $`\rho (p(x),x)`$. Both $`p`$ and $`\rho `$ evolve in accordance with classical equations of motion. For $`t>\tau _2`$ the function $`p(x)`$ is multivalued with an exponentially large number of branches $`\mathrm{exp}[\lambda (t\tau _2)]`$. The quantum wave function in this “mixing” regime has the form

$$\mathrm{\Psi }(x)=\underset{k}{}f_k(x)\mathrm{exp}\left(i\sigma _k(x)/\mathrm{}\right).$$
(9)
The summation over $`k`$ accounts for the different branches of the multivalued function $`p(x)`$. The two functions $`f`$ and $`\sigma `$ are related for $`\mathrm{}0`$ to $`p`$ and $`\rho `$ by the correspondence principle,

$$\frac{d\sigma }{dx}=p(x),f=\sqrt{\rho (p,x)}.$$
(10)
An explicit description of the evolution of the wave function (9) for quantum maps may be found in Ref. .

Near the point $`x_b`$ at which $`p(x)`$ bifurcates into two branches, one has $`p=p_b\pm a\sqrt{xx_b}`$, $`\rho =b/\sqrt{xx_b}`$. The wave function there is

$$\mathrm{\Psi }=\left(\mathrm{}/a\right)^{1/3}b^{1/2}\text{Ai}\left(\left(a/\mathrm{}\right)^{2/3}(xx_b)\right)e^{ip_bx},$$
(11)
up to an overall phase. The phase difference between the bifurcation points can be determined from Eqs. (9) and (10). Because the curve $`p(x)`$ is not closed, there is no analog of the Bohr-Sommerfeld quantization rule.

The Wigner function corresponding to the wave function (9), being quadratic in $`\mathrm{\Psi }`$, contains both diagonal ($`W_{kk}|f_k|^2`$) and oscillating nondiagonal ($`W_{km}f_k^{}f_m`$) contributions. Far from bifurcations the diagonal contributions to the Wigner function read

$`W_{kk}(x,p)={\displaystyle \mathrm{exp}\left(\frac{iy(\sigma ^{}p)}{\mathrm{}}+\frac{iy^3\sigma ^{\prime \prime \prime }}{24\mathrm{}}\right)\frac{|f(x)|^2dy}{2\pi \mathrm{}}}`$ (13)
$`={\displaystyle \frac{2}{\sqrt{\pi }}}\left({\displaystyle \frac{1}{\mathrm{}^2\sigma ^{\prime \prime \prime }}}\right)^{1/3}|f(x)|^2\text{Ai}\left({\displaystyle \frac{2(\sigma ^{}p)}{(\mathrm{}^2\sigma ^{\prime \prime \prime })^{1/3}}}\right).`$
We have made a Taylor expansion of $`\sigma (x\pm y/2)`$ and neglected the difference between $`f(x\pm y/2)`$ and $`f(x)`$.

If we parameterize time as in Eq. (7) we have for both $`l_{}`$ and $`l_{}`$ the same scaling with $`\mathrm{}`$ as in Eq. (8). The range of validity of Eq. (9) is limited by the condition that different branches should be distinguishable. This requires that the different parts of the curve $`p(x)`$ in phase space should not get closer than $`ł_{}`$. Their spacing is of order $`1/l_{}`$ (assuming a uniform filling of phase space), hence

$$l_{}l_{}1\mathrm{}^{7/6c}1.$$
(14)
The upper limit of $`7/6`$ for $`c`$ leads to the third Ehrenfest time,

$$\tau _3=\frac{7}{6\lambda }\mathrm{ln}\frac{1}{\mathrm{}}.$$
(15)
The third derivative $`\sigma ^{\prime \prime \prime }`$ in Eq. (13) vanishes at the points of inflection of the curve $`p(x)`$. In order to find the Wigner function there one should expand $`\sigma (x\pm y/2)`$ up to terms of order $`y^5`$. This leads to a different scaling $`l_{}\mathrm{}^{4/5}`$ of the width of the Wigner function near the inflection points. Because these are isolated points, they will not contribute to the matrix elements of nonsingular operators (containing only smooth functions of $`x`$ and $`p`$). This different scaling should therefore not affect the Ehrenfest time (15).

The nondiagonal contributions $`W_{km}`$ to the Wigner function lead to the “ghost curves” discussed in Ref. . The Wigner function near these curves is given by the same Airy function as in Eq. (13), but in addition acquires a strongly oscillating factor. Due to these oscillations the nondiagonal terms do not contribute to the matrix elements of nonsingular operators. (They may play a role in the decoherence by the environment .) At $`t\tau _3`$ the ghost curves merge with the (multivalued) curve $`p(x)`$ and become indistinguishable.

The time scale (15) for the breakdown of the WKB approximation is greater than the Ehrenfest time $`\frac{2}{3}\lambda ^1\mathrm{ln}(1/\mathrm{})`$ in the mixed regime obtained in Ref. . That shorter time scale may signal the breakdown of the series expansion $`\sigma _k(x)_{j=0}\sigma _{kj}(x)\mathrm{}^j`$. Then Eq. (10) would no longer hold, but for $`t<\tau _3`$ the representation (9) with a renormalized function $`\sigma _k(x)`$ would still be valid.

So far we have discussed a one-dimensional (1D) chaotic system, which in general can be represented by an area preserving map . A familiar example is the kicked rotator . For mesoscopic quantum dots, however, a more relevant model is the $`d`$-dimensional ($`d=2,3`$) Schrödinger equation with a smooth potential $`V(\stackrel{}{r})`$. The Gaussian wave packet then takes the form

$$\mathrm{\Psi }(\stackrel{}{r})\mathrm{exp}\left[\frac{i}{\mathrm{}}\left(S(\stackrel{}{r}_0(t))+\stackrel{}{p}_0\stackrel{}{x}+\frac{\zeta _{ln}}{2}x_lx_n\right)\right].$$
(16)
Here $`S`$ is the action for the classical trajectory $`\stackrel{}{r}_0(t)`$ and we have defined $`\stackrel{}{p}_0=m\dot{\stackrel{}{r}}_0`$, $`\stackrel{}{x}=\stackrel{}{r}\stackrel{}{r}_0,\zeta _{ln}=\beta _{ln}+i\alpha _{ln}`$. As before, we rescale momentum and coordinate such that the typical classical action $`A=1`$. Initially, $`\zeta _{ln}i\delta _{ln}`$. Similarly to the one-dimensional case, $`\alpha _{ln}`$ defines the form of the packet in coordinate space and $`\beta _{ln}=\mathrm{\Delta }p_l/\mathrm{\Delta }x_n`$ give the angles in phase space. Substituting the wave function (16) into the Schrödinger equation one finds Newton’s equation of motion for $`\stackrel{}{r}_0`$. The spreading of the wave packet in phase space is described by

$$\dot{\zeta }_{ln}=\frac{1}{m}\zeta _{lk}\zeta _{kn}+\frac{^2V}{r_lr_n}|_{\stackrel{}{r}=\stackrel{}{r}_0}.$$
(17)
This is the equation describing the spreading in phase space of a small Gaussian bunch of classical particles.

The Wigner function corresponding to the wave function (16) has the Gaussian form $`W\mathrm{exp}(Q_lM_{ln}Q_n/\mathrm{})`$, where $`\stackrel{}{Q}=(\stackrel{}{r}\stackrel{}{r}_0,\stackrel{}{p}\stackrel{}{p}_0)`$ is a vector in $`2d`$ dimensional phase space. The $`d`$ Lyapunov exponents $`\lambda _i`$ ($`i=1,2,\mathrm{}d`$) govern the large-time behavior of the eigenvalues $`m_i=1/m_{2di+1}\mathrm{exp}(2\lambda _it)`$ of the real symmetric matrix $`M`$. Because of energy conservation one Lyapunov exponent vanishes. We order the $`\lambda `$’s from large to small, so that $`\lambda _1`$ is the largest and $`\lambda _d=0`$.

The wave packet remains Gaussian (preserving the volume $`\mathrm{}^d`$ in phase space) until the curvature starts to play a role (via a cubic term in the action). The corresponding Ehrenfest time $`\tau _1=\frac{1}{6}\lambda _1^1\mathrm{ln}(1/\mathrm{})`$ is the same as in 1D, only now it is defined through the largest Lyapunov exponent $`\lambda _1`$. The second Ehrenfest time, when the length of the packet exceeds the size of the system, also has the same form $`\tau _2=\frac{1}{2}\lambda _1^1\mathrm{ln}(1/\mathrm{})`$.

The third time $`\tau _3`$ is different for $`d=2,3`$ from the 1D case. Instead of Eq. (8) one now has

$$l_{}^{(i)}\mathrm{}^{2/3},l_{}^{(i)}\mathrm{}^{1/2}e^{\lambda _it},i=1,2,\mathrm{}d1.$$
(18)
The longitudinal dimensions $`l_{}^{(i)}`$ correspond to eigenvalues $`m_i`$ with $`1id1`$, and the transverse dimensions $`l_{}^{(i)}`$ to $`m_i`$ with $`d+2i2d`$. The two unit eigenvalues $`m_d=m_{d+1}=1`$ contribute another factor $`\sqrt{\mathrm{}}`$ each to the total volume $`𝒱`$ in phase space covered by the wave packet:

$$𝒱=\mathrm{}\underset{i=1}{\overset{d1}{}}l_{}^{(i)}l_{}^{(i)}\mathrm{}^{7d/61/6}e^{\lambda _{\mathrm{tot}}t},\lambda _{\mathrm{tot}}=\underset{i=1}{\overset{d1}{}}\lambda _i.$$
(19)
The available area $`𝒱_{\mathrm{max}}`$ is restricted to a shell of constant energy with thickness $`\sqrt{\mathrm{}}`$, hence $`𝒱_{\mathrm{max}}\sqrt{\mathrm{}}`$. We require $`𝒱𝒱_{\mathrm{max}}`$ for the semiclassical approximation, which leads to the Ehrenfest time

$$\tau _3=\frac{7d4}{6\lambda _{\mathrm{tot}}}\mathrm{ln}\frac{A}{\mathrm{}},d2.$$
(20)
In conclusion, we examined different time scales $`\tau _i=c_i\lambda ^1\mathrm{ln}(1/\mathrm{})`$ for the breakdown of different types of semiclassical approximations. These Ehrenfest times differ only by a numerical coefficient $`c_i`$, which may seem insignificant. However, this difference is actually a signal of a different power law scaling with $`\mathrm{}`$ of the volume $`𝒱`$ in phase space covered by the wave packet. For short times Liouville’s theorem dictates $`𝒱\mathrm{}`$. For long times \[parameterized as $`t=(c/\lambda )\mathrm{ln}(1/\mathrm{})`$\] the WKB approximation gives $`𝒱\mathrm{}^{7/6c}`$ for a one-dimensional quantum map (such as the kicked rotator) and $`𝒱\mathrm{}^{7d/61/6c}`$ for a $`d`$-dimensional conservative system. These different power laws reflect the fundamental change in the structure of the wave function with increasing time and should therefore have observable consequences. Two possible applications are the Loschmidt echo  and the quantum shot noise , where the Ehrenfest time plays a key role.

This work was supported by the Dutch Science Foundation NWO/FOM and by the National Science Foundation under Grant No. PHY99-07949.