text:
# Supersonic turbulence and structure of interstellar molecular clouds

## Abstract

The interstellar medium (ISM) provides a unique laboratory for highly supersonic, driven hydrodynamics turbulence. We present a theory of such turbulence, confirm it by numerical simulations, and use the results to explain observational properties of interstellar molecular clouds, the regions where stars are born.

1. Introduction. Stars are formed as a result of gravitation (Jeans) collapse of dense clumps in interstellar molecular clouds. The structure of such clouds in a large interval of scales (from about $`100pc`$ to $`0.01pc`$) lacks any characteristic length, and can be understood as arising from supersonic hydrodynamic motions sustained on large scales by supernovae explosions . A fluid motion with characteristic large-scale relative velocities of order 1-10 km/sec compresses rapidly radiating and therefore relatively cold molecular gas ($`T10K`$) up to very high densities (above $`10^4cm^3`$). Instabilities of shock fronts create a hierarchy of gas clumps with broad mass, size, and velocity distribution controlled by the Mach number ($`M=v/c`$) and by the Alfvénic Mach number ($`M_a=v/v_a`$), where $`v_a`$ is the Alfvén velocity, $`v_a=B/\sqrt{4\pi \rho }`$, and $`c`$ is the sound speed. Depending on parameters of clouds, the Mach number can be about 30 on the largest scales, and the Alfvénic Mach number can exceed 1 as well . A systematic study of scaling properties of supersonic turbulence and of the structure of molecular clouds was initiated by the work of Larson , but despite the large number of observational and numerical investigations that appeared for the last 20 years, the theoretical understanding of the turbulence has been rather poor.

In the present paper we provide a new analytical model for a supersonic turbulent cascade and test the model against observations and numerical simulations. We find a very good agreement within the error bar uncertainty. The analytical approach is suggested by the following two numerical results. First of all, we establish that for large Mach numbers ($`M>2`$) the velocity field in the inertial interval is mostly divergence-free and shear-dominated, with the intensity of its potential component being only about 10% of the intensity of its solenoidal part. The second finding is that the most intense dissipative structures of the turbulence are two-dimensional sheets or shocks as oppose to the incompressible case where most of energy is dissipated in filaments, see also . Using these two ingredients in the framework of the so-called She-Lévêque model of turbulence we calculate two-point correlators of velocity and density fields. This allows us to construct the multifractal distributions of velocity and density fields, that statistically describe the structure of a turbulent molecular cloud. In the present paper we present the analytic derivation of the model and summarize the most important numerical results. The detailed numerical and observational analysis will appear elsewhere .

In the following section we construct the multifractal distribution of the velocity field, and numerically check the first 10 velocity-difference structure functions \[the definition is given below\]. In section 3 we proceed with the density distribution, derive a general formula for two-point density correlators and compare the result with the numerical simulations. Conclusions, applications, and future research are outlined in section 4.

2. Multifractal model of supersonic turbulence. In 1994 She and Lévêque suggested a model that turned out to be very successful in explaining the experimental findings for incompressible turbulence . The model represents a turbulent cascade as an infinitely-divisible log-Poisson process , and has three input parameters. The first two of them are the so-called naive, i.e., non-intermittent, scaling exponents of velocity fluctuation and of the eddy turnover time. The non-intermittent velocity of the eddy of size $`l`$ scales as $`v_ll^\mathrm{\Theta }`$, and the “eddy turnover” time scales as $`t_ll^\mathrm{\Delta }`$. The third parameter is the dimension of the most intense dissipative structure, $`D`$. In the incompressible case, the first two parameters take their Kolmogorov values, $`\mathrm{\Delta }=2/3`$, $`\mathrm{\Theta }=1/3`$, and the third parameter is $`D=1`$ since the dissipation mostly occurs in elongated filaments.

The model predicts the so-called structure functions of the velocity field that are defined as follows

$`S_n(l)=|v(x+l)v(x)|l^{\zeta (n)}.`$ (1)
The velocity components in this formula can be either parallel or perpendicular to vector $`l`$. In the former case the structure function is called longitudinal, in the latter case transversal. There is experimental evidence that both scale in the same way , therefore, we will not distinguish between them as far as the scaling is concerned. The She-Lévêque formula gives the following expression for the scaling exponents of the structure functions,

$`\zeta (n)=\mathrm{\Theta }(1\mathrm{\Delta })n+(3D)(1\mathrm{\Sigma }^{\mathrm{\Theta }n}),`$ (2)
where $`\mathrm{\Sigma }=1\mathrm{\Delta }/(3D)`$ . This expression has been experimentally checked to work for structure functions up to the 10th order, within an accuracy of a few percent .

To address the supersonic case, we note that in the inertial interval the turbulence is still mostly divergence-free (Fig. (1)). The effect of vorticity generation in the random 3D flow is analogous to the effect of magnetic dynamo, since magnetic field and vorticity obey the same dynamic equation. We therefore leave the Kolmogorov parameters $`\mathrm{\Theta }`$ and $`\mathrm{\Delta }`$ unchanged.

However, the most dissipative structures of a supersonic flow are different from the incompressible case. Instead of filaments they look like shocks or two-dimensional dissipative sheets, therefore, $`D=2`$ . With such input parameters, formula (2) is recast as

$`\zeta (n)=n/9+1\left(1/3\right)^{n/3}.`$ (3)
Quite remarkably, this formula works with good accuracy for the numerical simulations, fitting the first ten structure functions with an error of about 5%, see Fig.(2).

The observational data consistently indicate steeper than Kolmogorov velocity spectra, see e.g. . Our model (3) predicts $`|v_k|^2k^{1\zeta (2)}=k^{1.74}`$, which agrees with observations within error bars. Interestingly enough, the average velocity spectrum, originally inferred from observations by Larson, was $`k^{1.74}`$ on scales of order $`1pc<l<1000pc`$, and $`k^{1.76}`$ towards smaller scales, $`0.1pc<l<100pc`$ .

On the analytical side, the model (3) implies the so-called multi-fractal distribution of the turbulent fluctuations. Here we discuss the statistics of the velocity field, and in the next section, apply the results to the density field. To visualize the model, assume that the whole space (a molecular cloud or a simulation domain) contains turbulent structures of different (in general, fractal) dimensions. In the vicinity of a particular structure, the velocity difference scales with some particular exponent that we denote by $`h`$, i.e., $`v_ll^h`$. The dimension of this fractal structure will be denoted $`D(h)`$. If we divide the space into small boxes of size $`l`$, then the number of boxes covering the fractal structure of dimension $`D`$ is proportional to $`l^D`$, while the total number of boxes is proportional to $`l^3`$. The probability to find ourselves inside a box covering the fractal with dimension $`D(h)`$ is therefore $`p_h(l)l^{3D(h)}`$. To average the $`n`$’s moment of the velocity difference we just need to sum $`l^{nh}p_h(l)`$, which is the contribution of one particular fractal structure, over all the fractal structures. We get

$`S_n(l){\displaystyle \underset{h}{}}l^{nh+3D(h)}l^{\zeta (n)}.`$ (4)
Knowing $`D(h)`$ is equivalent to knowing $`\zeta (n)`$, these two functions are related by the Legendre transform as one immediately gets by assuming that $`l`$ is small, and by evaluating the sum in (4) by the steepest descent method . For example, knowing $`\zeta (n)`$ from (3), one can restore $`D(h)`$; we, however, will not need $`D(h)`$ for our present purposes.

3. Multifractal distribution of the density field. First we would like to derive an important relation of supersonic turbulence. We start with the Navier-Stokes and continuity equations:

$`_t𝐮+(𝐮)𝐮=(\eta /\rho )\mathrm{\Delta }𝐮c^2\rho /\rho +𝐟,`$ (5)
$`_t\rho +(\rho 𝐮)=0.`$ (6)
Let us introduce the density correlator $`R(t,x)=\rho (x_1,t)\rho (x_2,t)`$, and the density-weighted second-order velocity structure function, $`G^{ik}(t,x)=\rho (1)\rho (2)[u(1)u(2)]^i[u(1)u(2)]^k`$, where $`𝐱=𝐱_1𝐱_2`$. Averaging is performed over the random force. Differentiating $`R(t,x)`$ with respect to time, and using (5) and (6), one gets:

$`_t^2R+_i_kG^{ik}c^2\mathrm{\Delta }R+_i[f(1)f(2)]^i\rho (1)\rho (2)`$ (7)
$`=2\eta \mathrm{\Delta }_iu(1)^i\rho (2),`$ (8)
where the spatial derivatives are taken with respect to $`𝐱`$. In the inertial interval, the forcing and the viscous terms are small. In the supersonic regime, one can also neglect the $`c^2`$ term. Assuming now that the turbulence is in a steady state, we are left with a simple equation that must hold in the inertial interval, $`_i_kG^{ik}=0`$. Due to spatial isotropy we get:

$`G^{ik}(x)=A\delta ^{ik}+O(x^2/L^2)+\mathrm{},`$ (9)
where $`L`$ is the external force correlation length, and $`A`$ is some constant. In the inertial interval, $`xL`$, one gets $`G=const`$. To obtain the density distribution let us make a natural assumption that both the density and the velocity fields have fixed scalings in the vicinity of the same turbulent structures. In other words, close to a fractal structure where the velocity field has scaling $`h`$, the density field has some other, but also constant along the same structure, scaling $`\alpha (h)`$, i.e., $`\rho (l)l^{\alpha (h)}`$. The condition (9) now reads

$`G{\displaystyle \underset{h}{}}l^{2h+2\alpha (h)+3D(h)}=const.`$ (10)
The other restriction comes from the mass conservation law,

$`\rho {\displaystyle \underset{h}{}}l^{\alpha (h)+3D(h)}=const.`$ (11)
Strictly speaking, our constraint conditions (10) and (11) are to a certain extent phenomenological. However, we found them to be consistent with our simulations. Moreover, the theory based on them predicts density correlators rather successfully, which we are going to demonstrate now. Let us assume that the function $`\alpha (h)`$ is analytic and can be expanded as $`\alpha (h)=a+bh+gh^2+\mathrm{}`$. As a minimal model, consider the case, when $`\alpha (h)`$ is a linear function, $`\alpha (h)=a+bh`$. It turns out that such a linear ansatz is consistent with both restriction conditions (10) and (11). As follows from (4), the mass conservation condition (11) is satisfied with $`a=\zeta (p)`$ and $`b=p`$, where $`p`$ is arbitrary. The equation for $`p`$ is then derived from the dynamic constraint (10), which gives $`\zeta (2p+2)=2\zeta (p)`$. The solution of this equation will be denoted as $`p_0`$. If we use our formula (3), $`p_0`$ can be found exactly; we thus obtain $`b=p_0=2.28`$ and $`a=\zeta (p_0)=0.82`$. The multifractal distribution for the density field, $`D_\rho (\alpha )`$, is thus related to the multifractal distribution of the velocity field, $`D_\rho (\alpha )=D\left[(\alpha +\zeta (p_0))/p_0\right]`$. Fractal and multifractal distributions of density fields have indeed been inferred from observations, see, e.g., , and from numerical simulations .

By analogy with the velocity field, the quantities of practical interest are density correlators. They can be calculated with the aid of a formula analogous to (4),

$`[\rho (x+l)\rho (x)]^m{\displaystyle \underset{h}{}}l^{2m\alpha (h)+3D(h)}l^{\xi (m)}.`$ (12)
Upon substituting the linear expression for $`\alpha (h)`$ and using formula (4), one immediately gets $`\xi (m)=\zeta (2mp_0)2m\zeta (p_0)`$. This formula allows one to obtain the density scaling if the velocity structure functions are known either from theory or from experiment. Since $`\zeta (n)`$ is a concave function, $`\xi (m)`$ is negative for $`m>1/2`$. For $`m=1,2,3`$ the formula gives $`\xi (1)0.3`$, $`\xi (2)1.3`$, and $`\xi (3)2.4`$, values close to what is obtained in the numerics, see Fig. (3). We give here only the first three exponents since starting from $`m=3`$, the density exponents depend on velocity structure functions of order higher than 13, which cannot be reliably produced with our numerical resolution.

