format: text

text:
# Hyperbolic polynomials and spectral order

## Introduction and Main Results

Given a complex polynomial $`P`$ of degree $`n`$ we define $`Z(P)`$ to be the unordered $`n`$-tuple consisting of the zeros of $`P`$, where each zero occurs as many times as its multiplicity. We denote by $`\mathrm{}Z(P)`$ the (unordered) $`n`$-tuple consisting of the real parts of the points in $`Z(P)`$. The polynomial $`P`$ is said to be hyperbolic if all its zeros are real. Note that in this case $`\mathrm{}Z(P)=Z(P)`$. A hyperbolic polynomial whose zeros are simple is called strictly hyperbolic.

The main purpose of this paper is to study the behaviour of the $`n`$-tuple $`Z(P)`$ under the action of certain semigroups of differential operators. For this we shall use the following fundamental result from the theory of stochastic majorizations:

###### Theorem 1.

Let $`X=(x_1,x_2,\mathrm{},x_n)^t`$ and $`Y=(y_1,y_2,\mathrm{},y_n)^t`$ be two unordered $`n`$-tuples of vectors in $`^k`$. The following conditions are equivalent:

1. For any convex function $`f:^k`$ one has $`_{i=1}^nf(x_i)_{i=1}^nf(y_i)`$.
2. There exists a doubly stochastic $`n\times n`$ matrix $`A`$ such that $`\stackrel{~}{X}=A\stackrel{~}{Y}`$, where $`\stackrel{~}{X}`$ and $`\stackrel{~}{Y}`$ are $`n\times k`$ matrices obtained by some (and then any) ordering of the vectors in $`X`$ and $`Y`$.

If the conditions of Theorem 1 are satisfied then we say that $`X`$ is *majorized* by $`Y`$ or $`X`$ is less than $`Y`$ in the spectral order, and write $`XY`$. Theorem 1 is due to Schur and to Hardy, Littlewood and Pólya in the one-dimensional case (\[HLP\]), and to Sherman in the multivariate case (\[S\]). These cases are also known as classical and multivariate majorization, respectively. One can easily check that $`XY`$ implies that $`_{i=1}^nx_i=_{i=1}^ny_i`$.

Let $`_n`$ denote the manifold of monic hyperbolic polynomials of degree $`n`$. We may view $`(_n,)`$ as a partially ordered set, where the ordering relation $``$ is induced by the spectral order on $`n`$-tuples of real numbers (cf. Theorem 1). Thus, if $`P,Q_n`$ then $`PQ`$ if and only if $`Z(P)Z(Q)`$. Note that although the spectral order is only a preordering on $`n`$-tuples of points in $``$, Birkhoff’s theorem (\[MO, Theorem 2.A.2\]) implies that it actually induces a partial ordering on $`_n`$.

Define the following semigroups of differential operators:

$$𝒮=1\lambda \frac{d}{dx}|\lambda ,\stackrel{~}{𝒮}=\left(1\lambda \frac{d}{dx}\right)e^{\lambda \frac{d}{dx}}|\lambda .$$

Note that $`\stackrel{~}{𝒮}`$ is the (largest) subsemigroup of $`𝒮\times e^{\mu \frac{d}{dx}}\mu `$ which consists of operators that preserve the averages of the zeros of polynomials in $`_n`$. The operator $`\left(1\lambda \frac{d}{dx}\right)e^{\lambda \frac{d}{dx}}`$, $`\lambda `$, will be denoted by $`D_\lambda `$ throughout this paper. It follows from the well-known Hermite-Poulain theorem (see \[O\]) that the semigroups $`𝒮`$ and $`\stackrel{~}{𝒮}`$ act on $`_n`$. Our first main result asserts that these semigroups act in fact on $`_n`$ in an order-preserving fashion:

###### Theorem 2.

Let $`P,Q_n`$ be such that $`PQ`$. Then $`P+\lambda P^{}Q+\lambda Q^{}`$ for any $`\lambda `$.

We point out an interesting consequence of Theorem 2:

###### Corollary 1.

If $`P,Q_n`$ are such that $`PQ`$ then $`n^1P^{}n^1Q^{}`$.

The next theorem shows that any polynomial in $`_n`$ is the global minimum of its orbit under the action of the semigroup $`\stackrel{~}{𝒮}`$.

###### Theorem 3.

If $`P_n`$ then $`PD_\lambda P`$ for any $`\lambda `$.

A well-known theorem of Obreschkoff (see \[O\]) states that if $`P`$ and $`Q`$ are real polynomials then the linear pencil of polynomials $`P+\lambda Q`$, $`\lambda `$, consists of hyperbolic polynomials if and only if $`P`$ and $`Q`$ are hyperbolic and those of their zeros which are not common separate each other. The following converse to Theorem 3 shows that real lines in $`_n`$ of the form $`P+\lambda P^{}`$ are characterized by a local minimum property with respect to the partial ordering $``$ on $`_n`$.

###### Theorem 4.

Let $`P_n`$, $`\lambda `$, and let $`Q`$ be a complex polynomial of degree at most $`n1`$. Set $`R_\lambda (x)=P(x+\lambda )\lambda Q(x+\lambda )`$. If $`R_\lambda _n`$ and $`R_0R_\lambda `$ for all small $`\lambda `$ then $`Q=P^{}`$.

We also obtain a generalization of Theorem 3 which shows that real lines in $`_n`$ of the form $`P+\lambda P^{}`$ satisfy in fact a global monotony property:

###### Theorem 5.

If $`P_n`$ and $`\lambda _1,\lambda _2`$ are such that $`\lambda _1\lambda _20`$ and $`|\lambda _1||\lambda _2|`$ then $`D_{\lambda _1}PD_{\lambda _2}P`$.

Finally, we show that real lines in $`_n`$ of the form $`P+\lambda P^{}`$ satisfy an inequality à la Gårding (cf. \[G\]):

###### Theorem 6.

If $`P_n`$ then $`\lambda \mathrm{max}Z(D_\lambda P)`$ is a convex function with a global minimum at $`\lambda =0`$ while $`\lambda \mathrm{min}Z(D_\lambda P)`$ is a concave function with a global maximum at $`\lambda =0`$.

###### Remark 1.

It follows from Theorem 6 that the so-called spread function $`\lambda \mathrm{max}Z(D_\lambda P)\mathrm{min}Z(D_\lambda P)`$ is a convex function with a global minimum at $`\lambda =0`$.

The structure of the paper is as follows: in §1 we sketch the proofs of our main results and in §2 we present further questions and conjectures. The complete proofs will appear elsewhere.

## 1. Outline of the proofs

One of the key ingredients in the proofs of Theorems 2-6 is the following criterion for classical majorization due to Hardy, Littlewood and Pólya (cf. \[HLP\]). We should mention that there are no known analogues of this criterion for multivariate majorization.

###### Theorem 7.

Let $`X=(x_1x_2\mathrm{}x_n)`$ and $`Y=(y_1y_2\mathrm{}y_n)`$ be two $`n`$-tuples of real numbers. Then $`XY`$ if and only if the $`x_i`$’s and the $`y_i`$’s satisfy the following conditions:

1. $`_{i=1}^nx_i=_{i=1}^ny_i`$;
2. $`_{i=0}^kx_{ni}_{i=0}^ky_{ni}`$ for $`0kn2`$.

The proof of Theorem 2 is based on several auxiliary results. Let us first make the following definition:

###### Definition 1.

Let $`P(x)=_{i=1}^n(xx_i)_n`$, $`n2`$, and $`1k<ln`$. Assume that $`x_ix_{i+1}`$, $`1in1`$, and that $`x_kx_l`$. Let further $`t]0,\frac{x_lx_k}{2}]`$ and define $`Q_n`$ to be the polynomial with zeros $`y_i`$, $`1in`$, where $`y_k=x_k+t`$, $`y_l=x_lt`$, and $`y_i=x_i`$, $`ik,l`$. The polynomial $`Q`$ is called the contraction of $`P`$ of type $`(k,l)`$ and coefficient $`t`$. The contraction is called simple if $`l=k+1`$ and it is called non-degenerate if $`t\frac{x_lx_k}{2}`$.

###### Proposition 1.

Let $`P,Q_n`$ be two distinct strictly hyperbolic polynomials such that $`PQ`$. Then there exists a finite sequence $`P_1,\mathrm{},P_m_n`$ such that $`P_1=Q`$, $`P_m=P`$, and $`P_{i+1}`$ is a simple non-degenerate contraction of $`P_i`$, $`1im1`$.

###### Remark 2.

Proposition 1 is true even for polynomials with multiple zeros if the non-degeneracy condition is omitted.

###### Proposition 2.

Theorem 2 is true if $`P`$ and $`Q`$ are strictly hyperbolic polynomials and $`P`$ is a simple (non-degenerate) contraction of $`Q`$.

From Propositions 1 and 2 we deduce that Theorem 2 is true in the generic case when $`P`$ and $`Q`$ have simple zeros. If this is not the case then we let $`x_i`$, $`1in`$, and $`y_i`$, $`1in`$, denote the zeros of $`P`$ and $`Q`$, respectively, and we choose an arbitrary positive number $`\epsilon `$. Let $`P_\epsilon `$ and $`Q_\epsilon `$ be the polynomials with zeros $`x_i(ni)\epsilon `$, $`1in1`$, $`x_n+\frac{n(n1)}{2}\epsilon `$, and $`y_i(ni)\epsilon `$, $`1in1`$, $`y_n+\frac{n(n1)}{2}\epsilon `$, respectively. Note that $`P_\epsilon `$ and $`Q_\epsilon `$ are strictly hyperbolic and that $`P_\epsilon Q_\epsilon `$. The above arguments imply that $`P_\epsilon +\lambda P_\epsilon ^{}Q_\epsilon +\lambda Q_\epsilon ^{}`$ for any $`\lambda `$. Theorem 2 now follows by letting $`\epsilon 0`$.