4. Conclusions. We have suggested a self-consistent model that provides an explanation for numerical and observational scaling laws of supersonic ISM turbulence, the so-called Larson’s laws. We would like to conclude with the following remarks:

1. The She-Leveque approach was also applied to incompressible MHD turbulence in , where different scalings were suggested. Most successful was the approach of Müller and Biskamp , where the energy cascade was assumed to be Kolmogorov-like, but the dissipation occurred in micro current sheets. In this case, the same formula (3) gave a good agreement with numerical simulations for structure functions up to order 8. Our results together with those by Müller and Biskamp support the ideas put forward in  and , that completely different turbulent systems can belong to the same class of universality, i.e. have the same velocity scaling exponents.
2. Turbulence with small pressure is usually referred to as Burgers turbulence, the theory of which has been rapidly developing in recent years . However, the Burgers velocity field is usually assumed to be potential, which is true in one and two dimensions, but inconsistent with the 3D case due to strong vorticity generation. However, our general relation, $`\rho (1)\rho (2)|v(1)v(2)|^{2p_0}/l^{2\zeta (p_0)}`$, which is valid inside any correlation function, can be useful for the closure problems of Burgers turbulence.
3. Our model is consistent with available observational results, although the error bars of observed velocity structure functions are too large for a precise comparison. Moreover, only projected quantities (i.e., integrated along the line of sight) are observationally available, and therefore the 3D results should be reformulated for these projected fields — this is a subject of future work .


format: text

text:
# Comment on “Analysis of chaotic motion and its shape dependence in a generalized piecewise linear map”

## Abstract

Rajagopalan and Sabir \[Phys. Rev. E 63, 057201 (2001)\] recently discussed deterministic diffusion in a piecewise linear map using an approach developed by Fujisaka et al. We first show that they rederived the random walk formula for the diffusion coefficient, which is known to be the exact result for maps of Bernoulli type since the work of Fujisaka and Grossmann \[Z. Physik B 48, 261 (1982)\]. However, this correct solution is at variance to the diffusion coefficient curve presented in their paper. Referring to another existing approach based on Markov partitions, we answer the question posed by the authors regarding solutions for more general parameter values by recalling the finding of a fractal diffusion coefficient. We finally argue that their model is not suitable for studying intermittent behavior, in contrast to what was suggested in their paper.

The study of deterministic diffusion in simple chaotic maps on the line appears to have originated about twenty years ago (see, e.g., Ref.  and further references therein). Already in the seminal work by Fujisaka and Grossmann , a variety of piecewise linear models was defined and analyzed by means of stochastic modeling. All these maps are of the form $`x_{n+1}=M_h(x_n)`$, where $`hN`$ is a control parameter, and $`x_n`$ is the position of a point particle at discrete time $`n`$. $`M_h(x)`$ is continued periodically beyond the interval $`[0,1)`$ onto the real line by a lift of degree one, $`M_h(x+1)=M_h(x)+1`$. The map defined in Ref. , which is sketched again in Fig. 1, provides a straightforward generalization of the one introduced in Ref. , which is recovered at $`h=1/2`$. For this type of maps, indeed a vast literature exists on how to obtain exact analytical results at specific cases of parameter values; Refs.  sumarize some of these methods, with more complete lists of references therein. It is furthermore well-known that the calculations are particularly simple if the parameter is such that the map exhibits the Bernoulli property .

We first wish to present a considerable shortcut to the diffusion coefficient calculations published in Ref. . Based on a theory which appears to be a precursor of what was called “Fujisaka’s characteristic function method” in Ref. , Fujisaka and Grossmann have shown  that the diffusion coefficient formula

$$D=\frac{<j^2(x_n)>}{2}$$
(1)
provides the exact solution for types of maps as the one studied in Ref. , i.e., if they share the Bernoulli property. Here $`j(x_n)`$ is the jump velocity defined as $`j(x_n):=[x_{n+1}][x_n]`$ with $`[x]`$ being the largest integer less than $`x`$, and $`<\mathrm{}>`$ denotes the average over the invariant probability density. This expression is just identical to the familiar random walk formula for diffusion on a one-dimensional lattice, where the length of jumps squared is weighted with the probability to perform such jumps (see Ref.  and further references therein). Evaluating this equation for the map under consideration leads to

$$D(h,r)=\underset{j=1}{\overset{h}{}}j^2\delta (j,r)$$
(2)
with $`hN,0<r<1`$, where $`2\delta (j,r)`$ denotes the probability to jump over a distance of $`j`$ steps and is easily calculated to

$$\delta (j,r)=\frac{2}{m_j},$$
(3)
$`m_j`$ being the slopes of the map. Combining the above two equations yields Eq. (17) of Ref. . We conclude that Rajagopalan and Sabir have confirmed again Eq. (1) of Fujisaka and Grossmann as applied to their specific map. We now focus on the author’s special case of the map defined by the relation for the slopes

$$m_0=3+\frac{4(1r^h)}{r^h(1r)}$$
(4)
with $`m_i/m_{i1}=r`$. Solutions for Eqs. (2), (3) under this constraint are shown in Fig. 2 for different $`h`$. This figure corrects the erroneous result shown in Fig. 3 of Ref. , which only includes a few data points and appears to indicate a rather irregular curve for the diffusion coefficient at $`h=2`$. Below we will explain why all the curves shown in Fig. 2 must indeed be simple functions of $`r`$.

However, first we would like to recall a second method which is not restricted to special cases of parameters such as integer heights, in contrast to the one outlined in Refs. . The basic idea of this method is to directly solve the Frobenius-Perron equation of the dynamical system,

$$\rho _{n+1}(x)=𝑑y\rho _n(y)\delta (xM_h(y)),$$
(5)
where $`\rho _n(x)`$ is the probability density for points on the real line. There exists a dense set of parameter values $`h`$ for which one can construct Markov partitions of the map, and for each of these parameter values this equation can be written as a matrix equation ,

$$𝝆_{n+1}=T(h,r)𝝆_n.$$
(6)
$`𝝆_n`$ represents a column vector of the probability densities defined on each part of the Markov partition at time $`n`$, and $`T(h,r)`$ is a topological transition matrix which can be constructed from the Markov partition. This setup provides two ways of solution: one way is to solve the eigenvalue problem of $`T(h,r)`$ and to relate the diffusion coefficient to its eigenvalues. As is shown in detail in Refs. , in special cases all calculations can be performed analytically. For the simple map defined in Ref.  these calculations are straightforward and confirm again Eqs. (2), (3). For more general cases, the matrix equation can simply be iterated  yielding numerically exact solutions for the probability density vector $`𝝆_n`$ at any time step $`n`$, as well as for any other dynamical quantity based on probability density averages. Both such methods were previously applied to various examples of piecewise linear maps . Fig. 1 presents analogous results for the map studied in Ref.  at $`h=2`$ and $`r=0.5`$, cp. to Fig. 3.1 on p.54 of Ref. . The probability density is a Gaussian on a coarse scale, whereas the fine scale is determined by the invariant density of the map on the unit interval. These deviations from an exact Gaussian can quantitatively be evaluated, e.g., by calculating the curtosis of the respective map density; for a more detailed discussion of such aspects we refer to Chapter 3 of Ref. . This interplay between fine and coarse structure of the probability densities was furthermore discussed in terms of the spectrum of eigenmodes of the Frobenius-Perron operator, see . These known results appear to be recovered in Ref.  by a respective analysis of the fluctuation spectrum, which provides an alternative way to look at the probability density of the map.

In their outlook to further work, the authors of Ref.  raised the question of how to compute the diffusion coefficient for maps with fractional heights $`h`$, and how it may look like. The application of the arsenal of methods outlined above has already given a full answer to this problem. As a central result, it was found that the diffusion coefficient for these maps is a fractal function of the parameter $`h`$. To present an example, Fig. 3 depicts the result for a mirrored zigzag-map with uniform slope, which has some similarities with the one studied in Ref. ; for details see Refs. . Knowing these results, it is straightforward to conclude that, for arbitrary height $`h`$, the map studied in Ref.  will just yield another fractal diffusion coefficient; further evidence for that statement is provided by the numerical and analytical data presented in Ref. .

So why can the diffusion coefficient of the map in Ref.  not be fractal as a function of $`r`$ at integer values of the height? One way to look at this problem is to inquire how the topology of the map is affected by parameter variation. A fundamental tool providing detailed information about the topology of a dynamical system are Markov partitions. Varying $`r`$ at integer heights does not change the Markov partition, thus the topology of the map does not change, and any quantity resulting from an average over the invariant density is a simple function of the parameter . However, changing the height changes the Markov partition in a complicated way and reflects the topological instability of the map under this type of parameter variation. This topological instability results in fractal transport coefficients.

Finally, we comment on the conclusion of Rajagopalan and Sabir that the map studied in their paper is “suited in describing diffusion systems showing intermittency”. In this aspect the authors appear to follow Ref. , where the map shown in Fig. 1 at $`h=1/2`$ was introduced for the purpose of modeling “strong correlations between successive steps… as realized in Brownian motion with directional persistence”. Indeed, Grossmann and Thomae revealed a persistent dynamics which they characterized as “intermittent-like” behavior. They linked these correlations to deviations from a pure Gaussian probability density such as the ones discussed above.

In the following we use the term “intermittency” in the sense of Pomeau and Manneville (see, e.g., Ref.  for a tutorial about their results). Particularly, we wish to distinguish it from the denotation “intermittency-like” in the sense of Grossmann and Thomae. Extensive studies of diffusion in one-dimensional intermittent maps led to the conclusion that, generally, in this case a diffusion coefficient does not exist . Furthermore, all maps studied in these references are inherently nonlinear. For piecewise linear expanding maps which are uniquely ergodic if restricted to compact spaces, such as the one of Refs. , there is no evidence for intermittency nor for anomalous diffusion . Applying the concept of conjugacy enables to transform piecewise linear maps onto nonlinear ones. However, the diffusive dynamics is invariant under conjugacy , thus the corresponding nonlinear map is again non-intermittent and normal diffusive. To our knowledge the only piecewise linear map exhibiting intermittency was introduced in Ref. , and it belongs to a very different class than the one of Refs. .

In summary, by relating the piecewise linear map studied in Ref.  to intermittent behavior the authors confuse the meaning of intermittency, in the sense of Pomeau and Manneville, with the existence of intermittent-like behavior, in the sense of persistence in the diffusive motion. Intermittency generally leads to anomalous diffusion, whereas persistence in piecewise linear maps shows up in form of local extrema of the fractal diffusion coefficient at integer and half-integer heights, see Fig. 3. We conclude that the analysis of chaotic motion and its shape dependence as performed in Ref.  has nothing to do with intermittency, but instead recovers features of the parameter-dependent normal diffusion coefficient as studied in Refs. .

The author thanks N.Korabel and J.R.Dorfman for helpful remarks.


text:
# On the Degeneracy Inherent in Observational Determination of the Dark Energy Equation of State

## I Introduction

Although observations of acoustic peaks in the Cosmic Microwave Background (CMB) fluctuation spectra are sensitive to the curvature of spacetime, and appear to require some form of cosmological dark energy , they may be less useful for discriminating among different equations of state for the dark energy. Various studies (e.g. ) have concluded that the signatures of different types of spatially smooth, evolving dark energy on the CMB spectra are relatively undiscriminating, in part because CMB fluctuations were formed at high redshift, when the dark energy was not a prominent constituent, although spatial fluctuations in the dark energy could lift the degeneracy considerably . By contrast, observations of sources at small to moderate redshifts ($`z\mathrm{a}\mathrm{few}`$) would probe epochs where the dark energy is prominent, leading to suggestions that observations of Type Ia supernovae  or of galaxy counts  could be used to determine the nature of the dark energy most effectively. However, the efficacy of these observational programs is controversial, as there is considerable degeneracy among the predictions of different dark energy models for, for example, the luminosity distance as a function of redshift.