Let $`P(x)=_{i=1}^n(xx_i)_n`$, $`n2`$, and $`\lambda `$. Assume that $`x_i<x_{i+1}`$, $`1in1`$, and let $`x_i(\lambda )`$, $`1in`$, denote the zeros of $`D_\lambda P`$. If these are labeled so that $`x_i(0)=x_i`$, $`1in`$, then one can show that $`x_i(\lambda )<x_{i+1}(\lambda )`$ and that by varying $`x_n`$ and keeping $`\lambda `$ fixed each $`x_i(\lambda )`$, $`1in1`$, is an increasing function of $`x_n`$. This makes it possible to prove Theorem 3 by induction on $`n`$ in the generic case. If $`P_n`$ has multiple zeros we notice that $`\left(1\epsilon \frac{d}{dx}\right)^{n1}P`$ has simple zeros for any $`\epsilon 0`$. Since Theorem 3 is true for the latter polynomials, we get that Theorem 3 holds in the general case by letting $`\epsilon 0`$.

To prove Theorem 4 we first use Theorem 7 in order to show that if $`P`$ is strictly hyperbolic and $`Q`$ satisfies the assumptions of Theorem 4 then $`Q(x_i)=P^{}(x_i)`$, where $`x_i`$, $`1in`$, are the zeros of $`P`$. Using again the fact that $`\left(1\epsilon \frac{d}{dx}\right)^{n1}P`$ has simple zeros if $`P_n`$ and $`\epsilon 0`$ and also that the operator $`\left(1\epsilon \frac{d}{dx}\right)^{n1}`$ preserves the ordering on $`_n`$ (cf. Theorem 2) we deduce that Theorem 4 holds for any $`P_n`$.

Let now $`P(x)=_{i=1}^n(xx_i)_n`$, $`n2`$, $`\lambda `$, and assume that $`x_i<x_{i+1}`$, $`1in1`$. Denote the zeros of $`D_\lambda P`$ by $`x_i(\lambda )`$, $`1in`$, those of $`P^{}`$ by $`w_j`$, $`1jn1`$, and those of $`D_\lambda P^{}`$ by $`w_j(\lambda )`$, $`1jn1`$. If these are labeled so that $`x_i(0)=x_i`$, $`1in`$, and $`w_j(0)=w_j`$, $`1jn1`$, then one can show that $`x_i(\lambda )<w_i(\lambda )<x_{i+1}(\lambda )`$, $`1in1`$, and that for any $`\lambda >0`$ and $`1mn1`$ one has

$$\underset{i=1}{\overset{m}{}}x_i^{}(\lambda )=\lambda \underset{j=1}{\overset{n1}{}}\frac{P^{\prime \prime }(w_j(\lambda )+\lambda )}{D_\lambda P^{\prime \prime }(w_j(\lambda ))}\left(\underset{i=1}{\overset{m}{}}\frac{1}{x_i(\lambda )w_j(\lambda )}\right)<0.$$

Thus $`\lambda _{i=1}^mx_i(\lambda )`$, $`1mn1`$, are decreasing functions on $`]0,\mathrm{}[`$, which combined with Theorem 7 proves Theorem 5 for strictly hyperbolic polynomials (the case $`\lambda <0`$ is similar). The same density arguments as those used for Theorem 4 show that Theorem 5 is true for all $`P_n`$.

The second part of the statement in Theorem 6 follows from the first part by noticing that $`P(x)_n`$ if and only if $`(1)^nP(x)_n`$. Keeping the same notations as above, one can show that if $`P_n`$ is strictly hyperbolic then

$$x_n^{\prime \prime }(\lambda )=2(x_n^{}(\lambda )+1)^2\underset{i=1}{\overset{n1}{}}\frac{w_ix_i(\lambda )\lambda }{(x_n(\lambda )+\lambda w_i)(x_n(\lambda )x_i(\lambda ))}>0$$

for any $`\lambda `$, which proves Theorem 6 in the generic case. If $`P`$ has multiple zeros then one can use strictly hyperbolic polynomials of the form $`\left(1\frac{1}{s}\frac{d}{dx}\right)^{n1}P`$, $`s_+`$, in order to approximate the function $`\lambda \mathrm{max}Z(D_\lambda P)`$ uniformly on compact intervals by convex $`C^2`$-functions. This proves Theorem 6.

## 2. Remarks and open questions

The manifold $`𝒞_n`$ of monic complex polynomials of degree $`n`$ is a natural context for discussing possible extensions of the above results to the complex case. By analogy with the hyperbolic case we may view $`(𝒞_n,)`$ as a partially ordered set, where the ordering relation $``$ is now induced by the spectral order on $`n`$-tuples of vectors in $`^2`$ (cf. Theorem 1). This means that the zero sets of polynomials in $`𝒞_n`$ are viewed as subsets of $`^2`$ and that if $`P,Q𝒞_n`$ then $`PQ`$ if and only if $`Z(P)Z(Q)`$.

The following example shows that if the partial ordering $``$ on $`𝒞_n`$ is defined as above then one cannot expect a complex analogue of Theorem 3.

###### Proposition 3.

Let $`P(z)=z^n1`$ and $`\lambda `$. If $`n3`$ and $`|\lambda |`$ is small enough then $`D_\lambda P`$ and $`P`$ are incomparable with respect to the partial ordering $``$ on $`𝒞_n`$.

We also note that the results of the previous section (for the hyperbolic case) are valid only for real values of the parameter $`\lambda `$:

###### Proposition 4.

Let $`P_n`$ be a strictly hyperbolic polynomial and let $`\lambda `$. If $`n2`$ and $`|\lambda |`$ is sufficiently small then $`D_\lambda P`$ and $`P`$ are incomparable with respect to the partial ordering $``$ on $`𝒞_n`$.

These examples suggest that complex generalizations of Theorem 3 – if any – should involve only classical majorization and real values of the parameter $`\lambda `$. Based on extensive numeric calculations, we make the following

###### Conjecture 1.

If $`P𝒞_n`$ then $`\mathrm{}Z(P)\mathrm{}Z(D_\lambda P)`$ for any $`\lambda `$.

We end with a few questions related to the semigroups $`𝒮`$ and $`\stackrel{~}{𝒮}`$. Let $`_n`$ denote the set of all monic real polynomials of degree $`n`$. There is reason to believe that any two $`𝒮`$-orbits in $`_n`$ have a non-empty intersection:

###### Conjecture 2.

If $`P_1,P_2_n`$ then there exist differential operators $`\mathrm{\Lambda }_1,\mathrm{\Lambda }_2𝒮`$ such that $`\mathrm{\Lambda }_1P_1=\mathrm{\Lambda }_2P_2`$.

If true, Conjecture 2 would imply in particular that $`𝒮P_n\mathrm{}`$ for any $`P_n`$, which would answer in the affirmative a question of I. Krasikov.

Let finally $`P_n`$ and set $`P_{}=\{Q_nPQ\}`$. One can easily check that if $`P`$ is strictly hyperbolic and $`n3`$ then $`\stackrel{~}{𝒮}PP_{}`$. It would be interesting to know whether there exists a (semi)group $`𝒟`$ of differential operators (not necessarily with constant coefficients) such that $`𝒟\stackrel{~}{𝒮}`$ and $`P_{}=𝒟P`$ for any $`P_n`$. This would give a completely new way of describing classical majorization.

Acknowledgements. The results announced in this paper and their proofs were all obtained by the first author (JB). The second author (BS) only made some computer tests in connection with Remark 1 and Theorem 3. Both authors are grateful to Harold Shapiro and Ilia Krasikov for stimulating discussions.


text:
# Optimal representations of quantum states by gaussians in phase space

## 1 Introduction

In many situations there is a need to represent a given quantum state or a probability distribution of states in terms of a simple basis which can be easily handled in further processing, e.g., for integration in matrix elements. Gaussians form such a basis although it is overcomplete. Gaussians have the additional advantage that they can adopt the form of coherent states in phase space, thereby supplying a direct link to classical mechanics. Hence, the formulation in terms of gaussians is also very useful for a mixed quantum-classical or quasiclassical description of a system. In particular the evolution of systems with many degrees of freedom can often be described within classical molecular dynamics only. Here, one is interested in obtaining optimal initial conditions by translating the initial quantum state of the system into a suitable classical phase space distribution which can be classically propagated . While all quantum information on an amplitude level is lost in this approximation, we will introduce a representation which can preserve basic quantum interferences by construction even in the case of the purely classical propagation. Finally, it is of general and basic interest to reconstruct amplitudes or the density matrix, i.e., coherent information, from phase space distributions. This has been done recently in terms of the moments in time of the Wigner distribution . Although there is no principal need to do so we will concentrate on the Wigner distribution due to the attention it has received recently in different areas. In several experiments the Wigner distribution of microscopic systems has been measured . Beyond the quantum context, in the time-frequency phase space of signal processing, analysis of the Wigner distribution  is an important tool , as for optics , speech analysis  or the monitoring of machinery conditions .

The decomposition of a given state $`|\mathrm{\Psi }`$ into coherent states can be regarded as a particular form of phase space representation. Labeling positions by $`q`$ or $`Q`$ and momenta by $`p`$ or $`P`$, we translate the ground state of the harmonic oscillator by $`X=(P,Q)`$ to obtain the coherent state in position representation, i.e. its wavefunction is