One approach to the analysis of low redshift data would be to presume particular classes of models, perhaps parametrized by a series expansion of $`w(z)=p(z)/\rho (z)`$, and then attempt to constrain the model parameters by a likelihood or Bayesian method (e.g. ). A second approach is to fit the data by a parametric representation of, for example, the luminosity distance as a function of redshift, and then analyze the results to constrain the properties of the dark energy indirectly (e.g. ). Although the first approach is preferable (because it generally avoids difficulties associated with differentiating data, and allows a simpler assessment of uncertainties in derived parameters, in addition to outlining its implicit assumptions more clearly), the second approach is more useful for understanding whether or not such observations can ever yield a unique solution for the dark energy equation of state.

In this short note, the question considered is: Suppose one could analyze data from $`z\mathrm{a}\mathrm{few}`$ to determine the expansion rate of the Universe as a function of scale factor, $`H(a)`$, exactly. Could one then determine the equation of state of the dark energy exactly from these measurements alone? The answer, as we shall see, is no: a separate determination of the fraction of the Universe in the form of nonrelativistic particles is needed .

## II A Simple Illustrative Model

The necessity of separately determining the density of nonrelativistic particles in order to deduce the equation of state of the dark energy can be demonstrated most easily by constructing a specific example. Suppose the data show that the Universe is expanding according to the simple relationship

$$H^2=H_0^2\left(\mathrm{\Omega }_1+\frac{\mathrm{\Omega }_2}{a^3}\right),$$
(1)
where $`a(1+z)^1`$. One (tempting) interpretation of this result would be that the Universe consists of two components, nonrelativistic matter contributing a fraction $`\mathrm{\Omega }_2`$ of the closure density, and a cosmological constant contributing a fraction $`\mathrm{\Omega }_1=1\mathrm{\Omega }_2`$. However, this interpretation is not unique: there are models involving a scalar field $`\varphi `$ with nonconstant effective potentials $`V(\varphi )`$ that can lead to Eq. (1).

To see this, we can construct an explicit model with both a scalar field and nonrelativistic matter, so that

$$H^2=\frac{8\pi G\rho _\varphi }{3}+\frac{H_0^2\mathrm{\Omega }_M}{a^3}$$
(2)
with $`\mathrm{\Omega }_M\mathrm{\Omega }_2`$ in general. The scalar field energy density is therefore

$$\rho _\varphi =\rho _0\left(\mathrm{\Omega }_1+\frac{\mathrm{\Omega }_2\mathrm{\Omega }_M}{a^3}\right),$$
(3)
where $`\rho _0=3H_0^2/8\pi G`$. Differentiate $`\rho _\varphi `$ with respect to $`a`$ to find

$$\frac{d\rho _\varphi }{da}=\frac{3\rho _0(\mathrm{\Omega }_2\mathrm{\Omega }_M)}{a^4}=\frac{3(\rho _\varphi +P_\varphi )}{a}$$
(4)
where the second equality follows from conservation of energy for the scalar field. Thus, we find simply that $`P_\varphi =\rho _0\mathrm{\Omega }_1`$, and

$`\dot{\varphi }^2`$ $`=`$ $`\rho _\varphi +P_\varphi ={\displaystyle \frac{\rho _0(\mathrm{\Omega }_2\mathrm{\Omega }_M)}{a^3}}`$ (5)
$`V(\varphi )`$ $`=`$ $`{\displaystyle \frac{\rho _\varphi P_\varphi }{2}}=\rho _0\left(\mathrm{\Omega }_1+{\displaystyle \frac{(\mathrm{\Omega }_2\mathrm{\Omega }_M)}{2a^3}}\right).`$ (6)
Note that the first of Eqs. (6) requires that $`\mathrm{\Omega }_2\mathrm{\Omega }_M`$. We can determine $`V(\varphi )`$ by solving the first of Eqs. (6) for $`\varphi (a)`$, inverting to find $`a(\varphi )`$, and substituting into the second of Eqs. (6). Combining Eqs. (1) and the first of Eqs. (6) implies

$$a^2\left(\frac{d\varphi }{da}\right)^2=\frac{\dot{\varphi }^2}{H^2}=\frac{3(\mathrm{\Omega }_2\mathrm{\Omega }_M)}{8\pi G\mathrm{\Omega }_2\left(1+\mathrm{\Omega }_1a^3/\mathrm{\Omega }_2\right)}.$$
(7)
Define

$$x=\frac{\mathrm{\Omega }_1a^3}{\mathrm{\Omega }_2}\mathrm{and}\varphi _g^2=\frac{\mathrm{\Omega }_2\mathrm{\Omega }_M}{24\pi G\mathrm{\Omega }_2};$$
(8)
then if $`\psi \varphi /\varphi _g`$

$$x^2\left(\frac{d\psi }{dx}\right)^2=\frac{1}{1+x},$$
(9)
which has the general solution

$$\psi =\psi _{\mathrm{}}\pm 2\mathrm{sinh}^1\left(\frac{1}{\sqrt{x}}\right),$$
(10)
where $`\psi _{\mathrm{}}`$ is a constant of integration. The sign in this solution can be chosen arbitrarily; chosing the $``$ sign gives a solution in which $`\psi `$ increases with increasing $`a`$. If we let $`\psi =\psi _0`$ today, then we find

$$\psi =\psi _0+2\mathrm{sinh}^1\sqrt{\frac{\mathrm{\Omega }_2}{\mathrm{\Omega }_1}}2\mathrm{sinh}^1\sqrt{\frac{\mathrm{\Omega }_2}{\mathrm{\Omega }_1a^3}},$$
(11)
which we may invert to find

$$\frac{1}{a^3}=\left[\mathrm{cosh}\left(\frac{\widehat{\psi }}{2}\right)\sqrt{\frac{1}{\mathrm{\Omega }_2}}\mathrm{sinh}\left(\frac{\widehat{\psi }}{2}\right)\right]^2,$$
(12)
where $`\widehat{\psi }\psi \psi _0`$. The second of Eqs. (6) then determines the potential to be

$$V(\varphi )=\rho _0\left\{\mathrm{\Omega }_1+\frac{1}{2}(\mathrm{\Omega }_2\mathrm{\Omega }_M)\left[\mathrm{cosh}\left(\frac{\widehat{\psi }}{2}\right)\sqrt{\frac{1}{\mathrm{\Omega }_2}}\mathrm{sinh}\left(\frac{\widehat{\psi }}{2}\right)\right]^2\right\}.$$
(13)
Note that $`\widehat{\psi }2\mathrm{coth}^1\sqrt{1/\mathrm{\Omega }_2}`$ and $`V(\varphi )\rho _0\mathrm{\Omega }_1`$ asymptotically. Also, when $`\mathrm{\Omega }_2=1\mathrm{\Omega }_11`$, the potential becomes exponential: $`V(\varphi )\frac{1}{2}\rho _0(1\mathrm{\Omega }_M)e^{\widehat{\psi }}`$.

The key question here revolves around prior assumptions. The “simplest” interpretation of Eq. (1) is that the energy density of the Universe is the sum of contributions due to nonrelativistic matter and a cosmological constant. However natural this interpetation may seem to be, it only follows if one assumes either that the dark energy is in the form of a cosmological constant, or that the mass density in nonrelativistic form has closure parameter $`\mathrm{\Omega }_2`$. If one is somewhat more agnostic on either of these points, then a range of possibile interpretations accounts for the data equally well. The degeneracy is not just a property of models for which $`H^2(a)`$ is given by Eq. (1). If the data were to favor a model of the form

$$H^2(a)=H_0^2\left[\mathrm{\Omega }_1g(a)+\frac{\mathrm{\Omega }_2}{a^3}\right],$$
(14)
with $`g(a)`$ some function, and $`\mathrm{\Omega }_1g(1)+\mathrm{\Omega }_2=1`$, then without additional prior assumptions the most one can say for sure is that $`\rho _\varphi =\rho _0\left[\mathrm{\Omega }_1g(a)+(\mathrm{\Omega }_2\mathrm{\Omega }_M)/a^3\right]`$, which leaves open the possibility of a range of different solutions for $`V(\varphi )`$.

This model illustrates that if one makes no prior assumptions about the form of the quintessence potential $`V(\varphi )`$, then it is impossible to fit low redshift data alone to find $`V(\varphi )`$ uniquely unless $`\mathrm{\Omega }_M`$ is known independently – in this example, a range of models, parametrized by the value of $`f\mathrm{\Omega }_M/\mathrm{\Omega }_2`$, will fit equally well. Constraints on the value of $`f`$ limit the range of acceptable models. Although minimally $`f1`$, we know that this model would be unable to account for the formation of large scale structure if $`f`$ is too small. If we assume that only the nonrelativistic matter can clump (i.e. that the dark energy and relativistic matter remain smooth), then, in linear theory, the density contrast $`D(a)`$ obeys the equation (e.g. )

$$\frac{d^2D}{da^2}+\left(3+\frac{d\mathrm{ln}H}{d\mathrm{ln}a}\right)\frac{1}{a}\frac{dD}{da}\frac{4\pi G\rho _M}{H^2a^2}D=0;$$
(15)
Eq. (15) has both growing and shrinking modes, $`D_+(a)`$ and $`D_{}(a)`$, respectively. If we include radiation (which is a known component of the Universe, and subdominant today, but important in determining linear theory growth rates), then Eq. (1) must be altered to

$$H^2(a)=H_0^2\left(\mathrm{\Omega }_1+\frac{\mathrm{\Omega }_2}{a^3}+\frac{\mathrm{\Omega }_{rad}}{a^4}\right),$$
(16)
where $`\mathrm{\Omega }_1+\mathrm{\Omega }_2+\mathrm{\Omega }_{rad}=1`$ and, for three low mass or massless neutrinos, $`\mathrm{\Omega }_{rad}h_{0.7}^28.5\times 10^5`$ if $`h_{0.7}=H_0/(70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1)`$. The growing mode mainly amplifies perturbations between $`a_{eq}=\mathrm{\Omega }_{rad}/\mathrm{\Omega }_2`$ and $`a_1=(\mathrm{\Omega }_2/\mathrm{\Omega }_1)^{1/3}`$. During this phase, $`H^2a^3`$, and $`D_+(a)a^{\sigma _+}`$, where

$$\sigma _+(f)=\frac{1}{4}\left(\sqrt{1+24f}1\right).$$
(17)
The overall linear theory growth factor is $`(a_1/a_{eq})^{\sigma _+(f)}`$, so the ratio of the growth factor for $`f1`$ to its value for $`f=1`$ is approximately

$$\left(\frac{a_1}{a_{eq}}\right)^{\frac{1}{4}\left(\sqrt{1+24f}5\right)}\left(\frac{a_1}{a_{eq}}\right)^{3(1f)/5},$$
(18)
where the approximation assumes $`1f1`$. Thus, we estimate that for the growth factor to be within a factor of two of its value for $`f=1`$, we must have $`1f5\mathrm{ln}2/3\mathrm{ln}(a_1/a_{eq})0.12`$.

The solid line in Fig. 1 shows the ratio of the growth factor for $`\mathrm{\Omega }_M\mathrm{\Omega }_2`$ to its value for $`\mathrm{\Omega }_2=\mathrm{\Omega }_M`$, that is

$$R_{\delta \rho /\rho }\frac{D_+(1;\mathrm{\Omega }_M,\mathrm{\Omega }_2)}{D_+(1;\mathrm{\Omega }_2,\mathrm{\Omega }_2)}$$
(19)
as a function of $`\mathrm{\Omega }_M`$ for $`\mathrm{\Omega }_2=0.5`$, where $`D_+(a;\mathrm{\Omega }_M,\mathrm{\Omega }_2)`$ is the linear theory growth factor evaluated at scale factor $`a`$, given $`\mathrm{\Omega }_M`$ and $`\mathrm{\Omega }_2`$. These results were found by solving Eq. (15) numerically from $`aa_{eq}`$ to the present day, assuming identical initial perturbation amplitudes. (Relativistic particles were included as a smooth component via Eq. .) The numerical results show that for $`R_{\delta \rho /\rho }0.5`$, then we must require $`\mathrm{\Omega }_M0.44`$, corresponding to $`f0.88`$, approximately the same as the analytic estimate of the previous paragraph. (The estimate is accurate to about 20% for $`\mathrm{\Omega }_20.1`$.)

Linear theory also predicts velocity fluctuations proportional to the peculiar velocity factor $`T_+(a;\mathrm{\Omega }_M,\mathrm{\Omega }_2)=d\mathrm{ln}D_+(a;\mathrm{\Omega }_M,\mathrm{\Omega }_2)/d\mathrm{ln}a`$. The dashed line in Fig. 1 shows

$$R_{\delta v}=\frac{T_+(1;\mathrm{\Omega }_M,\mathrm{\Omega }_2)}{T_+(1;\mathrm{\Omega }_2,\mathrm{\Omega }_2)},$$
(20)
the ratio of the peculiar velocity factor at $`a=1`$ for $`\mathrm{\Omega }_M\mathrm{\Omega }_2`$ to its value for $`\mathrm{\Omega }_M=\mathrm{\Omega }_2`$. Although there is some dependence of $`R_{\delta v}`$ on $`\mathrm{\Omega }_M`$, it is not as extreme as for $`R_{\delta \rho /\rho }`$.

Thus, primarily from the requirement that linear density perturbations can grow substantially from small values in the early Universe, we know that if Eq. (1) were truly exact, $`f`$ must be close to one. Imposing limits on the range of plausible $`\mathrm{\Omega }_M`$ from independent, physical and phenomenological arguments amounts to using prior information to restrict the possibilities to be compared with the data from observations at $`z\mathrm{a}\mathrm{few}`$. (The use of prior information from large scale structure and CMB anisotropies has been discussed in e.g. .) This does not eliminate the degeneracy implicit in Eq. (1) altogether, but does diminish its significance by constraining the range of acceptable values of $`f`$ severely. However, the fact remains that if one were to try to determine the equation of state of the dark energy from observations that, hypothetically, yield Eq. (1), then one can never hope to find a unique result from an analysis that ignores constraints from the development of large scale structure (or any other considerations that yield independent information about $`\mathrm{\Omega }_M`$, such as peculiar velocities, CMB fluctuations and weak lensing). Moreover, even if such restrictions are imposed, unless $`\mathrm{\Omega }_M`$ is determined precisely by them, some degeneracy remains in the construction of $`V(\varphi )`$ from observations at $`z\mathrm{a}\mathrm{few}`$.

FIG. 1 Relative growth factor (solid) and velocity factor (dashed) of linear perturbations for models with $`\mathrm{\Omega }_M\mathrm{\Omega }_2=0.5`$ that expand according to Eq. (1).

## III Can $`\mathrm{\Omega }_M`$ and $`\mathrm{\Omega }_2`$ Differ Substantially?

If Eq. (1) with $`\mathrm{\Omega }_2`$ substantially different from $`\mathrm{\Omega }_M`$ is an excellent approximation at low to moderate redshifts, it can never be exact. In this case, there must be a transition in $`H^2(a)`$, so that the component of the density that scales like $`a^3`$ is predominantly due to nonrelativistic particles for a period sufficient to grow large scale structure. The implicit deviations from Eq. (1) would contain information on $`\mathrm{\Omega }_M`$ as well as $`\mathrm{\Omega }_2`$, and probably additional parameters as well. In this situation, it might be possible to ascertain the equation of state of the dark energy from measurements at low to moderate redshifts alone, provided that the imprint of $`\mathrm{\Omega }_M`$ and the additional parameters in $`H^2(a)`$ can be discerned from these observations. The question is, if Eq. (1) is an accurate (but inexact) fit at $`z\mathrm{a}\mathrm{few}`$, and $`\mathrm{\Omega }_2`$ differs substantially from $`\mathrm{\Omega }_M`$, how well could we determine $`\mathrm{\Omega }_M`$, and hence break the degeneracy, observationally?

To make these issues more concrete, suppose in reality Eq. (1) is actually the large $`a`$ limit of a more general, exact relationship, say

$$H^2(a)=H_>^2(a)F\left(\frac{a}{a_t}\right)+H_<^2(a)\left[1F\left(\frac{a}{a_t}\right)\right],$$
(21)
where $`F(q)0`$ as $`q0`$ and $`F(q)1`$ as $`q\mathrm{}`$, with $`a_t<1`$ a transition value of the scale factor, $`H_>^2(a)`$ given by Eq. (1), and, for example,

$$H_<^2(a)=H_0^2\left(\mathrm{\Omega }_<+\frac{\mathrm{\Omega }_M}{a^3}\right);$$
(22)
since the model is flat,

$$1=(\mathrm{\Omega }_1+\mathrm{\Omega }_2)F_0+(\mathrm{\Omega }_<+\mathrm{\Omega }_M)(1F_0),$$
(23)
where $`F_0F(1/a_t)`$. Such a model could be consistent with the growth of large scale structure, but might still be hard to distinguish from Eq. (1) based on observations at low to moderate redshift, depending on the value of $`a_t`$, and the form of $`F(a/a_t)`$. Eq. (21), although admittedly ad hoc, is useful for studying the extent to which observations of various sorts could discern that $`\mathrm{\Omega }_2\mathrm{\Omega }_M`$, and pin down the difference sufficiently to break the degeneracy discussed above.

For a flat cosmology, Eq. (21) implies that the radial coordinate for a source at $`a=(1+z)^1`$ is

$`r_S(a;\mathrm{\Omega }_M,\mathrm{\Omega }_2,F)`$ $`=`$ $`\eta (1)\eta (a)={\displaystyle _a^1}{\displaystyle \frac{db}{b^2H(b)}}`$ (24)
$`=`$ $`H_0^1{\displaystyle _a^1}{\displaystyle \frac{db}{\sqrt{b\left\{\mathrm{\Omega }_eb^3+1\mathrm{\Omega }_e+[(\mathrm{\Omega }_<\mathrm{\Omega }_1)b^3+\mathrm{\Omega }_M\mathrm{\Omega }_2][F_0F(b/a_t)]\right\}}}}`$ (25)
where $`\mathrm{\Omega }_e\mathrm{\Omega }_1F_0+\mathrm{\Omega }_<(1F_0)`$, and $`d\eta dt/a=da/a^2H(a)`$. (Here, the notation $`r_S(a;\mathrm{\Omega }_M,\mathrm{\Omega }_2,F)`$ means the radial coordinate as a function of $`a`$ given $`\mathrm{\Omega }_M`$ and $`\mathrm{\Omega }_2\mathrm{\Omega }_M`$, as well as a set of parameters required to define $`F(a/a_t)`$.) If $`F(a/a_t)`$ does not vary much over the range of $`a`$ covered by the observations, then the data would only determine $`\mathrm{\Omega }_e`$ accurately, and would yield little useful information on $`\mathrm{\Omega }_2`$ or $`\mathrm{\Omega }_M`$: only the combination $`\mathrm{\Omega }_2F_0+\mathrm{\Omega }_M(1F_0)=1\mathrm{\Omega }_e`$ could be deduced. This situation is possible provided that $`a_t`$ is fairly small, below the range covered by observations, and $`F(a/a_t)`$ varies relatively slowly as a function of $`a/a_t`$ for large values of $`a/a_t`$. In this case, we cannot expect to eliminate the degeneracy among models directly from observations of low to moderate redshift sources alone. Constraints from other sorts of observations that probe $`\mathrm{\Omega }_M`$ directly would be needed to deduce the quintessence model as accurately as possible.

More quantitative statements only can be made in terms of specific $`F(a/a_t)`$, and require some exploration of parameters. Let us examine one choice: $`\mathrm{\Omega }_1+\mathrm{\Omega }_2=1=\mathrm{\Omega }_<+\mathrm{\Omega }_M`$, in which case Eq. (23) is satisfied exactly, and Eq. (25) becomes

$$r_S(a;\mathrm{\Omega }_2,\mathrm{\Omega }_M,F)=H_0^1_a^1\frac{db}{\sqrt{b\left\{\mathrm{\Omega }_2+(1\mathrm{\Omega }_2)b^3(\mathrm{\Omega }_2\mathrm{\Omega }_M)(1b^3)\left[1F(b/a_t)\right]\right\}}}.$$
(26)
FIG. 2 Magnitude differences among models with $`\mathrm{\Omega }_2=0.5`$ for $`a_t=0.2`$ (solid) and $`a_t=0.5`$ (dashed) as functions of $`z`$ for, from the top for each $`a_t`$, $`\mathrm{\Omega }_M=0.1,0.2,0.3`$ and $`0.4`$.

Eq. (26) may be used to compute magnitude differences among models with given $`\mathrm{\Omega }_2`$, $`\mathrm{\Omega }_M`$ and $`F(a/a_t)`$; for illustrative purposes, let us choose $`1F(a/a_t)=e^{a/a_t}`$. Fig. 2 depicts the magnitude difference

$$\mathrm{\Delta }m(z;\mathrm{\Omega }_2,\mathrm{\Omega }_M,a_t)=5.0\mathrm{log}_{10}\left[\frac{r_S((1+z)^1;\mathrm{\Omega }_M,\mathrm{\Omega }_2,a_t)}{r_S((1+z)^1;\mathrm{\Omega }_2,\mathrm{\Omega }_2,a_t)}\right]$$
(27)
relative to the model with $`\mathrm{\Omega }_M=\mathrm{\Omega }_2=0.5`$ as a function of $`z`$ assuming $`a_t=0.2`$ (solid) and $`a_t=0.5`$ (dashed) for (from top to bottom in each case) $`\mathrm{\Omega }_M=0.1,0.2,0.3`$ and $`0.4`$. Not surprisingly, the magnitude differences are smaller for the larger value of $`a_t`$, and decrease as $`\mathrm{\Omega }_M\mathrm{\Omega }_2`$. Nevertheless, the largest magnitude difference shown for $`z3`$ is 0.29 mag, for $`(z,a_t,\mathrm{\Omega }_M)=(3,0.5,0.1)`$; for $`(z,a_t,\mathrm{\Omega }_M)=(3,0.2,0.1)`$, the difference is 0.077 mag.

FIG. 3 Fractional deviations of the angular size of the (sonic) horizon at $`a=10^3`$ from the $`\mathrm{\Omega }_M=\mathrm{\Omega }_2=0.5`$ model as functions of $`\mathrm{\Omega }_M\mathrm{\Omega }_2`$. The solid line is for $`a_t=0.2`$, and the dashed line for $`a_t=0.5`$. For computing the angular sizes, $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ was adopted.

Fig. 3 shows results for