$`q|X=\left({\displaystyle \frac{m\omega }{\pi \mathrm{}}}\right)^{1/4}\mathrm{exp}\left({\displaystyle \frac{\omega }{2\mathrm{}}}(qQ)^2+i{\displaystyle \frac{P}{\mathrm{}}}\left(q{\displaystyle \frac{Q}{2}}\right)\right).`$ (1)
For simplicity, we chose unit frequency ($`\omega =1`$) and mass ($`m=1`$) for the harmonic oscillator without loss of generality. (Often the complex variable $`Z=(Q+iP)/\sqrt{\mathrm{}}`$ is used and the constant phase $`PQ/2\mathrm{}`$ is omitted). Even though the coherent state basis is overcomplete, the exact representation,

$`|\mathrm{\Psi }={\displaystyle \frac{1}{\pi }}{\displaystyle 𝑑X|XX|\mathrm{\Psi }},`$ (2)
is unique . However, such a decomposition into coherent states is monstrously inefficient, as shown by the overlap of coherent states,

$`X|X^{}`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{(QQ^{})^2}{4\mathrm{}}}{\displaystyle \frac{(PP^{})^2}{4\mathrm{}}}+{\displaystyle \frac{i}{2\mathrm{}}}(QP^{}Q^{}P)\right)`$ (3)
$``$ $`\mathrm{exp}\left({\displaystyle \frac{(XX^{})^2}{4\mathrm{}}}+{\displaystyle \frac{i}{2\mathrm{}}}XX^{}\right),`$
so that

$`|X|X^{}|^2=\mathrm{exp}\left({\displaystyle \frac{(XX^{})^2}{2\mathrm{}}}\right).`$ (4)
We have used the wedge product,

$`XX^{}=PQ^{}QP^{}=(𝐉X)X^{},`$ (5)
to shorten the notation. (With the second equation we also define the symplectic matrix $`𝐉`$). It follows that one can approximate any given state by a finite expansion

$`|\mathrm{\Psi }{\displaystyle \underset{j}{}}a_j|X_j`$ (6)
with arbitrary accuracy. One possibility would be to place the coherent states on a grid and then determine the complex coefficients $`a_j`$ by minimizing the error, defined in some convenient way. Of course, it is much better to use some features of $`|\mathrm{\Psi }`$ as a preliminary guide to where the $`|X_j`$ should be located.

An obvious necessary condition for a good fit is that the coherent state intensity $`|X|\mathrm{\Psi }|^2`$ in $`(\text{4})`$ be well approximated by

$`|X|\mathrm{\Psi }|^2`$ $`=`$ $`|{\displaystyle \underset{j}{}}a_jX|X_j|^2`$ (7)
$`=`$ $`{\displaystyle \underset{j}{}}|a_j|^2|X|X_j|^2+{\displaystyle \underset{jk}{}}a_ja_k^{}X|X_jX|X_k^{}`$
for each phase space point $`X`$. But note that, if the diagonal terms in this expression are only appreciable for $`X_j`$ very close to $`X`$, a nondiagonal term decays as a product of gaussians centered on a pair of separate phase space vectors. Such a term is negligible for all $`X`$ unless $`X_j`$ is very close to $`X_k`$. Moreover, such nondiagonal terms are not positive definite, which justifies our performing the preliminary fit

$`\rho _\mathrm{\Psi }(X)={\displaystyle \frac{1}{\pi }}|X|\mathrm{\Psi }|^2{\displaystyle \underset{j}{}}|a_j|^2{\displaystyle \frac{1}{\pi }}|X|X_j|^2={\displaystyle \underset{j}{}}|a_j|^2G_j(X),`$ (8)
where

$`G_j(X)=\mathrm{exp}\left({\displaystyle \frac{(XX_j)^2}{2\mathrm{}}}\right).`$ (9)
The factor $`1/\pi `$ has here been introduced so that the left hand side becomes the Husimi function , i.e., the Husimi representation of the state $`|\mathrm{\Psi }`$, also known as the $`Q`$-function in quantum optics . Here, $`\rho _\mathrm{\Psi }`$ is then decomposed into gaussians phase space packets (gpps) which are Husimi representations of the coherent states, specified by $`(\text{4})`$.

The Husimi function in its own right uniquely defines a quantum state, though part of the quantum phase information is locked up into its delicate analytic properties. These are hard to translate into other representations, but we propose to sidestep this problem. Our method is to make a preliminary fit of the Husimi function $`(\text{8})`$, which is real and positive definite, and therefore much simpler than fitting the complex coherent states directly in $`(\text{6})`$. However we can now identify the coefficients of the Husimi fitting as being proportional to the square modulus of the coefficients of the superposition of coherent states $`(\text{6})`$. The phases of the coefficients are then determined by relating the Husimi function to the coarse-graining of the Wigner function, i.e. we bring in yet a third phase space representation of the state $`|\mathrm{\Psi }`$. It should be pointed out that, once the state $`|\mathrm{\Psi }`$ has been fitted by a finite number of coherent states, it can be viewed in any representation. The special attention given here to the Wigner function is justified by the way that it highlights quantum interferences. In fact such interferences are sometimes mimicked by extra gaussians in the context of numerical work with classical densities . Thus, in the next section we present the Wigner function, its relation to the Husimi function and its specific form for an arbitrary superposition of gaussians. In addition, we discuss and establish the differences between the representation of quantum states such as the wavefunction, the Wigner distribution and the Husimi distribution using these gaussians. The Wigner distribution sets the stage for the delicate procedure of fitting the phases, i.e., to determine the appropriate interference fringes between the gaussians as detailed in section 3. In this section we propose a method for determining the phase of the expansion coefficients by Fourier analysing the Wigner function. Actually, this resulting chord function or characteristic function (in quantum optics ) can be derived directly from the initial wave function, without obtaining first the Wigner function. The Husimi function, on the other hand, supplies a robust overall picture of the quantum states which can be used to place the gaussians. This will be discussed in section 4 using explicit examples for illustration. Emphasis will be put on general symmetry criteria for placing the gaussians and numerical demonstrations based on Monte Carlo fitting of sampled wavefunctions. Finally in section 5, we semiclassically interpret the results of section 4, namely the increase of gaussians needed to represent quantum states $`\psi _n`$ of a system with increasing quantum number $`n`$. A summary in section 6 concludes the paper.

For simplicity, all our formulae are presented for systems with a single degree of freedom, but they are easily generalized to high dimensional systems.

## 2 Wigner functions and their relations

The state $`|\mathrm{\Psi }`$ is uniquely described by the real function of phase space points $`x=(p,q)`$,

$`W_\mathrm{\Psi }(x)={\displaystyle \frac{dq^{}}{2\pi \mathrm{}}q+\frac{q^{}}{2}|\mathrm{\Psi }\mathrm{\Psi }|q\frac{q^{}}{2}\mathrm{exp}\left(\frac{i}{\mathrm{}}pq^{}\right)}`$ (10)
known as the Wigner function . Inversion of this symmetrized Fourier transform leads back to the wavefunction, while its intensity results from the projection

$`|q|\mathrm{\Psi }|^2={\displaystyle 𝑑pW_\mathrm{\Psi }(x)}.`$ (11)
Though the Wigner function can be negative for some points $`x`$, its projection is always positive. Actually, this is just a particular case of the general relation,

$`|\mathrm{\Phi }|\mathrm{\Psi }|^2=2\pi \mathrm{}{\displaystyle 𝑑xW_\mathrm{\Phi }(x)W_\mathrm{\Psi }(x)},`$ (12)
for arbitrary states $`|\mathrm{\Phi }`$ and $`|\mathrm{\Psi }`$. Thus, introducing the Wigner function for the coherent state $`|X`$,

$`W_X(x)={\displaystyle \frac{1}{\pi \mathrm{}}}\mathrm{exp}\left({\displaystyle \frac{(xX)^2}{\mathrm{}}}\right),`$ (13)
we obtain the Husimi representation of $`|\mathrm{\Psi }`$ as

$`{\displaystyle \frac{1}{\pi }}|X|\mathrm{\Psi }|^2`$ $`=`$ $`2\mathrm{}{\displaystyle 𝑑xW_X(x)W_\mathrm{\Psi }(x)}`$ (14)
$`=`$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle 𝑑x\mathrm{exp}\left(\frac{(xX)^2}{\mathrm{}}\right)W_\mathrm{\Psi }(x)}.`$
This shows that the positive definite Husimi function results from the gaussian smoothing of the Wigner function. Inserting the coherent state expansion $`|\mathrm{\Psi }`$ given by $`(\text{6})`$ into the Wigner transform (10), we obtain the well known form of the corresponding Wigner function,

$`W_\mathrm{\Psi }(x)={\displaystyle \underset{j}{}}|a_j|^2W_{X_j}+2Re\left({\displaystyle \underset{jk}{}}a_ja_k^{}W_{jk}(x)\right),`$ (15)
where we define the crossed Wigner function (or Moyal bracket )

$`W_{jk}(x)`$ $``$ $`{\displaystyle \frac{dq^{}}{2\pi \mathrm{}}q+\frac{q^{}}{2}|X_jX_k|q\frac{q^{}}{2}\mathrm{exp}\left(\frac{i}{\mathrm{}}pq^{}\right)}`$
$`=`$ $`{\displaystyle \frac{1}{\sqrt{2}\pi \mathrm{}}}\mathrm{exp}\left({\displaystyle \frac{(x\overline{X}_{jk})^2}{\mathrm{}}}\right)\mathrm{exp}\left({\displaystyle \frac{i}{\mathrm{}}}\left(x\delta X_{jk}+{\displaystyle \frac{1}{2}}X_kX_j\right)\right).`$
Here we use the abbreviated relations $`\overline{X}_{jk}=(X_k+X_j)/2`$ and $`\delta X_{jk}=X_jX_k`$, as well as the skew product $`(\text{5})`$. Comparing $`(\text{15})`$ with $`(\text{7})`$ we find the same gaussian form for the diagonal terms of the Wigner function as for the Husimi function, but the offdiagonal terms are no longer negligible. For the Wigner function, any pair of coherent states determines a third gaussian centered halfway between their two centres, of comparable amplitude to the diagonal contributions of the individual gaussians. This interference term is only cancelled by the smoothing employed to obtain the Husimi function, because of the phase oscillations. Their wavelengths in phase space are $`2\pi \mathrm{}/|\delta X|`$, so that the Husimi function is increasingly damped with growing separation of these gpps. This is an excellent illustration of the way that the Husimi function highlights the classical structure in a quantum state, while it hides quantum phase information. This question will be further developed in section 4. Of course the phase information is still contained in the small offdiagonal terms in $`(\text{7})`$, because the Husimi function is fully quantum mechanical, but this information is certainly more accessible in the Wigner representation.