$$\frac{\mathrm{\Delta }\theta }{\theta }(\mathrm{\Omega }_M;\mathrm{\Omega }_2,a_t)=\left|\frac{\eta (10^3;\mathrm{\Omega }_M,\mathrm{\Omega }_2,a_t)r_S(10^3;\mathrm{\Omega }_2,\mathrm{\Omega }_2,a_t)}{r_S(10^3;\mathrm{\Omega }_M,\mathrm{\Omega }_2,a_t)\eta (10^3;\mathrm{\Omega }_2,\mathrm{\Omega }_2,a_t)}1\right|,$$
(28)
the fractional difference between the angular size of the (sonic) horizon at $`a=10^3`$ as a function of $`\mathrm{\Omega }_M\mathrm{\Omega }_2`$ and its value for $`\mathrm{\Omega }_M=\mathrm{\Omega }_2`$ for the two different values of $`a_t`$. (Deviations of the sound speed from $`1/\sqrt{3}`$ have not been included.) Larger fractional deviations arise for smaller $`\mathrm{\Omega }_M`$ and $`a_t`$. The largest deviation shown in the two figures is 12.7% for $`(\mathrm{\Omega }_M,a_t)=(0.1,0.2)`$; the largest fractional deviation for $`a_t=0.5`$ is about 1.3%, and occurs near $`\mathrm{\Omega }_M=0.225`$.

Although this is only one example, based on a particular (ad hoc) choice of $`F(a/a_t)`$, it illustrates that while data at both low to moderate $`z`$ as well as CMB observations can yield information that lifts the degeneracy between $`\mathrm{\Omega }_2`$ and $`\mathrm{\Omega }_M`$, considerable precision may be needed to determine the equation of state of the dark energy accurately. Note that no attempt has been made to fit simulated data based on Eq. (28) to figure out how well $`H(a)`$ could be determined in practice; the comparisons made in Figs. 2 and 3 are based on exact functions. The observational challenge of distinguishing among models at these levels could be even greater than these relatively small differences would indicate. On the other hand, Figs. 2 and 3 were constructed based on a particular (ad hoc) model for $`H^2(a)`$ designed for a phenomenological study of how easily information on $`\mathrm{\Omega }_M`$ could be gleaned from luminosity distances and CMB acoustic peaks. The task could be easier or harder in cosmologies based on particular physical theories for the dark energy, depending on the details of the model.

## IV Discussion

To conclude, the point of this paper is not to suggest that it is absolutely inconceivable that the equation of state of the dark energy could be measured, but rather to show, via a specific example, how measurements at $`z3`$ admit a continuum of interpetations in terms of evolving dark energy fields unless $`\mathrm{\Omega }_M`$ can be determined separately somehow. The specific example, Eq. (1), can be intepreted “naturally” in terms of a Universe containing a mixture of nonrelativistic particle dark matter plus a cosmological constant, but can be interpreted equally well in terms of a quintessence model with a range of potentials, Eq. (13). This specific example sheds some light on why simulated analyses of low to medium redshift data appear to allow degenerate interpretations: it is always possible that the component of the total energy density of the Universe that evolves like $`a^3`$ is only partially due to nonrelativistic particles, with the rest arising from quintessence. Unique interpretations can only be obtained in terms of specific models for the quintessence that forbid such a conspiracy a priori, by either specifying the form of $`V(\varphi )`$ or fixing the value of $`\mathrm{\Omega }_M`$.

Realistically, the extent to which the $`\rho a^3`$ constituent must be due to nonrelativistic particles is constrained by the requirement that large scale structure formation evolve “normally,” that is, unimpeded by the existence of a component with density proportional to $`a^3`$ that is incapable of clustering. If Eq. (1) were truly exact, then $`\mathrm{\Omega }_2`$ would have to be close to $`\mathrm{\Omega }_M`$ for large scale structure to grow, but even in this case, a limited degeneracy remains in the determination of $`V(\varphi )`$. If $`\mathrm{\Omega }_2`$ and $`\mathrm{\Omega }_M`$ differ substantially, then Eq. (1) cannot be exact, raising the question of how well one could discern the two parameters separately from observations. To examine this issue, a modified (phenomenological) expansion law, Eq. (21), that encodes information on both $`\mathrm{\Omega }_2`$ and $`\mathrm{\Omega }_M`$ was introduced. In the context of a particular model that is consistent with large scale structure formation, but reduces to Eq. (1) with $`\mathrm{\Omega }_2\mathrm{\Omega }_M`$ at large enough $`a`$, we have seen that information about $`\mathrm{\Omega }_M`$ could be gleaned both from measurements of the positions of CMB peaks and from luminosity distance determinations. However, the deviations among models with various $`\mathrm{\Omega }_M\mathrm{\Omega }_2`$ can be quite small, which would still pose a substantial challenge for programs that aim to determine the equation of state of the dark energy.

Of course, in the context of specific models for the quintessence field, embodied in particular forms for $`V(\varphi )`$, the analyses may not encounter pronounced degeneracies. At present, though, there is little compelling reason to assume any particular $`V(\varphi )`$ a priori. When analyzing data for a particular form of $`V(\varphi )`$, one is primarily engaged in estimating the parameters of the model, as well as $`\mathrm{\Omega }_M`$. Thus, one might be able to find the best fit model of a particular type (e.g. $`V(\varphi )=\mathrm{constant}`$), without being able to tell if the underlying model would be favored by the data if other possibilities were admitted. If the allowed ranges of parameters for a particular quintessence model do not shrink with the accumulation of data, one can conclude that the model is inadequate with confidence. However, there is no guarantee that a given model is correct even if the data seem to converge on a unique set of parameters, as Eqs. (1) and (13) illustrate.

###### Acknowledgements.

I thank Eanna Flanagan, Daniel Holz, and Paul Steinhardt for helpful comments on this paper.


text:
# The few-body problem for trapped bosons with large scattering length

\[

## Abstract

We calculate energy levels of two and three bosons trapped in a harmonic oscillator potential with oscillator length $`a_{\mathrm{o}sc}`$. The atoms are assumed to interact through a short-range potential with a scattering length $`a`$, and the short-distance behavior of the three-body wave function is characterized by a parameter $`\theta `$. For large positive $`a/a_{\mathrm{o}sc}`$, the energies of states which, in the absence of the trap, correspond to three free atoms approach values independent of $`a`$ and $`\theta `$. For other states the $`\theta `$ dependence of the energy is strong, but the energy is independent of $`a`$ for $`|a/a_{\mathrm{o}sc}|1`$.

\]

By exploiting the properties of Feshbach resonances, it has become possible to tune atomic scattering lengths to essentially any value, either repulsive or attractive . This has led to renewed interest in the properties of systems with large scattering length $`a`$, a problem which has also been studied extensively in the context of nuclear physics, because the scattering lengths for nucleon-nucleon interactions are large in magnitude compared with the spatial scale $`r_0`$ of the interparticle interaction.

A novel possibility for exploring experimentally the properties of few-body systems is to trap atoms in a periodic potential created by standing-wave light beams, a so-called optical lattice. For example, in the experiments described in Ref. , up to $`200,000`$ 87Rb atoms were distributed over more than $`150,000`$ lattice sites, thereby creating a large number of few-body systems with up to 2.5 atoms each on average. This provides one motivation for theoretical studies of few-body problems in traps. A second motivation is to understand theoretically the properties of bulk matter when the magnitude of the scattering length is large compared with the interparticle spacing $`n^{1/3}`$, where $`n`$ is the particle density. In this regime, the usual low density approximations, which are valid only if $`n|a|^31`$, fail. Recently, one has sought a universal description, with $`n|a|^3`$ as the only parameter, of the many-body problem for Fermi  and Bose  systems interacting via short-range potentials ($`r_0n^{1/3}`$), but with $`n|a|^31`$. Studies of three-body recombination in ultracold Bose gases have, however, revealed that the rate depends on an additional three-body parameter , and in a recent work, it has been suggested that the equation of state of a Bose gas in the limit $`n^{1/3}|a|`$ may depend on the three-body parameter . Study of the few-body system in a trap is therefore of interest in order to obtain insights into the properties of bulk matter.

In this Letter we calculate the energy of systems with two or three identical bosons of mass $`m`$ in a harmonic trapping potential $`V_{\mathrm{t}rap}=m\omega ^2r^2/2`$, where $`\omega `$ is the frequency of oscillations of a particle about the minimum of the potential. The length associated with the zero-point motion of a particle in the trap is $`a_{\mathrm{o}sc}=\sqrt{\mathrm{}/m\omega }`$. The bosons are assumed to interact via a short-range potential. To take this into account, we impose the boundary condition $`\mathrm{ln}(r_{ij}\mathrm{\Psi })/r_{ij}|_{r_{ij}0}=1/a`$. Here $`\mathrm{\Psi }(𝐫_1,\mathrm{}𝐫_N)`$ is the wave function of the $`N`$-body system, $`𝐫_i`$ being the coordinates of the particles, and $`𝐫_{ij}=𝐫_i𝐫_j`$.

We first consider the two-body problem. Solving the Schrödinger equation for two bosons in a harmonic oscillator potential, one finds that the eigenenergies $`E`$, excluding the center of mass contribution, are given by

$$\frac{a}{a_{\mathrm{o}sc}}=\frac{1}{\sqrt{2}}\frac{\mathrm{\Gamma }[1/4E/(2\mathrm{}\omega )]}{\mathrm{\Gamma }[3/4E/(2\mathrm{}\omega )]}.$$
(1)
For $`a=0`$ the $`s`$-wave states have energies $`E=\mathrm{}\omega (2\nu +3/2)`$, ($`\nu =0,1,\mathrm{}`$) at which the gamma function in the denominator diverges, while for $`a/a_{\mathrm{o}sc}\pm \mathrm{}`$, the energies are $`E=\mathrm{}\omega (2\nu +1/2)`$. Hence, when the bosons are confined to a volume smaller than $`|a|^3`$, the scattering length becomes unimportant. Interestingly, the energy does not increase without limit as $`a`$ becomes large, as one would predict from the standard low-density result for the energy.

We now investigate the analogous problem for three trapped bosons, and adopt the method of Ref. . The wave function of three identical bosons may be decomposed in the usual Faddeev fashion as

$$\mathrm{\Psi }(𝐫_1,𝐫_2,𝐫_3)=\psi (𝐱,𝐲)+\psi (𝐱^{},𝐲^{})+\psi (𝐱^{\prime \prime },𝐲^{\prime \prime }).$$
(2)
Here the Jacobi coordinates $`𝐱=𝐫_{12}`$, and $`𝐲=𝐫_3(𝐫_1+𝐫_2)/2`$ have been used. If interactions between atoms are purely $`s`$-wave, $`\psi `$ is a function of the hyperradius $`\rho =(x^2/2+2y^2/3)^{1/2}`$ and the hyperangle $`\alpha =\mathrm{arctan}(\sqrt{3}x/2y)`$ only. Variables related to the original ones by cyclic permutations of the particles are denoted by primes and double primes.

The wave function $`\psi `$ may be represented exactly by an infinite sum

$$\psi (\rho ,\alpha )=\frac{1}{\sqrt{3}}\underset{i}{}\frac{f_i(\rho )}{\rho ^{5/2}}\frac{\varphi _i(\rho ,\alpha )}{\mathrm{sin}\alpha \mathrm{cos}\alpha },$$
(3)
where $`\varphi _i(\rho ,\alpha )`$ are eigensolutions to the hyperangular part of the problem, which for a short-range potential reduce to

$$\varphi _i(\rho ,\alpha )=\sqrt{N(\rho )}\mathrm{sin}\left[\sqrt{\lambda _i}(\alpha \pi /2)\right].$$
(4)
The eigenvalues $`\lambda _i`$ are determined by the boundary condition as $`\alpha 0`$, i.e. when two of the bosons approach each other, and are given by the equation

$$\sqrt{\lambda _i}\mathrm{cos}\left[\sqrt{\lambda _i}\frac{\pi }{2}\right]\frac{8}{\sqrt{3}}\mathrm{sin}\left[\sqrt{\lambda _i}\frac{\pi }{6}\right]=\sqrt{2}\frac{\rho }{a}\mathrm{sin}\left[\sqrt{\lambda _i}\frac{\pi }{2}\right].$$
(5)
The hyperradial wave functions $`f_i`$ are the solutions of an infinite number of coupled differential equations