We can now understand the differences between using the position representation, i.e. the wave function, the Wigner representation or the Husimi representation, to fit a given state by a set of coherent states. The position representation has the advantage that there are no interference terms between pairs of coherent states. However, if the state is highly oscillatory, such as a highly excited eigenstate of a Hamiltonian, this needs to be fitted by coherent states with large momenta $`P`$, which also have narrow oscillations. Clearly the fitting procedure becomes very unstable with respect to small errors in the initial function if large momenta are necessary. In contrast, the diagonal terms in the expansion of the Wigner function are smooth gaussians distributed over phase space, but again we obtain increasingly narrow oscillations for the interferences around the midpoint of each pair. Of course, the oscillations of the wave function and the Wigner function are not unrelated, since the latter have the same wavelength as the wave intensity in the case of a pair of gaussians located at ($`\pm P,Q`$).

Finally, the Husimi function washes out oscillations, the more effectively the tighter they happen to be. We shall review in section 4 how the smooth phase space distribution singles out the basic classical structure within the quantum state, which can be stably fitted by real gpps.

The cost is that only the location of the gaussians in phase space and the square modulus of their coefficients are then determined. In the next section we present a method for retrieving the phases of the coefficients by taking the Fourier transform of the Wigner function.

## 3 Fitting the phases

Both, the direction and the wavelength of the fringes that modulate the interference peak halfway between two gpps at $`X_j`$ and $`X_k`$ in the Wigner representation, are entirely determined by $`\delta X_{jk}`$. However, these fringes are translated by changing the phases of the coefficients $`a_j`$ and $`a_k`$, which we wish to determine. Clearly, the best way to recover these phases is to take the Fourier transform of the Wigner function:

$`\stackrel{~}{W}_\mathrm{\Psi }(\xi )`$ $`=`$ $`{\displaystyle \frac{dpdq}{2\pi \mathrm{}}W_\mathrm{\Psi }(p,q)\mathrm{exp}\left(\frac{i}{\mathrm{}}(p\xi _qq\xi _p)\right)}`$ (17)
$`=`$ $`{\displaystyle \frac{dx}{2\pi \mathrm{}}W_\mathrm{\Psi }(x)\mathrm{exp}\left(\frac{i}{\mathrm{}}(x\xi )\right)}.`$
It is also possible to express this directly from the wave function

$`\stackrel{~}{W}_\mathrm{\Psi }(\xi )={\displaystyle 𝑑qq+\frac{\xi _q}{2}|\mathrm{\Psi }\mathrm{\Psi }|q\frac{\xi _q}{2}\mathrm{exp}\left(\frac{i}{\mathrm{}}\xi _pq\right)},`$ (18)
which resembles the Wigner transform $`(\text{10})`$. This function is known as the Woodward ambiguity function in communications theory  or the characteristic function in quantum optics . Usually the sign of $`\xi _p`$ is reversed in the definition of $`\stackrel{~}{W}`$, but here we follow , where it is simply called the chord function. A justification for this term appears immediately if we consider the case of a superposition of gaussians that is our concern :

$`\stackrel{~}{W}_\mathrm{\Psi }(\xi )={\displaystyle \underset{j}{}}|a_j|^2\stackrel{~}{W}_{jj}(\xi )+{\displaystyle \underset{jk}{}}a_ja_k^{}\left(\stackrel{~}{W}_{jk}(\xi )+\stackrel{~}{W}_{kj}(\xi )\right),`$ (19)
where

$`\stackrel{~}{W}_{jj}(\xi )={\displaystyle \frac{1}{2\pi \mathrm{}}}\mathrm{exp}\left({\displaystyle \frac{\xi ^2}{\mathrm{}}}+{\displaystyle \frac{i}{\mathrm{}}}X_j\xi \right)`$ (20)
and