$$\begin{array}{c}\left(\frac{^2}{\rho ^2}+\frac{\lambda _i(\rho /a)1/4}{\rho ^2}Q_{ii}(\rho )+\frac{\rho ^2}{a_{\mathrm{o}sc}^4}\frac{2mE}{\mathrm{}^2}\right)f_i(\rho )\hfill \\ \hfill =\underset{ji}{}\left(Q_{ij}(\rho )+2P_{ij}(\rho )\frac{}{\rho }\right)f_j(\rho ).\end{array}$$
(6)
The effective potential for the hyperangular motion therefore depends on the eigenvalue $`\lambda _i`$ associated with the hyperangular degrees of freedom. The coupling matrices are given by the kinetic-energy operator acting on the hyperangular wave functions:

$$P_{ij}(\rho )=\mathrm{\Phi }_i\left|\frac{}{\rho }\right|\mathrm{\Phi }_j,Q_{ij}(\rho )=\mathrm{\Phi }_i\left|\frac{^2}{\rho ^2}\right|\mathrm{\Phi }_j,$$
(7)
where $`\mathrm{\Phi }_i(\rho ,\alpha ,\alpha ^{},\alpha ^{\prime \prime })=\varphi _i(\rho ,\alpha )+\varphi _i(\rho ,\alpha ^{})+\varphi _i(\rho ,\alpha ^{\prime \prime })`$. An approximation which has worked well in a number of applications  is to neglect the coupling terms $`P`$ and $`Q`$. This is an adiabatic approximation in the sense that the hyperangular functions $`\mathrm{\Phi }_i(\rho ,\alpha ,\alpha ^{},\alpha ^{\prime \prime })`$ are assumed to vary much more slowly with hyperradius than with the hyperangles. The method can be refined by including coupling between a finite subset of states.

In Fig. 1 we plot $`\lambda _i(\rho /|a|)`$ calculated from Eq. (5) for some adiabatic states. For both $`a>0`$ and $`a<0`$ we show the state approaching $`\lambda _i=4`$ for large $`\rho /|a|`$. This would, in the absence of the trap, correspond at large $`\rho `$ to three free atoms, and it is the state of greatest interest in the present work. Also shown for $`a>0`$ is the result for the state which would correspond to a diatomic molecule and an atom.

Let us now examine the form of the adiabatic hyperspherical wave functions. The functions of interest are ones which remain finite as $`\rho 0`$. The effective potential for small $`\rho /|a|`$ varies as $`(\lambda _i(0)1/4)\rho ^2`$. Consequently for $`\lambda _i(0)>0`$ the hyperradial wave function $`f_i`$ tends to zero as $`\rho ^{\sqrt{\lambda _i(0)}+1/2}`$ for small $`\rho /|a|`$, and it does not depend on the nature of three-body correlations when no two atoms are separated by more than the typical range of the interaction. The situation is quite different for $`\lambda _i(0)<0`$, since the effective potential is then attractive at small $`\rho `$. The two linearly independent solutions oscillate rapidly as $`\rho 0`$, but do not diverge, and the general solution for small $`\rho `$ may be written

$$f(\rho )\sqrt{\rho }\mathrm{sin}[|s_0|\mathrm{ln}(\rho /a_{\mathrm{o}sc})+\theta ],$$
(8)
where $`s_0=\sqrt{\lambda _0(0)}1.00624i`$ and $`\theta `$ is a three-body parameter analogous to the phase shift for the two-body problem.

We turn now to the energy eigenvalues. In Fig. 2 we exhibit numerical results for the solution of Eq. (6), where the only states included are those shown in Fig. 1. For $`a>0`$, coupling between the two hyperangular functions has been included. Since this is weak, it is still possible to separate states that for large $`\rho `$ correspond to three free atoms from states that correspond to a molecule and an atom. We show the energy only of the states that for large $`\rho `$ correspond to three free atoms, since these are the states of interest here.

For small $`|a/a_{\mathrm{o}sc}|`$ the energy of the lowest harmonic oscillator state can be expanded as

$$E=\mathrm{}\omega \left[3+\frac{6}{\sqrt{2\pi }}\frac{a}{a_{\mathrm{o}sc}}+c\left(\frac{a}{a_{\mathrm{o}sc}}\right)^2+\mathrm{}\right].$$
(9)
A lower bound to the energy may be found by omitting $`Q_{ii}(\rho )`$ in Eq. (6). An upper bound may be obtained by evaluating the energy for the trial wave function $`f_i^{(0)}\mathrm{\Phi }_i`$, where $`f_i^{(0)}`$ is the solution of Eq. (6) in the absence of the coupling terms . Including the same states as before, we find $`0.55<c<0.25`$. For comparison, if three-body correlations are neglected, the correlation energy (the difference between the total energy and that of the corresponding state in the absence of interactions) is three times the correlation energy for two bodies, as obtained from Eq. (1). The resulting approximation for the energy has the same term linear in $`a/a_{\mathrm{o}sc}`$, but the quadratic term has a different coefficient, $`c=6(1\mathrm{ln}2)/\pi 0.59`$. Hence, three-body effects are to leading order proportional to $`(a/a_{\mathrm{o}sc})^2`$.

With increasing $`|a|/a_{\mathrm{osc}}`$ the dependence of the energy on $`\theta `$ grows. For $`a<0`$ this is due to the sensitivity of the adiabatic state to the inner boundary condition. The adiabatic state for $`a>0`$ is insensitive to $`\theta `$, and the $`\theta `$ dependence arises due to coupling to states having different hyperangular functions. Energies for $`\theta =0`$ and $`\theta =\pi /2`$ are compared in Fig. 2, and the difference grows with $`a/a_{\mathrm{o}sc}`$ until $`a/a_{\mathrm{o}sc}1`$, where it is $`5\%`$ of the correlation energy. At larger values of $`a/a_{\mathrm{o}sc}`$ the $`\theta `$ dependence decreases. This is an effect of the confining potential. The hyperradial wave function $`f(\rho )`$ varies on a scale $`a_{\mathrm{o}sc}`$. As shown above, the hyperangular wave function $`\mathrm{\Phi }_i`$ does, however, depend on the hyperradius through the combination $`\rho /a`$. Hence the couplings $`P`$ and $`Q`$ Eq. (7) are suppressed by factors $`a_{\mathrm{o}sc}/a`$ and $`(a_{\mathrm{o}sc}/a)^2`$ respectively when $`a>a_{\mathrm{o}sc}`$.

The $`\theta `$ dependence of the energy is shown in more detail in Fig. 3 for $`a>0`$. Physical conditions only determine $`\theta (\mathrm{mod}\pi )`$, but to exhibit the structure we have continued the energy curves over a larger range of $`\theta `$. Hence, given values of $`a/a_{\mathrm{o}sc}`$ and $`\theta `$, $`E(\theta )`$, $`E(\theta +\pi )`$, $`E(\theta +2\pi )`$ etc. along the same curve represent the energies of different physical states. States corresponding mainly to free atoms at large $`\rho `$ are represented by the nearly $`\theta `$-independent segments of the curves. In the adiabatic approximation, i.e. neglecting all couplings, these states have energies shown by the dotted lines. Energies of other states, i.e. states associated with molecular bonding, are given by the dashed lines. The difference between the solid and dashed/dotted lines is the effect of the couplings between the two hyperangular functions. The couplings are most important when $`aa_{\mathrm{o}sc}`$, and are almost entirely absent for $`0<a<0.1a_{\mathrm{o}sc}`$ and $`a>10a_{\mathrm{o}sc}`$. The results for $`a=10a_{\mathrm{o}sc}`$ indicate that for $`aa_{\mathrm{o}sc}`$ the couplings vanish and our solution becomes exact.

In the limit $`a/a_{\mathrm{o}sc}\pm \mathrm{}`$ Eq. (6) may be solved exactly. The energy tends to a constant value, which for the lowest state with a given $`\lambda _i(0)>0`$ is $`E=\mathrm{}\omega (1+\sqrt{\lambda _i(0)})`$. For the state which at large $`\rho `$ corresponds to free atoms, and for $`a>0`$, one has $`\lambda _i(0)19.94`$ and $`E5.47\mathrm{}\omega `$. The corresponding energy if only two-body correlations are included is $`E=6\mathrm{}\omega `$, and therefore the correlation energy is reduced by approximately $`18\%`$ by three-body effects. For $`\lambda _i(0)<0`$ the energies are functions of the boundary condition for $`\rho 0`$, i.e. of $`\theta `$, and are given by the solutions of the equation

$$\theta =\mathrm{arg}\left\{\frac{\mathrm{\Gamma }[1/2E/(2\mathrm{}\omega )s_0/2]}{\mathrm{\Gamma }[1s_0]}\right\}.$$
(10)
For $`E>0`$ this equation gives harmonic oscillator states shifted by the interaction, with spacings $`E_{n+1}E_n2\mathrm{}\omega `$, while for $`E<0`$ it gives bound three-body states, so-called Efimov states with energies scaling like $`E_{n+1}/E_n=\mathrm{exp}(2\pi /|s_0|)515`$  perturbed by the oscillator potential. For example $`\theta =0`$ gives energies $`E\mathrm{}`$, $`291649\mathrm{}\omega `$, $`566\mathrm{}\omega `$, $`0.85\mathrm{}\omega `$, $`2.27\mathrm{}\omega `$, $`4.48\mathrm{}\omega `$, $`\mathrm{}`$ .

Let us now examine the implications of our results for the properties of bulk matter. On physical grounds, we would expect the energy per particle in trapped few-body systems to be similar to that for bulk matter if the oscillator length is replaced by the typical particle separation. For $`|a|a_{\mathrm{osc}}`$, the energy per particle of two- and three-body states for which $`\lambda _i(0)>0`$ is $`\mathrm{}\omega `$ greater than in the absence of interactions. Since the density of particles is $`n1/a_{\mathrm{osc}}^3`$, this corresponds to an energy per particle $`(\mathrm{}^2/m)n^{2/3}`$. There is no dependence on $`\theta `$ for these states. The state corresponding to unbound atoms in a trap falls into this class for positive $`a`$. This suggests that the energy per particle of bulk matter has the universal dependence $`Cn^{2/3}`$, where $`C`$ is a constant that does not depend on the three-body parameter $`\theta `$. In particular, it does not have density-dependent contributions that are periodic in $`\mathrm{ln}n`$, contrary to the possibility suggested in Ref.. There are other states, such as those associated with molecular states bound by the interatomic potential, which do exhibit a strong dependence on $`\theta `$. These latter states are the final states in recombination processes, where diatomic molecules are formed from free atoms, and consequently recombination rates depend strongly on $`\theta `$.

The expansion Eq. (9) of the energy of three trapped bosons, may be compared to the low-density expansion for the energy of $`N`$ bosons in a homogenous gas

$$E=\frac{2\pi \mathrm{}^2an}{m}(N1)\left[1+\frac{128}{15}\left(\frac{na^3}{\pi }\right)^{1/2}+\mathrm{}\right].$$
(11)
For comparison, our result Eq. (9) for the energy of three trapped bosons with $`|a/a_{\mathrm{o}sc}|1`$, may be rewritten as

$$E=3\mathrm{}\omega +\frac{2\pi \mathrm{}^2a\overline{n}}{m}(N1)\left[1+\frac{c\sqrt{2\pi }}{6}\frac{a}{a_{\mathrm{o}sc}}+\mathrm{}\right].$$
(12)
Here the average density, weighted by the particle number, is defined by $`\overline{n}=n(𝐫)𝑑N(𝐫)/𝑑N(𝐫)=n^2(𝐫)𝑑𝐫/N`$, where $`dN(𝐫)=n(𝐫)d𝐫`$, and therefore for a homogeneous condensate $`\overline{n}=n`$. The three-body and many-body results agree for the term linear in $`a`$, which usually is the only term retained in, e.g. the Gross-Pitaevskii equation for condensates. Beyond the linear term the two expansions differ. The next term in Eq. (11) is of order $`(na^3)^{1/2}`$ relative to the first one, while in Eq. (12) it is proportional to $`a/a_{\mathrm{o}sc}`$. If $`Na_{\mathrm{o}sc}/a`$ the most important correction to the linear dependence of energy on $`a`$ will be the second term in Eq. (12), which arises from the trapping potential.

To determine when the three-body parameter $`\theta `$ becomes important, we note that, in the absence of the trap, the S-matrix element between the adiabatic state corresponding to free atoms and that corresponding to a diatomic molecule and an atom is proportional to $`(ka)^2`$, where $`k`$ is the wave number of the colliding atoms. Therefore the recombination rate is proportional to $`a^4`$, in agreement with earlier work . The constant of proportionality depends on $`\theta `$, because of the $`\theta `$ dependence of the final state. The same result is valid in a trap with $`a_{\mathrm{o}sc}|a|`$, where due to the zero point motion $`k`$ is of order $`1/a_{\mathrm{o}sc}`$. Therefore, the first term in the small $`a/a_{\mathrm{o}sc}`$ expansion of the energy \[Eq. (9)\] that depends explicitly on the coupling between different adiabatic states (and therefore on $`\theta `$) is of order $`(a/a_{\mathrm{o}sc})^4`$. In the expansion Eq. (12) this corresponds to the cubic term, since there one power of $`a/a_{\mathrm{o}sc}`$ has been factored out. To this order the result agrees with the results for the energy density of the homogeneous Bose gas in Ref. . As we saw above, for $`aa_{\mathrm{o}sc}`$ the couplings $`P`$ and $`Q`$ between hyperangular functions are reduced, which will tend to suppress recombination rates.

In conclusion, we have found that there are states of the trapped three-boson system with very weak dependence of the energies on the three-body parameter $`\theta `$. This dependence vanishes completely in the limit of an infinitely large scattering length. The fact that the energy of these states is independent of the three-body parameter in the limit of large scattering lengths gives some grounds to hope that the energy of trapped systems with a larger number of bodies is insensitive to the parameters describing the four- and higher-body problems, and that similar conclusions apply to the energy of bulk matter. Topics for future study include higher-body correlations, and the lifetime of a condensate with a large scattering length.

We are grateful to Eric Braaten and Oleg I. Kartavtsev for valuable discussions.


text:
# Statistical properties of stock order books: empirical results and models

## Abstract

We investigate several statistical properties of the order book of three liquid stocks of the Paris Bourse. The results are to a large degree independent of the stock studied. The most interesting features concern (i) the statistics of incoming limit order prices, which follows a power-law around the current price with a diverging mean; and (ii) the shape of the average order book, which can be quantitatively reproduced using a ‘zero intelligence’ numerical model, and qualitatively predicted using a simple approximation.

Commissariat à l’Energie Atomique, Orme des Merisiers
91191 Gif-sur-Yvette cedex, France

Science & Finance, CFM, 109-111 rue Victor Hugo
92 353 Levallois cedex, France

<sup>∗∗</sup> Laboratoire de Physique Théorique et Modèles Statistiques
Université Paris Sud, Bat. 100, 91 405 Orsay cedex, France

Financial markets offer an amazing source of detailed data on the collective behaviour of interacting agents. It is possible to find many reproducible patterns and even to perform experiments, which bring this atypical subject into the realm of experimental science. The situation is simple and well defined, since many agents, with all the same goal, trade the very same asset. As such, the statistical analysis of financial markets also offers an interesting testing ground not only for economic theories, but also for more ambitious theories of human activities. One may indeed wonder to what extent it is necessary to invoke human intelligence or rationality to explain the various universal statistical laws which have been recently unveiled by the systematic analysis of very large data sets.

Many statistical properties of financial markets have already been explored, and have revealed striking similarities between very different markets (different traded assets, different geographical zones, different epochs) . For example, the distribution of price changes exhibits a power-law tail with an apparently universal exponent . The volatility of most assets shows random variations in time, with a correlation function which decays as a small power of the time lag, again quite universal across different markets .

Here, we study the statistics of the order book, which is the ultimate ‘microscopic’ level of description of financial markets. The order book is the list of all buy and sell ‘limit’ orders, with their corresponding price and volume, at a given instant of time. A limit order specifies the maximum (resp. minimum) price at which an investor is willing to buy (resp. sell) a certain number of shares (volume). At a given instant of time, all limit buy orders are below the best buy order called the bid price, while all sell orders are above the best sell order called the ask price. When a new order appears (say a buy order), it either adds to the book if it is below the ask price, or generates a trade at the ask if it is above (or equal to) the ask price (we call all these ‘market orders’ even if technically they could also be ‘marketable limit orders’). The price dynamics is therefore the result of the interplay between the order book and the order flow. The study of the order book is very interesting both for academic and practical reasons. It provides intimate information on the processes of trading and price formation, and reveals a highly non trivial structure of the agents’ expectations (see below): as such, it is of importance to test some basic notions of economics and models of market microstructure. The practical motivations are also clear. Regulatory bodies want to set the rules of exchange markets such as to produce fair, orderly trading that maximize flow. For market participants, issues such as the market impact or the relative merit of limit versus market orders are determined by the structure and dynamics of the order book.

Complete data on the order book of certain markets, such as the Paris Bourse (now Euronext), is now available, and contains a particularly abundant information. All orders (representing hundreds of thousands orders per month on liquid stocks) on all stocks are stored, which makes possible the reconstitution of the full order book at any instant of time, including orders that are not directly observable on traders’ screens.

Many questions can be studied; the systematic investigation of these data sets is only very recent  and has motivated a number of interesting theoretical work . Here, we mostly focus on ‘static’ properties of the order book, such as the distribution of incoming limit orders, the average shape of the order book in the moving reference frame of the price, or the distribution of volume at the bid/ask. Many other ‘dynamical’ properties can also be analyzed, such as the response of the price to order flow, the full temporal correlation of the book, etc. In order to keep our study reasonably focused, we leave these dynamical aspects, that we have not fully understood yet, for a subsequent publication.

Our main results are as follows: (a) the price at which new limit orders are placed is, somewhat surprisingly, very broadly (power-law) distributed around the current bid/ask; (b) the average order book has a maximum away from the current bid/ask, and a tail reflecting the statistics of the incoming orders; (c) the distribution of volume at the bid (or ask) follows a Gamma distribution. Notably, we find exactly the same statistical features, for the three liquid stocks studied. We then study numerically a ‘zero intelligence’ model of order book which reproduces most of these empirical results. Finally, we show how the characteristic shape of the average order book can be analytically predicted, using a simple approximation.

The data provided by Paris Bourse gives the history of all transactions, with their price, volume and time stamp, of all quotes (bid and ask prices) with the corresponding volumes, and of all orders, with their price, volume and time stamp. This information allows in principle to reconstruct the whole order book at any instant of time. However, some limit orders (roughly $`10\%`$ of them) are modified or cancelled before being executed; unfortunately the data base only contains this qualitative information, but not the time at which a given order is modified/cancelled, nor the resulting new price. When reconstructing the order book at a given instant of time, we have therefore made two extreme assumptions, which lead to nearly identical conclusions. Either we discard these orders altogether, or we keep them until the time where we can be sure that they have been previously modified/cancelled, otherwise they would have been executed since the transaction price was observed to be below (for buy orders) or above (for sell orders) the corresponding limit price. The results given below uses the second procedure. We had available the data corresponding to all stocks during February 2001, from which we have extracted three of the most liquid stocks: France-Telecom (F.T.), Vivendi and Total. Some basic information is summarized in Table I for France-Telecom and Total. We expect that our findings will not depend on the particular month considered, but that some differences may appear when one studies stocks with smaller capitalisation (for example the average bid-ask spread is much larger): these questions will be investigated in the near future. Let us call $`a(t)`$ the ask price at time $`t`$ and $`b(t)`$ the level of bid price at time $`t`$. The midpoint $`m(t)`$ is the average between the bid and the ask: $`m(t)=[a(t)+b(t)]/2`$. The smallest possible change of these quantities, called the ‘tick’, was $`0.05`$ Euros for France-Telecom and Vivendi, and $`0.10`$ Euros for Total.

We denote by $`b(t)\mathrm{\Delta }`$ the price of a new buy limit order, and $`a(t)+\mathrm{\Delta }`$ the price of a new sell limit order. Notice that $`\mathrm{\Delta }`$ can be negative (this is the only case where the spread $`g(t)=a(t)b(t)`$ can be narrowed), but is always larger than $`g(t)`$ (otherwise it would be a market order). A first interesting question concerns the distribution density of $`\mathrm{\Delta }`$, i.e. the distance between the current price and the incoming limit order. We find that $`P(\mathrm{\Delta })`$ is identical for buy and sell orders (up to statistical fluctuations); the shape of $`P(\mathrm{\Delta })`$ is found to be very well fitted (see Fig. 1) by a power-law:

$$P(\mathrm{\Delta })\frac{\mathrm{\Delta }_0^\mu }{(\mathrm{\Delta }_1+\mathrm{\Delta })^{1+\mu }},\mathrm{\Delta }1$$
(1)
with an exponent $`\mu 0.6`$ for all three stocks. This power-law extends from $`1`$ tick to over $`100`$ ticks (sometimes even $`1000`$ ticks), corresponding to a relative change of price of $`5\%`$ to $`50\%`$ . There are also orders placed at the bid (or ask) or within the spread. One finds that $`P(\mathrm{\Delta }=1)P(\mathrm{\Delta }=0)P(\mathrm{\Delta }=1)`$, and these orders add up to roughly half of the total number of orders. The distribution of orders Eq. (1) has a maximum around the current price (this was already noted in ). But the fact that $`\mu <1`$ means that the average $`\mathrm{\Delta }`$ is formally infinite (although of course the distribution is ultimately cut-off for large $`\mathrm{\Delta }`$’s).

It is quite surprising to observe such a broad distribution of limit order prices, which tells us that the opinion of market participants about the price of the stock in a near future could be anything from its present value to $`50\%`$ above or below this value, with all intermediate possibilities. This means that market participants believe that large jumps in the price of stocks are always possible, and place orders very far from the current price in order to take advantage of these large potential fluctuations. (Note that placing a limit order is free of charge.) A naive argument would then suggest that the probability to place an order at distance $`\mathrm{\Delta }`$ should be proportional to the probability that the price moves more than $`\mathrm{\Delta }`$ in order to meet the order. Since the tail of the distribution of price increments is a power-law with an exponent $`\mu _{\delta p}3`$ , this would indeed lead to a power-law for $`P(\mathrm{\Delta })`$, but with a value $`\mu =\mu _{\delta p}12`$ larger than the observed one. This however does not take into account the fact that market participants have different time horizons – large $`\mathrm{\Delta }`$’s presumably correspond to more patient investors. The value of $`\mu `$ should result from an interplay between the perceived distribution of future price changes and the distribution of time horizons, which could well be itself a power-law. We feel that more theoretical work is needed to fully account for this striking result.