$`\stackrel{~}{W}_{jk}(\xi )`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{2}\pi \mathrm{}}}\mathrm{exp}\left({\displaystyle \frac{(\xi \delta X_{jk})^2}{4\mathrm{}}}\right)`$ (21)
$`\times `$ $`\mathrm{exp}\left({\displaystyle \frac{i}{\mathrm{}}}\left((\xi \delta X_{jk})\overline{X}_{jk}+{\displaystyle \frac{X_kX_j}{2}}\right)\right).`$
Thus all the diagonal terms collapse onto gaussians centered on $`\xi =0`$, whereas each interference term is a gaussian centered on the chord joining $`X_j`$ and $`X_k`$. (Both directions $`\pm \delta X_{jk}`$ are present, leading to symmetric contributions in chord space. Indeed, we must have $`\stackrel{~}{W}(\xi )=\stackrel{~}{W}^{}(\xi )`$). Let us now suppose that a given chord $`\delta X_{jk}`$ is sufficiently far from all the other chords, so that at $`\xi =\delta X_{jk}`$ the chord function is dominated by the single nondiagonal contribution $`\stackrel{~}{W}(\delta X_{jk})`$, i.e.

$`\stackrel{~}{W}_\mathrm{\Psi }(\delta X_{jk}){\displaystyle \underset{l}{}}|a_l|^2\stackrel{~}{W}_{ll}(\delta X_{jk})`$ $``$ $`a_ja_k^{}\stackrel{~}{W}_{jk}(\delta X_{jk})`$ (22)
$`=`$ $`{\displaystyle \frac{a_ja_j^{}}{2\sqrt{2}\pi \mathrm{}}}\mathrm{exp}\left({\displaystyle \frac{i}{\mathrm{}}}{\displaystyle \frac{X_kX_j}{2}}\right).`$
Then, defining

$`a_j=|a_j|\mathrm{exp}(i\theta _j)`$ (23)
and recalling that we have already fitted $`|a_j|`$ from the Husimi function, we obtain $`\theta _{jk}=\theta _j\theta _k`$ from

$`\theta _{jk}{\displaystyle \frac{1}{2\mathrm{}}}X_jX_k=phase\left(\stackrel{~}{W}_\mathrm{\Psi }(\delta X_{jk}){\displaystyle \underset{l}{}}|a_l|^2\stackrel{~}{W}_{ll}(\delta X_{jk})\right).`$ (24)
Of course, the modulus of both sides of $`(\text{22})`$ should be approximately equal, which provides a check on the previous Husimi fitting. The diagonal sum that is substracted from $`\stackrel{~}{W}_\mathrm{\Psi }(\delta X_{jk})`$ in $`(\text{22})`$ decays exponentially with $`\delta X_{jk}^2`$, so it will only affect the phases of the smaller chords. It might appear to be only consistent with our previous approximation of the Husimi function to neglect the diagonal sum in (22), but this would perturb the phases of small chords.

## 4 Placing the gaussians: examples

So far, nothing has been said about how to place the gaussians in phase space that are meant to approximate a given Husimi function. An obvious criterion is that any relative maximum of this smooth distribution should also receive a coherent state. Further knowledge of the state should also be used. For instance, if it is known to be the eigenstate of a given Hamiltonian $`\widehat{H}`$ with energy $`E_n`$, then the semiclassical considerations in the following section allow the restriction of the gaussian centers to $`|H(X)E_n|\mathrm{}`$. Of course, the fitting procedure is only useful for cases in which $`\widehat{H}`$ is unknown, or else $`|\mathrm{\Psi }`$ may itself be a linear superposition of states, such as result from a two slit experiment. Here we immediately feel the advantage that the decomposition $`(\text{6})`$ is defined in terms of the state itself, rather than the density operator $`|\mathrm{\Psi }\mathrm{\Psi }|`$. Thus, the linearity of $`|\mathrm{\Psi }=a_1|\mathrm{\Psi }_1+a_2|\mathrm{\Psi }_2`$ allows us to fit $`|\mathrm{\Psi }_1`$ and $`|\mathrm{\Psi }_2`$ independently and then to superpose the fitted state as $`|\psi =a_1|\psi _1+a_2|\psi _2`$. A further simplication results from possible symmetries. For instance, for the eigenstate of an even potential, one should place the gaussians symmetrically at $`(P_j,\pm Q_j)`$ with $`a_j=\pm a_{j+}`$, so as to guarantee even or odd states. Again, this is a very particular case, but time reversal invariance is much more common. If the state $`|\mathrm{\Psi }`$ resulted from initial real wave functions and the evolution proceeded through interactions that preserve time reversal symmetries (in usual practice, if there are no magnetic fields) then the final state can also be real. In this case, one must choose symmetric pairs of gaussians at $`(\pm P_j,Q_j)`$ with equal real coefficients, $`a_j=a_{j+}`$, so as to ensure that $`|\mathrm{\Psi }`$ is also real. This halves the number of independent coherent states to be fitted and reduces the choice of phase to either zero, or $`\pi `$. Any further knowledge should generally be used to reduce the randomness of the positions of the coherent states to be fitted to $`|\mathrm{\Psi }`$. The semiclassical considerations in the following section may be a further guideline. When this knowledge is exhausted, the best course is to optimize random guesses by Monte Carlo fitting of arbitrary gaussians, within the given constraints.

Specifically, in order to measure the quality of the fitting, we define the relative error

$`\sigma (x,A,N)={\displaystyle \underset{l=1}{\overset{M}{}}}\left(\rho _\psi (x_l){\displaystyle \underset{k=1}{\overset{N}{}}}A_kG_k(x_l)\right)^2/{\displaystyle \underset{l=1}{\overset{M}{}}}\rho _\psi ^2(x_l)`$ (25)
as the mean square deviation between the known Husimi function $`\rho _\psi (x)`$ of a quantum state $`|\psi `$ and its fitting to the superposition of gpps defined in (9). This error $`\sigma `$ has to be minimized for a set of $`M`$ grid points $`x_l=(p_l,q_l)`$. To fit this superposition of gpps as closely as possible to $`\rho _\psi (x)`$, we proceed as follows: (i) the centers $`x_k=(p_k,q_k)_{k=1,..,N}`$ of the coherent states are obtained by Monte-Carlo sampling the region where the distribution is significant; that is the region where the distribution exceeds a certain threshold $`\delta `$; (ii) the centers must not be too close, this means a minimum distance $`\mathrm{\Delta }`$ between centers is fixed a priori, thereby avoiding overlaps and reducing the total number $`N`$ of gpps used for the representation; (iii) starting with $`N=1`$ we increase N by one and repeat the process until the desired accuracy is reached. Since only the coefficients $`A_k`$ are fitted, minimizing the error $`\sigma `$ reduces to a problem of linear optimization which is equivalent to solving the real matrix equation

$`GA=\rho _\psi `$ (26)
where the phase space points $`x_l=(p_l,q_l)_{l=1,..,N}`$ that define the matrix elements $`G_{lk}=G_k(x_l)=\mathrm{exp}\left(\frac{(x_lx_k)^2}{2\mathrm{}}\right)`$ are randomly selected. Here, $`A=(A_1,A_2,\mathrm{},A_N)`$ contains the coefficients to be determined and $`\rho _\psi =(\rho _1,\rho _2,\mathrm{},\rho _N)`$ those of the known Husimi function. In general, the matrix $`G`$ is expected to be sparse and a variety of special algorithms can be used to solve (26) efficiently . In principle one may also consider to optimize this fitting with respect to the widths of the gpps; but since the Husimi function is always smooth and positive definite, it is convenient to keep them frozen. In what follows, we have tested our method on the fitting of the eigenstates of the harmonic oscillator (HO) and that of the Morse oscillator (MO) supporting 18 bound states . We have not used the knowledge that these states have time reversal symmetry in the way mentioned earlier, so we should obtain results of comparable quality even for states with complex wave functions.

Figure 1 shows, on the top, the known Husimi function for the 3rd excited state of HO (left) and for the 5th excited state of MO (right), with parameters $`\omega =m=\mathrm{}=1`$. On the bottom the corresponding fitted distributions are shown. We found that $`N=7`$ is sufficient to reproduce quite well the Husimi function of HO whereas 14 coherent states (gpps) are needed in the MO example, to achieve a global relative error of $`\sigma 0.01`$. The respective values for the threshold $`\delta `$ and the minimum distance $`\mathrm{\Delta }`$ are given in the caption of figure 1. This preliminary fit on the Husimi function provides us with the total number $`N`$ of gpps, their location in phase space $`(p_k,q_k)`$ and the expansion coefficients $`A_k`$. These coefficients are related to those of the coherent states $`a_k`$ and their phases $`\theta _k`$ by

$`a_k=\sqrt{2\pi \mathrm{}A_k}\mathrm{exp}(i\theta _k).`$ (27)
Up to this point, the phases $`\theta _k`$ are still unknown. We only need to determine the phases differences $`\theta _{jk}`$ in (15) which can be easily extracted from (24). This says that, for $`N`$ gpps, $`N(N1)/2`$ chords ($`\theta _{jk}`$) are required to reproduce the Wigner function. These $`N(N1)/2`$ pairs of gaussians (chords) are symmetrically spaced about $`\xi =0`$, so that, being more dense, they overlap more than the original set of $`N`$ gaussians in phase space. Figure 2 shows on the top the Wigner function of the 3rd excited state of the HO together with its corresponding fit obtained with $`7`$ gpps. The present fitting reproduces the Wigner function everywhere in phase space, even for negative regions which are signatures of quantum interference. This is clearly illustrated in figure 3 where the deviations of the Monte Carlo sampled Wigner function, less than $`10^4`$ everywhere in phase space, are depicted for both, the HO (left, $`n=3`$) and the MO (right, $`n=5`$). The global relative error $`\sigma `$ is subsequently found not to exceed $`0.01`$. Moreover, we have computed the number $`N`$ of gpps needed to represent the Wigner functions of the eigenstates of our test systems, as can be seen in figure 4. One sees that $`N`$ grows proportionally to $`\sqrt{n}`$ as the vibrational quantum number $`n`$ increases. In the following section, we provide a semiclassical justification of this result by means of the WKB quantization.

## 5 Semiclassical approach

A semiclassical state $`|\mathrm{\Psi }_{SC}`$ is supported by a curve in phase space, in the case of one degree of freedom, or in general by a surface with half the phase space dimension . If the curve or surface is closed, it must be Bohr-quantized and it will be symmetric about $`p=0`$ in the case of time reversal symmetry. In the case of the position representation, for each branch of the curve (surface), $`p_j(q)`$, one defines the action

$`𝒮_j(q)={\displaystyle _{q_0}^q}p_j𝑑q,`$ (28)
leading to the generalized WKB wavefunction,

$`q|\mathrm{\Psi }_{SC}={\displaystyle \underset{j}{}}a_j(q)\mathrm{exp}\left({\displaystyle \frac{i}{\mathrm{}}}𝒮_j(q)\right).`$ (29)
The amplitudes $`a_j(q)`$ can also be expressed in terms of the actions, but the important point is that they are purely classical so that the only $`\mathrm{}`$-dependence occurs in the exponential. Thus locally, for any small range of positions, the semiclassical wave function reduces to a superposition of plane waves characterized by the wave-vectors $`p_j(q)/\mathrm{}`$. The semiclassical Wigner function is also defined in terms of an action $`𝒮_j(x)`$ with respect to the classical curve (surface). But instead of the area between the curve and the $`q`$-axis, we are now concerned with the area sandwiched between the curve and one of its chords. The latter, $`\xi _k(x)`$, is selected by the property that it is centred on the point $`x`$, as shown in figure 5. Thus, for one degree of freedom ,

$`W_\mathrm{\Psi }(x)_{SC}={\displaystyle \underset{k}{}}A_k(x)\mathrm{cos}\left({\displaystyle \frac{𝒮_k(x)}{\mathrm{}}}{\displaystyle \frac{\pi }{4}}\right),`$ (30)
with straightforward generalizations . The important point is that

$`{\displaystyle \frac{𝒮_k}{x}}=𝐉\xi _𝐤,`$ (31)
where we use the symplectic matrix $`𝐉`$ defined in (5). Thus the semiclassical state is again represented by a superposition of waves,

$`W_\mathrm{\Psi }(x)_{SC}{\displaystyle \underset{k}{}}A_k(x)\mathrm{cos}\left({\displaystyle \frac{\delta x\xi _k(x)}{\mathrm{}}}{\displaystyle \frac{\pi }{4}}\right),`$ (32)
but now these are phase space waves. Comparison of this expression with (22) immediately reveals that these waves have exactly the same wave vector as a pair of gaussians at the tip of each chord $`\xi _k(x)`$. It follows that an important feature of a semiclassical state is automatically reproduced by fitting it with gaussians precisely placed on the corresponding classical curve (surface). Notice that, in fact, this is a major a priori obstacle to performing the fitting, because the oscillations near the midpoint of a long chord are very fine for $`\mathrm{}<<𝒮_k`$. One would need very narrow gaussians indeed to fit these phase space waves directly, rather than having them arise naturally as interferences. All that is left to determine the phase: $`𝒮_k/\mathrm{}\pi /4`$. How can we single out the region in which to place the gpps, if we have no a priori knowledge of the classical structures correponding to a quantum state? The obvious course is to smooth the Wigner function with a Gaussian window. In other words, we initially fit the Husimi function. Recalling that the gpps have linear width of order $`\sqrt{\mathrm{}}`$, while the wavelength of the Wigner oscillations are $`\mathrm{}/|\xi _k|`$, it follows that this gaussian window erases effectively a semiclassical Wigner function, except in the limit of very small chords. This is a simple explanation of the well known fact that only in the neighbourhood of the classical curve (surface) itself is the semiclassical Husimi function appreciable. In the simple case of one degree of freedom, the Husimi function is concentrated near the energy shell. Thus, by fitting gpps to a Husimi function, one automatically samples its relevant classical mainifold, if it happens to have a (possibly unknown) semiclassical structure. A basic assumption in the above argument is that the gpps along the curve (surface) are neither too crowded, to avoid confusing superpositions, nor too sparsely spaced, lest gaps should arise in the fitting. It is thus easy to estimate the growth of the number of coherent states required to fit the eigenstate $`|n`$ with the degree of excitation. Since the linear width of the gaussian scales as $`\sqrt{\mathrm{}}`$, the phase space area of the curve grows as $`n\mathrm{}`$ and hence its length grows as $`\sqrt{n\mathrm{}}`$, the number of equispaced gaussians needed to cover the curve grows as $`\sqrt{n}`$. Of course, we have assumed here that the shape of the curve does not change with $`n`$, as for the harmonic oscillator. If the eigencurve elongates for higher excitation, the number of gpps necessary for a good fit will grow as $`n^\alpha `$, with $`1/2\alpha 1`$. We can also estimate the growth in time of the number of gpps required to fit an evolving semiclassical state with arbitrary precision. According to the theory of van Vleck , it is sufficient to evolve the curve (surface) classically and then to reconstruct the wave function. This principle can also be applied to Wigner functions, according to Berry and Balazs . If the driving Hamiltonian is chaotic, then the curve will stretch at a rate depending on the Lyapunov exponent, $`\lambda `$. Clearly, this also determines the initial rate of growth of the number of gpps needed for an adequate fit. Ultimately, when the curve covers densely all the available phase space (the energy shell of the driving Hamiltonian), a steady state will be reached where the number of gpps saturates.

## 6 Conclusions

Optimal fitting of quantum states by phase space gaussians is achieved by first fitting the Husimi function and then determining the quantum phases from the chord function, i.e. the Fourier transform of the Wigner function. It is clear from the discussion in the introduction that the density of gaussians in phase space must be finely adjusted: If ther are too few, essential features of the state will be missed, whereas an excessive number of gaussians would introduce interference terms in the Husimi function itself, which could only be accomodated by a much more complex variation of this method. The numerical examples in section 4 indicate clearly that , as well as achieving excellent overall accuracy, all essential qualitative features of the states are captured by this method, using a basis of gaussians that grows more slowly than the excitation number of the fitted states. This last result and further insight into the fitting procedure, follows from the semiclassical analysis in the previous section. This can be generalized readily to states of quantum systems with higher degrees of freedom if these are eigenstates of integrable systems (i. e. if they are supported by a Lagrangean surface in phase space; see e. g. ). Of course, more gaussians will then be needed for the fitting, but the preliminary fit of the Husimi function should still provide optimal results. The potential for this method to deal with the eigenstates of chaotic systems is even more interesting. Though the chord and centre description of the Wigner function still applies to mixtures of eigenstates over narrow energy windows, no classical theory accounts for individual eigenstates at present. Therefore, it will be extremely useful to describe these as interfering superpositions of gaussians placed near the energy shell. Of course, the size of this basis would diverge at the classical limit as $`\mathrm{}0`$, but manageable approximations should be attainable for finite excitations.

AK gratefully acknowledges the financial support by Alexander von Humboldt (AvH) Foundation/Bonn-Germany, under the grant of Research fellowship No.IV.4-KAM 1068533 STP. AMOA thanks the MPIPKS-Dresden for a Martin Gutzwiller Fellowship, during which this work was initiated, and acknowledges support from CNPq and Pronex in Brazil.

## References


text:
# Planetary migration in protoplanetary disks

## 1. Introduction

In order to study the formation of planetary systems it is necessary to study the global evolution of solid material which constitutes, together with gas, the protoplanetary discs. In the following, we introduce a time-dependent accretion disc model that shall be used in the next sections to study planets migration.

It is usually assumed that the surface density in planetesimals is proportional to that of gas (Lin & Papaloizou (1980, page 47); Murray et al. (1998)) and that the spatial structure and the time evolution of the surface density, $`\mathrm{\Sigma }_\mathrm{s}(R,t)`$, is such that $`\mathrm{\Sigma }_\mathrm{s}(R,t)\mathrm{\Sigma }(R,t)`$. Although this approximation has been and is widely used in literature, a better way to proceed is to directly calculate the evolution of solids in the disc. We know that the evolution of the surface density of gas is described by a diffusive-type equation, while that of solid particles is an advection-diffusion equation:

$$\frac{\mathrm{\Sigma }_\mathrm{s}}{t}=\frac{3}{r}\frac{}{r}\left[r^{1/2}\frac{}{r}(\nu _\mathrm{s}\mathrm{\Sigma }_\mathrm{s}r^{1/2})\right]+\frac{1}{r}\frac{}{r}\left[\frac{2r\mathrm{\Sigma }_s\overline{v}_\varphi _s}{\mathrm{\Omega }_kt_s}\right]$$
(1)
where $`\nu _\mathrm{s}=\frac{\nu }{\mathrm{Sc}}`$, the Schmidt number, Sc, is given by:

$$\mathrm{Sc}=\left(1+\mathrm{\Omega }_\mathrm{k}t_\mathrm{s}\right)\sqrt{1+\frac{\overline{𝐯}^2}{V_\mathrm{t}^2}}$$
(2)
where $`𝐯`$ is the relative velocity between a particle and the gas, $`V_\mathrm{t}`$ the turbulent velocity, $`\mathrm{\Omega }_\mathrm{k}`$ is the Keplerian angular velocity, $`t_\mathrm{s}`$ the so called stopping time.

Stepinski & Valageas (1996, 1997) developed a method that, using a series of simplifying assumptions, is able to simultaneously follow the evolution of gas and solid particles due to gas-solid-coupling, coagulation, sedimentation, and evaporation/condensation for up to $`10^7\mathrm{yr}`$.

The equation of gas evolution is solved by means of an implicit scheme and the evolution of gas is computed independently from the evolution of particles. At every time step the quantities needed for evaluating the change in the mass distribution of solids are calculated and the change itself is computed from Eq. (1) using the operator splitting method. In such a method the advective term in Eq. (1) is treated by the numerical method of characteristics, whereas an implicit scheme is applied to the diffusion term. The obtained distribution of solid material is then modifed because of the existence of the evaporation radius, and the mass distribu- tion of the vapor is calculated using the implicit scheme. Finally, the new particle size distribution is calculated be- fore proceeding to the next time step.

We suppose that a single planet moves in a planetesimal disc under the influence of the gravitational force of the Sun. The equation of motion of the planet can be written as:

$$\ddot{𝐫}=𝐅_{}+𝐑$$
(3)
(Melita & Woolfson 1996), where the term $`𝐅_{}`$ represents the force per unit mass from the Sun, while $`𝐑`$ is the dissipative force (the dynamical friction term-see Melita & Woolfson 1996).

In order to take into account dynamical friction, we need a suitable formula for a disc-like structure such as the protoplanetary disc.

We assume that the matter-distribution is disc-shaped, then we have that:

$$𝐑=k_{}v_1𝐞_{}k_{}v_1𝐞_{}$$
(4)
where $`𝐞_{}`$ and $`𝐞_{}`$ are two versors parallel and perpendicular to the disc plane and $`k_{}`$ and $`k_{}`$ are given in Del Popolo & Ekşi (2002).

## 2. Results

Our model starts with a fully formed gaseous giant planet of $`1M_\mathrm{J}`$ at $`5.2`$ AU. According to several evidences showing that the disc lifetimes range from $`10^5`$ yr to $`10^7`$ yr (Strom et al. 1993; Ruden & Pollack 1991), we assume that the disc has a nominal effective lifetime of $`10^6`$ years (Zuckerman et al. 1995). We integrated the model introduced in the previous section for several values of disc masses: $`M_\mathrm{D}=0.1`$, 0.01, 0.001, 0.0001 $`M_{}`$, and several values of $`\alpha `$.

The results of the disc model are plotted in Fig.1. Fig. 1a shows the evolution of the density of the gas component of the disc in the case $`M=0.1M_{}`$ and $`\alpha =0.1`$ (see figure). As shown in Fig. 1b, in agreement with Stepinski & Valageas (1997) and Kornet et al. (2001), the most important result of the low-mass models calculation (discs with $`M<0.1M_{}`$, and angular momentum $`10^{52}\mathrm{gcm}^2/\mathrm{s}`$) is that such models lead to the survival of solid material, as can be seen from the evolution of the particle size, or by the emergence of the converged, nonvanishing surface density distribution of solids.

In the successive figures (2a-2d), we plot the evolution of semi-major axis of the planet (see figure). Summarizing figures 2a-2d, according to the final distribution of planets distances, the present model predicts that, unless the disc mass is very small $`M0.0001M_{}`$, planets tend to move close to the central star to distances of the order of $`0.03\mathrm{AU}`$. It is also evident that it is possible to find a planet at any distance from their locations of formation and very small distances from the parent star for peculiar values of the parameters $`\alpha `$ and $`M`$.

In Fig. 3a, in order to show the predictions of the model for what concerns the distribution of planets in the inner part of the disc, we plot the fraction of planets in the orbital period range 0-20 days calculated using the model of this paper. Fig. 3b, represents the same distribution obtained with the data given in www.exoplanets.org (see also Kuchner & Lecar 2002).

According to Kuchner & Lecar 2002, disc temperature determines the orbital radii of the innermost surviving planets, similarly to our model. As streseed in DP1, DP2, the model has not the drawback of Murray et al. 1998 model, namely that of requiring a too large disc mass for migration and at the same time has the advantage of Murray et al. 1998 model of having an intrinsic natural mechanism that provides halting of migration.


text:
# WMAP and the Generalized Chaplygin Gas

## 1 Introduction

The picture of a universe mostly filled by an unknown component, Dark Energy (DE henceforth), seems overwhelmingly justified by the latest observations of CMB anisotropies (Spergel et al. 2003).

While there is a strong degeneracy of the multitude of DE candidates for the background evolution – and therefore in the SNIa data (Perlmutter et al. 1999; Riess et al. 1998) –, additional information should be contained in the linear fluctuations. However, in order to differentiate among DE models, the avenue of studying DE fluctuations is not promising as it seems. Infact, in the simplest uncoupled quintessence models (Ratra & Peebles 1988; Frieman et al. 1995; Caldwell, Dave & Steinhardt 1995) based on scalar fields with ordinary Klein-Gordon Lagrangian, quintessence fluctuations are evanescent on sub-horizon scales. This is due to the presence of a Jeans scale for scalar field perturbations, of the order of the Hubble radius. This Jeans scale is related to a unitary sound speed (in $`c=1`$ units) for scalar field perturbations. For the same reason, isocurvature perturbations caused by an offset of the quintessence component (Abramo & Finelli 2001) are not observationally harmful as other isocurvature modes.

A possibility of having a smaller Jeans scale for dark energy fluctuations is offered by theories with non standard kinetic terms for the scalar field candidate. A scalar field with a Born-Infeld action is an example:

$$S=d^4x\sqrt{g}V\sqrt{1+\frac{_\mu \phi ^\mu \phi }{M^4}},$$
(1)
where $`M`$ is a fundamental scale and $`V`$ is a potential. Recently, the tachyon field (with $`V=V(\phi )`$) motivated from string theory has received a lot of attention (Sen 2002, Gibbons 2002, Padmanabhan 2002). A Chaplygin gas is a hydrodinamical description of the scalar field with the Born-Infeld Lagrangian in Eq. (1) with constant $`V`$ (Kamenshchik, Moschella & Pasquier 2001). A Chaplygin gas is characterized by a pressure which is inversely proportional to the energy $`p1/\rho `$. A generalized Chaplygin gas (henceforth GCG) is a perfect fluid with a polytropic equation of state (Bento, Bertolami & Sen 2002):

$$p_X=\frac{A}{\rho _X^\alpha }$$
(2)
where $`\alpha `$ is a parameter between $`0`$ and $`1`$ in order to have a sound speed at most luminal for perturbations. $`A`$ is a constant with dimensions $`[M^{4(1+\alpha )}]`$. A GCG is a phenomenological extension of the Chaplygin gas, which exhausts all the possibilities for a polytropic perfect fluid dark energy candidate, whose perturbations are stable on small scales (Fabris, Goncalves & Souza 2002a, Carturan & Finelli 2002 – henceforth CF –).

The behaviour in time of a GCG interpolates between dust and a cosmological constant, with an intermediate behaviour as $`p=\alpha \rho `$. The Jeans instability of GCG perturbations is first similar to CDM fluctuations (when the GCG has a negligible pressure) and then disappears (when the GCG behaves as a cosmological constant) (see CF). Both this late suppression of GCG fluctuations and the apperance of a non zero Jeans length leave a large integrated Sachs-Wolfe (ISW) imprint on the CMB anisotropies (CF). A model in which the GCG is the only dark component is referred to as a unified dark matter/dark energy (UDM) model.

The viability of GCG as DE and UDM model has been first analyzed in the context of Supernovae (Fabris et al. 2002b; Avelino et al. 2002; Makler, de Oliveira & Waga 2003) and other complementary data (Dev, Jain & Alcaniz 2003, Alcaniz, Jain and Dev 2003). The comparison with CMB data has been the next step (CF, Bean & Dore’ 2003) and it has led to much stronger constraints (Bean & Dore’ 2003).

The imprint of a GCG is also present on the matter power spectrum, since GCG perturbations affect both CMB anisotropies and structure formation. This last issue has been considered by Sandvik et al. (2002) in the context of UDM models, claiming that a GCG is ruled out as a UDM candidate at $`99.999\%`$ (Sandvik et al. 2002). However, their analysis does not take into account baryons, while we will show they are relevant for the shape of the total matter power spectrum. Moreover, the analysis by Sandvik et al. is based on a linear treatment of perturbations until the present time, neglecting any non linear effects which may be important and unexpected because of the time dependence of the GCG Jeans length. Such nonlinear effects should be more important for the present matter power spectrum than for the CMB spectrum. For this reason, we focus on the comparison with CMB data.

In this paper we compare the predictions of GCG models with the recently released WMAP data, which provide the most precise determination of the CMB spectrum so far (Hinshaw et al. 2003), and with SNIa data (Perlmutter et al. 1999, Riess et al. 1998).

## 2 The Model

In this section we review a brief and basic description of the model (for more details see CF). The conservation equation for the GCG in a Robertson-Walker metric is solved by

$$\rho _X=[A+\frac{B}{a^{3(\alpha +1)}}]^{\frac{1}{1+\alpha }},$$
(3)
where $`a`$ is the scale factor ($`a=1`$ today) and $`B`$ is an integration constant with dimension $`[M^{4(1+\alpha )}]`$.

The equation of state is therefore

$$w_X(a)=\frac{p_X}{\rho _X}=\frac{A}{\rho _X^{1+\alpha }},$$
(4)
while the present equation of state is

$$w_X=\frac{A}{A+B}.$$

Instead of $`A,B`$ we find convenient to employ the parameters $`\mathrm{\Omega }_X,w_X`$ , where $`\mathrm{\Omega }_X`$ is the present value of the density parameter of the GCG and

$`A`$ $`=`$ $`w_X(\mathrm{\Omega }_X\rho _c)^{1+\alpha },`$
$`B`$ $`=`$ $`(1+w_X)(\mathrm{\Omega }_X\rho _c)^{1+\alpha },`$
where $`\rho _c`$ is the present critical density. Beside the GCG we introduce for generality also baryons and CDM. Although the GCG is very interesting as an unified model of dark energy and dark matter, it is conceivable that it is in fact only a DE component additional to ordinary CDM and baryons.

Perturbations of GCG are stable on small scales since the sound speed is positive ($`c_X^2=p_X/\rho _X=\alpha w_X`$). An important consequence of assuming a perfect fluid instead of a scalar field as a candidate of dark energy is the absence of intrinsic non adiabatic pressure perturbations. Infact, for a GCG pressure perturbations are locked to density perturbations

$$\delta p_X=c_X^2\delta \rho _X.$$
(5)
Instead, in presence of non adiabatic pressure perturbations, we have in general:

$$\delta p_X=c_X^2\delta \rho _X+\mathrm{\Delta }p_X.$$
(6)
For example, in a scalar field model with standard kinetic term for quintessence $`c_X^2=1`$ and $`\mathrm{\Delta }p_X=2V_\varphi \delta \varphi `$, where $`V_\varphi `$ is the first derivative of the potential and $`\delta \varphi `$ is the field fluctuation.

A perfect fluid description of a GCG offers the chance to explore the possibility of the basic relation (5) for dark energy pressure perturbations.

In Figs. 1 and 2 we compare the $`C_{\mathrm{}}`$ temperature spectrum with the WMAP data varying $`w_X`$ and $`\alpha `$, respectively (see also CF). These figures show how strongly the $`C_{\mathrm{}}`$ spectrum depends on perturbations of GCG, when this is a DE candidate ($`\mathrm{\Omega }_c`$ is fixed to $`0.27`$). Fig. 1 shows how the resulting CMB spectrum strongly differ even for models which have $`w_X0.91`$. Notice that the 2nd and 3rd peak move to the left increasing $`w_X`$, because then the expansion is less and less accelerated, while the 1st peak remains more or less stable because of the concurring strong integrated Sachs-Wolfe (ISW) effect. Fig. 2 shows the strong dependence of CMB spectrum on $`\alpha `$, which regulates the Jeans length of GCG perturbations and therefore the amount of ISW effect. The peak-to-plateau ratio decreases for $`\alpha `$ increasing, because after the dust behaviour the GCG induces a strong ISW. The dependence of CMB spectrum and LSS on the Jeans length of DE perturbations is a generic effect and it has been also found recently in K-essence models (DeDeo, Caldwell & Steinhardt 2003), but it is less severe than in GCG models.

## 3 Likelihood analysis

We compare the models to the combined power spectrum estimated by WMAP (Hinshaw et al. 2003). To derive the likelihood we adopt a version of the routine described in Verde et al. (2003), which takes into account all the relevant experimental properties (calibration, beam uncertainties, window functions, etc). Since the likelihood routine employs approximations that work only for spectra not too far from the data, we run it only for models whose $`\chi ^2`$ per d.o.f. is less than 1.5 from the number of degrees of freedom.

Our theoretical model depends on two GCG parameters, four cosmological parameters and the overall normalization $`N`$:

$$\alpha ,w_X,n_s,h,\omega _b,\omega _c,N.$$
(7)
where $`\omega _{b,c}=\mathrm{\Omega }_{b,c}h^2`$ and $`n_s`$ is the primordial fluctuation slope. As anticipated, to save computing time we found necessary to restrict the analysis to a flat space; moreover, we fixed the optical depth to $`\tau =0.17`$, the mean value found by WMAP (Spergel et al. 2003). The overall normalization has been integrated out numerically. We calculate the theoretical $`C_{\mathrm{},t}`$ spectra by a modified parallelized CMBFAST (Seljak & Zaldarriaga 1996) code that includes the full set of perturbation equations (see CF). We do not include gravitational waves and the other parameters are set as follows: $`T_{cmb}=2.726K,`$ $`Y_{He}=0.24,N_\nu =3.04`$.

We evaluated the likelihood on a non-uniformly spaced grid of roughly $`50,000`$ models (for each normalization) with the following top-hat broad priors: $`w_X(1,0.5),`$ $`\alpha (0,.5),`$ $`n_s(0.8,1.2),`$ $`\omega _b(0.005,0.04),`$ $`\omega _c(0.0,0.2)`$. For the Hubble constant we adopted the top-hat prior $`h(0.5,0.9);`$ we also employed the HST result (Freedman et al. 2001) $`h=0.72\pm 0.08`$ (Gaussian prior). One reason to adopt a grid approach than a Markov chain method (as in Verde et al. 2003) is that the independent grid evaluation can be parallelized with maximum efficiency.

## 4 Results from comparison with WMAP

In Fig. 3 we plot the likelihood functions for each parameter, marginalizing in turn over the others. The horizontal lines mark the 68% and 95% CL; the vertical lines the 68% and 95% upper bounds. We obtain the bounds

$`\alpha `$ $`<`$ $`0.2(0.05)`$
$`w_x`$ $`<`$ $`0.8(0.92)`$ (8)
at 95%(68%) CL. These constraints are the most stringent ones obtained so far using only CMB data; the limit on $`\alpha `$ is 5–10 stronger than the previous one obtained with pre-WMAP data (Bean & Dore’ 2002). For the other parameters we obtain results not very different than the standard ones. We also show the results imposing the HST prior $`h=0.72\pm 0.08`$; they are almost identical to the one without prior on $`h`$. The best fit DE model is $`n_s=1`$, $`\alpha =0`$, $`w_X=0.98`$, $`h=0.71`$, $`\omega _b=0.024`$, $`\omega _c=0.12`$ has $`\chi ^2`$/d.o.f.=979/892.

## 5 Constraints on the Unified Dark Model

Any DE model which interpolates in time between a dust epoch and an accelerating stage has the possibility of playing the role of a UDM candidate, if one renounces to the presence of CDM. Naturally, in such a case one should explain also why 30% or so of the matter collapsed into galaxies and clusters while the remaining did not. The GCG has received a lot of attention also in this light of a UDM model. While SNIa data allow for this possibility, it was claimed that LSS data completely ruled it out (Sandvik et al. 2002). Here we find that WMAP data alone (for which a linear analysis is sufficient) show that a GCG as a unified candidate is much less likely than a model which includes a CDM component. The likelihood for $`\omega _c=0`$ is infact roughly 50 times smaller than for $`\omega _c=0.12`$; formally, this rules out the UDM at slightly more that 99.99% C.L. when compared with GCG as DE. The lowest value allowed by WMAP data is 0.1 (95C.L.). However, if we calculate the $`\chi ^2`$ for the best fit among the UDM cases (corresponding to a model with $`n_s=0.98`$, $`\alpha =0`$, $`w_X=0.8`$, $`h=0.82`$, $`\omega _b=0.026`$) we find $`\chi ^2`$/d.o.f.=986/893, indicating that the UDM model is an acceptable fit to WMAP data at at 2$`\sigma `$ level.

Note that our results have been obtained varying also $`n_s`$, which has been instead fixed in the analysis of Bean & Dore’ (2002) and of Sandvik et al. (2002) (who also fixed $`h`$).

## 6 Adding the Supernovae constraints

We now proceed with a further refinement of our result, by taking into account the constraint from the SN I a Hubble diagram. This has been evaluated by several groups (Fabris et al. 2002a; Avelino et al. 2002; Makler et al. 2003), but not in a suitable form for our purposes.

The luminosity distance in a GCG model in flat space is

$$d_L(z)=(1+z)_0^z\frac{dz^{}}{H(z^{})},$$

where the Hubble function is

$$H^2(z)=H_0^2\left[\mathrm{\Omega }_m(1+z)^3+(1\mathrm{\Omega }_m)((1+w_X)(1+z)^{3(1+\alpha )}w_X)^{\frac{1}{1+\alpha }}\right],$$

and $`\mathrm{\Omega }_m=\mathrm{\Omega }_b+\mathrm{\Omega }_c`$.

We compare the luminosity distance to the SCP data of Perlmutter et al. (1999) (their fit C), to which we add the supernova at $`z1.7`$ (Benitez et al. 2002). We show in Fig. 6 the 2-dimensional likelihood function marginalized over the “nuisance parameters” (see definitions in Perlmutter et al. 1999) $`\alpha ,M^{}`$ , and over $`\mathrm{\Omega }_m`$ with Gaussian prior $`\mathrm{\Omega }_m=0.3\pm 0.1`$. Notice that the contours are almost insensitive to $`\alpha `$ since $`d_L(z)`$ at small redshifts is independent of $`\alpha `$. This shows how CMB data constrain the GCG models much more than SNIa data. In Fig. 7 we multiply the SN and WMAP likelihood functions. The final constraints turn out to be almost identical to those in (8).

## 7 Considerations on the Mass Power Spectrum

The influence of the GCG on structure formation is another important observational imprint of this candidate for DE and UDM model. GCG fluctuations accumulate on all scales when the GCG behaves as CDM. As soon as the GCG equation of state drops below zero, the GCG fluctuations are no longer Jeans unstable at all scales, oscillating and decaying. This exotic behaviour of fluctuations occurs also when the GCG simulates $`\mathrm{\Lambda }`$. As already said in the introduction, this suppression of fluctuations is in part responsible for the large ISW effect in CMB anisotropies. It is clear that a similar behaviour of fluctuations differs both from CDM and standard quintessence, and it is more important in absence of CDM, i.e. in a UDM model.

Sandvik et al. (2002) used a comparison of the predicted power spectrum in GCG models with LSS data – 2dF survey (Colless et al. 2001; see also Tegmark, Hamilton & Xu 2002) – to rule out the GCG as a UDM candidate at 99.999% CL. Their approximate analysis neglects baryons (this has also been noted recently by Beca et al. 2003) and is based on a linear treatment of perturbations until the present time. However, the nonlinear clustering of GCG could be rather nontrivial because of the time dependence of its Jeans length and deserves more study.

In Fig. 8 we compare the predicted power spectra – total (solid line), GCG (dashed line) and baryon (dotted line) – for UDM model with the APM (Padilla & Baugh 2003) and 2dF (Colless et al. 2001) data when $`\alpha `$ varies. The power spectra are obtained through the transfer functions of the modified CMBFAST code (CF), i.e. by using the linear evolution. We do not show as target the CDM power spectrum because almost indistinguishable from the baryon one. The power spectra are normalized at $`k=0.01`$ h/Mpc as in Sandvik et al. (2002) and the results in Fig. 8 can be directly compared with those of Fig. 1 in Sandvik et al.. It is very interesting to note the role of baryons, which were neglected by Sandvik et al.: baryons keep on clustering on all scales, at all time after decoupling. Even if baryons are a subdominant component ($`\mathrm{\Omega }_b=0.05`$), they are very important when one considers the total matter power spectrum $`P_T(k)=\mathrm{\Omega }_b^2P_b(k)+\mathrm{\Omega }_X^2P_X(k)`$. The inclusion of baryons smooths out the oscillations of the GCG component in the total matter power spectrum, in particular for $`\alpha 10^3`$: when comparing to the LSS data, only the shape is relevant and not the features of the GCG component. We also note that for $`\alpha 10^2`$ the appearance of a finite Jeans length and/or of a suppression of DE component maybe very interesting in connection with the bump present in APM data. Therefore, it appears that an intermediate value of $`\alpha `$ ($`0.01\alpha 0.2`$) may fit both CMB and LSS data.

On concluding, in UDM models the comparison of the total matter power spectrum with LSS changes by including baryons, even within the linear approximation. In the context of DE models, this effect is even more important since CDM perturbations clustering on all scales are also present.

## 8 Conclusions

We have analyzed quantitatively the observational effects of a GCG as a DE candidate. We have restricted our analysis to a flat universe and to purely adiabatic pressure perturbations.

So far, most authors studied only the comparison of the Hubble law predicted in presence of a GCG with the SNIa (Fabris et al. 2002a; Avelino et al. 2002; Makler et al. 2003). However, as also our Figs. 6–7 show, SNIa data constrain GCG models very weakly. For instance, the Hubble law is almost insensitive to $`\alpha `$, as we see in Fig. 6. Instead, $`\alpha `$ is the maximum value which the sound speed can reach and it is related to the time at which the Jeans instability disappears and therefore fluctuations are very sensitive to this quantity: this results in a strong dependence of the CMB spectrum on $`\alpha `$, as we see from Fig. 2.

The combined result of CMB and SNIa data is presented in Fig. 7, which leads to $`\alpha <0.2`$ and $`w_X<0.8`$ at $`95\%`$ c.l. The bound imposed by WMAP data on $`\alpha `$ is 5–10 stronger than the previous data (Bean & Dore 2003), and the Chaplygin gas ($`\alpha =1`$) is ruled out as a DE candidate at more than the $`99.99\%`$ c.l.

The possibility of a unified dark model with a GCG (corresponding to $`\mathrm{\Omega }_c=0`$) seems disfavoured by the latest CMB anisotropy data in comparison to the GCG playing the role of a DE model, although the UDM on its own is not a bad fit to the data. The $`\chi ^2`$ statistics is 987/893 for a UDM model, whereas is 979/892 for a DE model. Our result is in contrast with the analysis of the peaks position recently performed by Bento, Bertolami and Sen (2003), in which a different region in the ($`\alpha ,w_X`$) plane is allowed. This disagreement may come from an inaccurate analytic estimate of the peaks position (already pointed out by CF).

We have also addressed the problem of the comparison with the LSS data. We have compared results on the matter power spectrum with 2dF and APM data, within the linear approximation. By considering the worst case, i.e. a UDM model, we find that LSS data may be complementary to CMB data. The main question about the validity of linear treatment on such short scales however remains.

On concluding, standard quintessence models seem the most economic variation of DE models with respect to $`\mathrm{\Lambda }`$CDM. Differentiating among DE models by studying the sound speed of its perturbations seems a promising avenue and this has already been shown in two different class of models (CF; DeDeo et al. 2003). We have shown that GCG models are models with high predictive power when compared with CMB anisotropies. Current and future 111http://astro.estec.esa.nl/Planck high precision CMB anisotropy data can finally constrain and/or rule out physical scenarios based on different models of DE.

## Acknowledgements

Most of the computations have been performed on a 128-processors Linux cluster at CINECA. We thank the staff at CINECA for support. F. F. would like to thank R. Abramo for many discussions on this topic. L. A. and F. F. would like to thank O. Bertolami and I. Waga for useful discussions.