Limit orders strongly vary in volume. We find that the unconditional limit order size, $`\varphi `$, is distributed uniformly in log-size, between $`10`$ and $`50,000`$ (both for buy or sell orders). One can study the correlation between the incoming volume $`\varphi `$ (number of shares for a given order) and the distance $`\mathrm{\Delta }`$ between the order and the current price. We find that the conditional averaged volume $`\varphi |_\mathrm{\Delta }`$ is roughly independent of $`\mathrm{\Delta }`$ between $`1`$ and $`20`$ ticks, but decays as a power law $`\mathrm{\Delta }^\nu `$, with $`\nu 1.5`$ beyond $`\mathrm{\Delta }^{}`$ (with $`\mathrm{\Delta }^{}20`$ ticks for FT, and $`50`$ ticks for the other two stocks). Not unexpectedly, extremely far limit orders tend to be of smaller volumes. Since the average volume arriving at $`\mathrm{\Delta }`$ is equal to the number of orders at $`\mathrm{\Delta }`$ multiplied by $`\varphi |_\mathrm{\Delta }`$, one predicts, using the shape of $`P(\mathrm{\Delta })`$, that the distribution of incoming volume should then decay as $`\mathrm{\Delta }^{1\mu }`$ for $`\mathrm{\Delta }<\mathrm{\Delta }^{}`$ and as $`\mathrm{\Delta }^{1\mu \nu }`$ for large $`\mathrm{\Delta }>\mathrm{\Delta }^{}`$, a feature that we have confirmed directly.

We now turn to the shape of the order book. The order flow is maximum around the current price, but an order very close to the current price has a larger probability of being executed and thus disappear from the book. It is thus not a priori clear what will be the average shape of the order book. Quite interestingly, we find that the (time-averaged) volume of the queue in the order book (that we will call for simplicity the ‘average order book’) is symmetrical, and has a maximum away from the current bid (ask): see Fig. 2. This was also noted in . This shape furthermore appears to be universal, up to a rescaling of both the $`\mathrm{\Delta }`$ axis and the volume axis, at least for the three stocks studied. In order to test the universality of this result, we are currently studying less liquid stocks, and in a different time interval. The empirical determination of the average order book is the central result of this paper, and simple models that explain this shape will be discussed below.

The next question concerns the volume fluctuations around this average shape. We first study the distribution $`R(V)`$ of volume at the bid (or ask). Again, the two are identical, and can be fitted by a Gamma distribution for the volume (Fig. 3):

$$R(V)V^{\gamma 1}\mathrm{exp}\left(\frac{V}{V_0}\right).$$
(2)
We find $`\gamma 0.70.8`$ for all three stocks. A Gamma distribution with $`\gamma 1`$ has its maximum for $`V=0`$. This shows that the most probable value of the volume at bid is very small, although its typical value is quite large ($`V_02700`$ for France-Telecom).

Therefore, the order book has strong fluctuations and at a given instant of time can look rather different from its average, Fig. 2. We have also studied the fluctuations of volume in the book as a function of $`\mathrm{\Delta }`$, defined as:

$$\sigma _V(\mathrm{\Delta })=\sqrt{V^2|_\mathrm{\Delta }V|_\mathrm{\Delta }^2}.$$
(3)
We find that $`\sigma _V(\mathrm{\Delta })`$ is of order $`1`$ for $`\mathrm{\Delta }=1`$ (as expected from the above Gamma distribution), and is roughly constant up to $`\mathrm{\Delta }=50`$. This shows that far from the current price, the instantaneous order book shows large relative fluctuations, which is indeed reasonable. We have also studied the full covariance matrix of the fluctuations $`C_V(\mathrm{\Delta },\mathrm{\Delta }^{})`$ (note that $`C_V(\mathrm{\Delta },\mathrm{\Delta })\sigma _V^2(\mathrm{\Delta })`$). For France-Telecom, its first eigenvalue corresponds to a dilation of the book, whereas the second one reflects an even/odd oscillation (orders are more numerous at prices in tenth of Euros than in twentieth of Euros). Higher eigenvalues are structureless, suggesting that the orderbook has no ‘elasticity’. This is expected, since most investors do not have access to the order book in real time.

We now turn to a quantitative interpretation of the data. First, we have simulated an artificial ‘zero intelligence’ order book, much in the same spirit as in , but more directly inspired from our empirical results. Limit orders, with a size distributed uniformly on a log scale, are launched at random with an a priori ‘sprinkling’ distribution $`P(\mathrm{\Delta })`$, which we take identical to the empirical one. Here we depart from , where all orders have unit size and the ‘sprinkling’ distribution is uniform between $`0`$ and a certain $`\mathrm{\Delta }_{\mathrm{max}}`$. The size distribution was actually found to play a minor role except in the far tail of the average order book, where the exponent $`\nu `$ introduced above becomes important. For example, we also studied the case of orders of equal size with very similar results.

On the other hand, choosing the correct ‘rain’ distribution $`P(\mathrm{\Delta })`$ is crucial to reproduce the shape of the order book. We also launch a certain fraction $`p_m`$ of market orders (chosen to be typically $`1/61/4`$th of the total) in order to trigger trades. The results we find (in particular the averge size of the bid-ask spread) are quite sensitive to the value of $`p_m`$. Finally, with a constant probability $`\mathrm{\Gamma }`$ (typically of order $`10^3`$) per unit time, independently of both size and position in the book, an order is cancelled. We have found that this simple model is able to reproduce quantitatively many of the features observed in the empirical data, such a the shape of the order book (see Fig. 2), or the Gamma distribution of volume at the bid/ask $`R(V)`$, with a similar value for the exponent $`\gamma `$. This model also reproduces other empirical properties, such as the short time dynamics of the ‘midpoint’ and the sublinear volume dependence of the response to an incoming market order . However, we note that the above simple model overestimates the average size of the queue close to $`\mathrm{\Delta }=0`$. This might be due either to the fact that our procedure to treat modified orders leads to an systematic bias which is large near the best price, or that orders close to the best price have a non negligible strategic content that our model fails to capture.

Finally, we discuss a simple analytical approximation which allows us to compute the average order book from the ingredients of the numerical model. This was also attempted in . Although very few details were given in , our method appears to be quite different from theirs. Consider sell orders. Those at distance $`\mathrm{\Delta }`$ from the current ask at time $`t`$ are those which were placed there at a time $`t^{}<t`$, and have survived until time $`t`$, that is, (i) have not been cancelled; (ii) have not been touched by the price at any intermediate time $`t^{\prime \prime }`$ between $`t^{}`$ and $`t`$. An order at distance $`\mathrm{\Delta }`$ at time $`t`$ in the reference frame of the ask $`a(t)`$ appeared in the order book at time $`t^{}`$ at a distance $`\mathrm{\Delta }+a(t)a(t^{})`$. The average order book can thus be written, in the long time limit, as:

$$\rho (\mathrm{\Delta },t)=_{\mathrm{}}^tdt^{}𝑑uP\left(\mathrm{\Delta }+u\right)𝒫\left(u|𝒞(t,t^{})\right)\mathrm{e}^{\mathrm{\Gamma }(tt^{})},$$
(4)
where $`𝒫\left(u|𝒞(t,t^{})\right)`$ is the conditional probability that the time evolution of the price produces a given value of the ask difference $`u=a(t)a(t^{})`$, given the condition that the path always satisfies $`\mathrm{\Delta }+a(t)a(t^{\prime \prime })0`$ at all intermediate times $`t^{\prime \prime }[t^{},t]`$. The evaluation of $`𝒫`$ requires the knowledge of the statistics of the price process. Because of the exponential cut-off $`\mathrm{e}^{\mathrm{\Gamma }(tt^{})}`$, only the short time behaviour of the process is relevant where the confinement effects of the order book on the price are particularly important . Nevertheless, to make progress, we will assume that the process is purely diffusive. In this case, $`𝒫`$ can be calculated using the method of images. One finds:

$$𝒫\left(u|𝒞(t,t^{})\right)=\frac{1}{\sqrt{2\pi D\tau }}\left[\mathrm{exp}\left(\frac{u^2}{2D\tau }\right)\mathrm{exp}\left(\frac{(2\mathrm{\Delta }+u)^2}{2D\tau }\right)\right],$$
(5)
where $`\tau =tt^{}`$ and $`D`$ is the diffusion constant of the price process.

After a simple computation, one finally finds, up to a multiplicative constant which only affects the overall normalisation of $`\rho _{\mathrm{st}}(\mathrm{\Delta })=\rho (\mathrm{\Delta },t\mathrm{})`$:

$$\rho _{\mathrm{st}}(\mathrm{\Delta })=\mathrm{e}^{\alpha \mathrm{\Delta }}_0^\mathrm{\Delta }duP(u)\mathrm{sinh}(\alpha u)+\mathrm{sinh}(\alpha \mathrm{\Delta })_\mathrm{\Delta }^{\mathrm{}}duP(u)\mathrm{e}^{\alpha u},$$
(6)
where $`\alpha ^1=\sqrt{D/2\mathrm{\Gamma }}`$ measures the typical variation of price during the lifetime of an order, and fixes the scale over which the order book varies. When the distribution of order flows $`P(\mathrm{\Delta })`$ has the power-law shape Eq. (1) with $`\mu <1`$, the parameter $`\alpha `$ can be rescaled away in the ‘continuous’ limit where $`\alpha ^1`$ is much larger than the tick size (which is the relevant limit for stocks, where $`\alpha ^110`$). In this case, the shape of the average order book only depends on $`\mu `$ and $`\widehat{\mathrm{\Delta }}=\alpha \mathrm{\Delta }`$, and is given by the following convergent integral:

$$\rho _{\mathrm{st}}(\widehat{\mathrm{\Delta }})=\mathrm{e}^{\widehat{\mathrm{\Delta }}}_0^{\widehat{\mathrm{\Delta }}}duu^{1\mu }\mathrm{sinh}(u)+\mathrm{sinh}(\widehat{\mathrm{\Delta }})_{\widehat{\mathrm{\Delta }}}^{\mathrm{}}duu^{1\mu }\mathrm{e}^u.$$
(7)
For $`\mathrm{\Delta }0`$, the average order book vanishes in a singular way, as $`\rho _{\mathrm{st}}(\mathrm{\Delta })\mathrm{\Delta }^{1\mu }`$, whereas for $`\mathrm{\Delta }\mathrm{}`$, the average order book reflects the incoming flow of orders: $`\rho _{\mathrm{st}}(\mathrm{\Delta })\mathrm{\Delta }^{1\mu }`$. We have plotted in Fig.4 the average order book obtained numerically from the ‘zero intelligence’ model (which is very close to empirical data, see Fig. 2) and compared it with Eq. (7), with $`\mu =0.6`$ and various choices of parameters. After rescaling the two axes, the various numerical models lead to very similar average order books, and the analytic approximation appears rather satisfactory, bearing in mind the roughness of our ‘diffusive’ approximation.

The shape of the average order book therefore reflects the competition between a power-law flow of limit orders with a finite lifetime, and the price dynamics that removes the orders close to the current price. These effects lead to a universal shape which will presumably hold for many different markets, provided the lifetime of orders is sufficiently long compared to the typical time between trades, and the volatility on the scale of the orders lifetime is somewhat larger than the ticksize. <sup>1</sup><sup>1</sup>1Preliminary results indeed show that the average order book in futures markets has a similar hump away from the current midpoint. On the other hand, the detailed shape of the book in the immediate vicinity of the best price (i.e., around $`\mathrm{\Delta }=0`$) is not expected to be universal.

In conclusion, we have investigated several ‘static’ properties of the order book. We have found that the results appear to be independent of the stock studied. The most interesting features concern (i) the statistics of incoming limit order prices, which follows a power-law around the current price with a diverging mean – suggesting, quite surprisingly, that market participants believe that very large variations of the price are possible within a rather short time horizon; and (ii) the shape of the average order book, which can be quantitatively reproduced using a ‘zero intelligence’ numerical model, and qualitatively predicted using a simple approximation. One of the most interesting open problems is, in our view, to explain the clear power-law behaviour of the incoming orders and the difference between our results and those of . The dynamical properties of the order book, and the way the price reacts to order flow, will be the subject of a further study .

Acknowledgements: We thank Jean-Pierre Aguilar, Jelle Boersma, Damien Challet, J. Doyne Farmer, Laurent Laloux, Andrew Matacz, Philip Seager and Denis Ullmo for stimulating and useful discussions.