text:
# Disentangling discrepancies between stellar evolution theory and sub-solar mass stars

## 1 Introduction

It is of the utmost importance for stellar evolution theory to match at least the best known objects (the Sun and non-interacting binary systems) before any attempt to derive proper information for star clusters or stellar populations in galaxies is made.

Some binary stars are known to provide uncomparable astrophysical laboratories to calibrate theoretical stellar evolutionary models. Since the individual components of well-detached binary systems can be assumed to be two single-like stars with a common origin, they share the same chemical composition and same age, and, therefore, the observed parameters of both components are expected to be matched by a single isochrone at the same chemical composition.

Whilst in general a good overall agreement is found, some studies (pioneered by Popper (1997), see also e.g. Pols et al. (1997), Lastennet et al. (1999), Clausen et al. (1999) and Lastennet & Valls-Gabaud (2002) for a systematic analysis) have unambiguously found that systems with components in the 0.7 to 1.1 $`M_{}`$ mass range raise serious difficulties.

There are many well-studied systems with components in this mass range which present puzzling discrepancies with the predictions of stellar evolution theory. For example, if the seismic observations of $`\alpha `$ Cen A are fitted to state-of-the-art structure models, the resulting masses differ from the ones derived from dynamical analyses (Pourbaix et al. 2002, Thévenin et al. 2002). The visual binary 85 Peg, some eclipsing binaries in the field, and in the Hyades open cluster are other clear examples (see Lastennet et al. (2002) for a brief review of the results obtained so far on these objects).

This situation is very puzzling because one would expect these stars to be reasonably well understood. Indeed, stars with masses larger than $``$0.6 M bypass difficulties in the treatment of the equation of state and the atmosphere, while stars with masses larger than about 1.1 M have a permanent convective core, introducing an additional parameter, the amount of overshooting, for their modelling. Hence, current stellar evolutionary models are expected to be able to match the basic properties of stars in the 0.7 – 1.1 M mass range, provided a good description of the convection in their envelopes is used.

One of these puzzling objects is the eclipsing binary UV Piscium (HD 7700, HIP 5980, hereafter UV Psc), a detached main-sequence binary (with components of types G and K), which is also a short period ($`<`$1 day) RS CVn111RS CVn are binary systems where one of the components presents an extreme form of solar-like activity with starspots and variable magnetic fields. See Montesinos et al. (1988) for a discussion. . This binary deserves particular attention ever since Popper (1997) derived age estimates for each component which differed by a factor larger than two, with an absolute age difference larger than 11 Gyr. Montesinos et al. (1988) classified its components as “apparently normal main-sequence stars”. What is the origin of this discrepancy?

Popper (1997) used the Geneva models (Schaller et al. 1992) for $`Z=`$0.02 and $`Y=`$0.30, i.e. with a helium abundance to metallicity ratio fixed to $`\mathrm{\Delta }`$$`Y`$/$`\mathrm{\Delta }`$$`Z`$$`=`$3, and $`\alpha _{\mathrm{MLT}}`$ fixed to a solar-calibrated value ($`\alpha _{\mathrm{MLT}}`$$`=`$1.63).

Since the RS Canum Venaticorum stars show activity, mass transfer episodes between active regions of both components are possible. In this case, the application of the single-star theoretical models used by Popper may not be relevant since the evolution of both components would not be independent any longer. However, for UV Psc, mass transfer does not account for the age difference because the computation of the Roche lobe radii (e.g. with the formulae given by Eggleton 1983) shows that none of these stars have overflown their Roche lobes. Therefore, we consider that UV Psc is a well-detached system and assume that the evolution of each component is independent.

Another concern about RS CVn stars is the possible distortion of light curves due to starspots, which could bias the derivation of their fundamental stellar parameters. Jassur & Kermani (1994) suggested the presence of cool spots on the surface of the primary component, but according to Popper (1997, and references therein), this should not affect the value of the quantities used in this paper. The period has been recently revised by Sowell et al. (2001) who didn’t find any evidence for a cyclic modulation or period change. Hence activity cannot account for the age discrepancy either.

To tackle this problem, we explore the influence of each physical parameter. A stellar evolutionary code allowing to compute tracks for different and independent physical parameters is needed for such a purpose, and so we choose the CESAM code (Morel 1997) which fulfills this condition and which has already been successfully tested in a variety of astrophysical contexts (e.g. Zahn et al. 1997, Suran et al. 2001, Cordier et al. 2002).

In this paper, we present a detailed study of the system UV Psc with the CESAM models. Our goal is twofold: 1) to check if there exists a solution giving consistent ages for both components for a set of ($`Z`$, $`Y`$), and 2) to constrain $`\alpha _{\mathrm{MLT}}`$ for both components and check if these values are consistent with the generally assumed solar value, $`\alpha _{\mathrm{MLT}}`$.

The paper is organised as follows: Sect. 2 deals with the description of the CESAM models used and the grid computed for the purpose of this paper. Sect. 3 presents the results which fit both components of the eclipsing binary UV Psc, with particular emphasis on the constraints derived from the $`T_{\mathrm{eff}}`$-radius diagram on the $`\alpha _{\mathrm{MLT}}`$ parameters. We also discuss the robustness of the results against changes in the $`T_{\mathrm{eff}}`$ scale and in diffusion. Finally, Sect. 4 provides a summary of our results.

## 2 Evolutionary stellar models

The stellar evolution calculations were computed with the CESAM code (Version 3, Morel morel97 (1997)), running at the Coimbra Observatory. Details on the physics of these models can be found in Lebreton et al. (lebreton99 (1999)). Here we just provide a short summary of the main CESAM physical inputs: the CEFF equation of state is used, including Coulomb corrections to the pressure (Eggleton et al. eggleton73 (1973), Christensen-Dalsgaard jcd91 (1991)); the nuclear reactions rates are from Caughlan and Fowler (caughlan88 (1988)); the solar mixture is from Grevesse and Noels (grevesse93 (1993)); the OPAL opacities (Iglesias and Rogers iglesias96 (1996)) are used and complemented at low temperatures by opacity data from Alexander and Ferguson (alexander93 (1993)) following a prescription by Houdek & Rogl (houdek96 (1996)); the atmosphere is described with an Eddington $`T(\tau )`$-law; the convection is treated according to the mixing-length theory from Böhm-Vitense (bohm-vitense58 (1958)), with the formalism of Henyey et al. (1965), giving the mixing-length scale ($`\alpha _{\mathrm{MLT}}`$$`\times H_p`$) in terms of a free parameter times the local pressure scale height $`H_p`$.

With these physical ingredients, we computed a solar model that fits the observed luminosity and radius with an accuracy of $`10^4`$, with $`\alpha _{\mathrm{MLT}}`$$`=1.63`$, helium abundance $`Y_{}`$$`=0.268`$ and $`Z_{}`$$`=0.0175`$ for the commonly accepted solar age of $``$4.6 Gyr (Dziembowski et al. dziembowski99 (1999)) and solar abundances as derived by Grevesse & Noels (grevesse93 (1993)).

Given this calibration, we computed a large grid of CESAM models from the Zero Age Main Sequence (ZAMS) to the terminal Age Main Sequence (TAMS, central hydrogen exhaustion), or for ages $`t`$ younger than 10 Gyr222UV Psc being a binary belonging to the close solar neighbourhood ($`d`$ 63 pc), ages older than the Galactic disc (likely to be 8-10 Gyr old, at most, see Carraro 2000) should be irrelevant here.. This was done for each component of UV Psc by fixing the masses to their accurately measured values: $`M_A=`$ 0.975 $`M_{}`$ and $`M_B=`$ 0.760 $`M_{}`$ (Popper 1997). This approximation is justified since the relative mass accuracy is better than 1% in each star.
We considered different values of the metallicity $`Z`$ (0.01, 0.02, 0.03, 0.04 and 0.05), helium content $`Y`$ (0.25, 0.28, 0.31 and 0.34) and mixing length parameter $`\alpha _{\mathrm{MLT}}`$ (0.8, 1.3 and 1.7) in a large three-dimensional grid, since we do not want to introduce correlations between these parameters. We also computed additional models for the secondary component with $`\alpha _{\mathrm{MLT}}`$$`=`$0.1 for $`Z=`$ 0.01 and $`Z=`$ 0.02333 The internal structure of such low-$`\alpha _{\mathrm{MLT}}`$ models does not present any convective region at all (this can be explained by the fact that the convective efficiency is proportional to the square of the mixing length and hence to $`\alpha _{\mathrm{MLT}}^2`$). The discussion on how realistic these models are is beyond the scope of this paper, but we would like to point out that we only used them to interpolate models with $`\alpha _{\mathrm{MLT}}`$ slightly lower than 0.8. We stress that the quantities $`Z`$ and $`Y`$ are independent, unlike most of the extant publicly available theoretical models where $`Y`$ is derived from a fixed law once $`Z`$ is determined. All the models computed with the input physics described above are displayed on Fig. 1.

## 3 Analysis of the results

UV Psc is an eclipsing binary with very accurate masses and radii derived for both components. The fundamental data used in this paper (from Popper, 1997) are listed in Table 1 and are a key point of this study: the measures by Popper represent such a significant revision and improvement that a comparison with older data is irrelevant.

For instance, Montesinos et al. (1988) provided a comprehensive work on RS CVn EBs and derived a metal-rich solution for UV Psc ($`Z=`$0.04 and $`Y=`$0.25) from the log g-$`T_{\mathrm{eff}}`$ diagram, but the data they used were very significantly revised by Popper ($`\mathrm{\Delta }`$$`M_A`$$`=`$0.3 $`M_{}`$, etc) hampering any useful comparison.

First, before exploring all possible solutions in the ($`Z`$, $`Y`$, $`\alpha _{\mathrm{MLT}}`$, $`t`$) parameter space, interesting qualitative conclusions can be derived from a simple inspection of Fig. 1. Whatever the values of $`Y`$ and $`\alpha _{\mathrm{MLT}}`$ are, all models with $`Z`$$``$0.04 are unable to match simultaneously the location of both components, strongly suggesting that metal-rich solutions can be ruled out. In addition, while a fit to both components seems to be possible in some panels –e.g., the upper left panel at $`(Z,Y)=`$(0.01,0.25)– inconsistent ages appear, excluding these possible solutions. This is also the case for the solution found by Popper at $`(Z,Y)=`$(0.02,0.30) using the Geneva models.

Altogether, this means that we expect to find solutions with consistent ages for chemical compositions in the range $`Z=`$0.01-0.02 and $`Y=`$0.28-0.34. At this level of qualitative discussion, it appears difficult to infer strong conclusions about the $`\alpha _{\mathrm{MLT}}`$ parameter but the panels of Fig. 1 suggest that $`\alpha _{\mathrm{MLT},\mathrm{B}}`$ may be systematically smaller that $`\alpha _{\mathrm{MLT},\mathrm{A}}`$.

To obtain more detailed conclusions, a systematic work matching some selected relevant criteria has to be performed. For this purpose, we explore by interpolation444$`Z`$, $`Y`$ and $`\alpha _{\mathrm{MLT}}`$ in steps of 0.001, 0.005 and 0.05, respectively. the ($`Z`$, $`Y`$, $`\alpha _{\mathrm{MLT}}`$, $`t`$) parameter space described in the previous section to find out the best mathematical solutions matching the location of UV Psc in the $`T_{\mathrm{eff}}`$-radius diagram for both components ($`i=`$ A, B). This is performed by minimizing the $`\chi ^2`$-functional

$`\chi ^2(Z_i,Y_i,\alpha _{\mathrm{MLT}i},t_i)`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{log}T_{\mathrm{eff}}(track)\mathrm{log}T_{\mathrm{eff}i}}{\sigma (\mathrm{log}T_{\mathrm{eff}i})}}\right)^2`$ (1)
$`+\left({\displaystyle \frac{\mathrm{log}R(track)\mathrm{log}R_i}{\sigma (\mathrm{log}R_i)}}\right)^2`$
We assume that both components share the same chemical composition, i.e. $`Z_A=Z_B`$ and $`Y_A=Y_B`$. Moreover, they should also have equal (or at least consistent) ages, hence $`t_At_B`$ should be minimal.

In Fig. 2, we show all the solutions in three projections: $`Z`$-$`Y`$, $`Z`$$`\alpha _{\mathrm{MLT},\mathrm{A}}`$ and $`Z`$$`\alpha _{\mathrm{MLT},\mathrm{B}}`$ diagrams. The solutions (small dots on Fig. 2) which have a good $`\chi ^2`$ for both components ($`\chi `$$`{}_{i}{}^{}{}_{}{}^{2}`$ $``$ $`\chi `$$`{}_{i,min}{}^{}{}_{}{}^{2}`$ $`+`$ 1 with $`i=`$A, B) show that i) high metallicities ($`Z>`$0.032) are indeed ruled out, ii) 1.7 $``$$`\alpha _{\mathrm{MLT},\mathrm{A}}`$ $``$ 0.9, and iii) 1.0 $``$$`\alpha _{\mathrm{MLT},\mathrm{B}}`$ $``$ 0.6.

Adding the constraint that the ages are consistent with each other ($`\mathrm{\Delta }`$t/t$``$20%) provides the subsample shown with crosses. It defines a much smaller region (Fig. 2) and shows that i) the metallicity is solar or sub-solar ($`Z<`$0.019) and the helium content is larger than the solar value ($`Y>`$0.28), ii) 1.1 $``$ $`\alpha _{\mathrm{MLT},\mathrm{A}}`$ $``$ 0.9, and iii) 0.7 $``$ $`\alpha _{\mathrm{MLT},\mathrm{B}}`$ $``$ 0.6.

The best solution (filled circle in Fig. 2) is obtained with a $`Z=`$0.012, $`Y=`$0.31 chemical composition and provides an excellent agreement in age (0.9% discrepancy). We predict an age of about 1.9 Gyr for both components of UV Psc A (1.92$`{}_{0.4}{}^{}{}_{}{}^{+0.6}`$ Gyr) and B (1.94$`{}_{0.4}{}^{}{}_{}{}^{+3.0}`$ Gyr). While a unique ($`Z`$, $`Y`$, $`t`$) combination is needed to simultaneously match the (masses, radii, $`T_{\mathrm{eff}}`$s) of both components, this solution implies a different value for the convection parameter for the two components of the binary: $`\alpha _{\mathrm{MLT}}`$(UV Psc A) $`=`$ 0.95$`\pm `$0.12 and $`\alpha _{\mathrm{MLT}}`$(UV Psc B) $`=`$ 0.65$`\pm `$0.07. The fit provided by this solution is shown in Fig. 3.

The non-solar values for $`\alpha _{\mathrm{MLT}}`$ parameters that we derive obviously deserve further comments, but the robustness of these solutions needs to be discussed first.

### 3.1 Influence of the $`T_{\mathrm{eff}}`$ scaling

Since our main concern was to solve the age discrepancy found by Popper (1997), we used the same values for the mass, radius and $`T_{\mathrm{eff}}`$. We have shown in the previous section that this discrepancy can be solved within the same data set. However, whilst the individual masses and radii are very accurately defined, the $`T_{\mathrm{eff}}`$s are not direct quantities, and depend on colour indices and photometric calibrations. We would like to point out that according Popper (popper1998 (1998)) the presence of dark spots has relatively little effect on the color index of a star. To test the possible influence of the $`T_{\mathrm{eff}}`$ scaling on our results, we used as an alternative the metallicity-dependent empirical calibrations of Alonso et al. (alonso96 (1996)), $`T_{\mathrm{eff}}`$ $`=`$f($`BV`$, \[Fe/H\]) and used the $`BV`$ indices of both components from Table 6 of Popper (1997) assuming that the reddening can be neglected because UV Psc is in the inner solar neighbourhood ($``$ 63 pc according to its Hipparcos parallax). Since there is no direct determination of the iron abundance \[Fe/H\] available for UV Psc, we have computed both $`T_{\mathrm{eff}}`$s from the Alonso et al.’s formula (1) assuming three different metallicities: $`T_{\mathrm{eff},\mathrm{A}}`$$`=`$ 5574 K and $`T_{\mathrm{eff},\mathrm{B}}`$$`=`$ 4506 K ($`Z=`$ 0.01), $`T_{\mathrm{eff},\mathrm{A}}`$$`=`$ 5682 K and $`T_{\mathrm{eff},\mathrm{B}}`$$`=`$ 4586 K ($`Z=`$ 0.02), $`T_{\mathrm{eff},\mathrm{A}}`$$`=`$ 5754 K and $`T_{\mathrm{eff},\mathrm{B}}`$$`=`$ 4638 K ($`Z=`$ 0.03). The uncertainty on the primary (secondary) component is less than 150 K (100 K), and so – to be conservative – we consider these upper limits as the uncertainty. For further comparisons, we take the $`T_{\mathrm{eff}}`$s derived for $`Z=`$0.02. They are cooler by about 100 K (170 K) for the primary (secondary) component in comparison to the $`T_{\mathrm{eff}}`$s that we have used previously (see Table 1).

The influence of adopting the $`T_{\mathrm{eff}}`$s we have derived from the Alonso et al. calibrations can be seen in Fig. 4. A comparison with Fig. 2 shows that the global solutions (dots) are fully consistent with each other. However, the best solution according to the criteria defined in Section 3 is slightly different: ($`Z`$, $`Y`$) $`=`$ (0.015, 0.295), the ages predicted being about twice older (UV Psc A: $``$ 4.07 Gyr and B: $``$ 4.14 Gyr). In spite of these shifts, we note that a different value is still needed for the convection parameter of the two components ($`\alpha _{\mathrm{MLT},\mathrm{A}}`$$`=`$ 1.15 and $`\alpha _{\mathrm{MLT},\mathrm{B}}`$$`=`$ 0.65). The difference $`\mathrm{\Delta }`$$`\alpha _{\mathrm{MLT}}`$ is even increased in comparison to the results inferred from the Popper (1997)’s $`T_{\mathrm{eff}}`$ scale, and hence the conclusions of the previous sub-section remain qualitatively unchanged. Quantitatively, a systematic error of about +0.20 in $`\alpha _{\mathrm{MLT},\mathrm{A}}`$ (and a negligible one in $`\alpha _{\mathrm{MLT},\mathrm{B}}`$) may be due to the uncertainties in the temperature scaling.

### 3.2 Influence of diffusion

Another potential source of systematic error is the influence of diffusion. We have treated this by computing CESAM models incorporating the microscopic diffusion as described by the simplified formalism proposed by Michaud & Proffitt (1993) with metal elements as trace elements. Note that the radiative acceleration can be neglected for sub-solar stars (Turcotte et al., 1998). As shown on Fig. 5, diffusion almost mimics an increase in $`\alpha _{\mathrm{MLT}}`$, even though the slopes of the tracks are obviously changed. At fixed radius, the effect of diffusion is to change the effective temperatures by some 65 K (55 K), well within the 1-$`\sigma `$ uncertainties in $`T_{\mathrm{eff}}`$ discussed above. Alternatively, taking these uncertainties at face value, one would have to increase $`\alpha _{\mathrm{MLT}}`$by about 0.11 (0.07) at most.

In summary, even though the absolute values of $`\alpha _{\mathrm{MLT}}`$ may change due to these combined effects, the trend of increased $`\alpha _{\mathrm{MLT}}`$with mass remains robust.

### 3.3 Universality of the $`\alpha _{\mathrm{MLT}}`$ parameter ?

We have quantitatively confirmed the suggestion made by Clausen et al. (1999) that the adoption of a significantly smaller $`\alpha _{\mathrm{MLT}}`$ for the secondary component may remove (or at least decrease) the general problem of age discrepancy presented in the Introduction. While it is obvious that only observational constraints on the chemical composition of UV Psc would allow a proper test, it seems that this problem can be solved theoretically for this binary. The solution obtained for this system suggests that i) the components have a different $`\alpha _{\mathrm{MLT}}`$ ($`\mathrm{\Delta }`$$`\alpha _{\mathrm{MLT}}`$$`=`$0.3$`\pm `$0.14) and ii) these values are clearly different from the solar one.

This implies that the solar value of the $`\alpha _{\mathrm{MLT}}`$ is not universal, an hypothesis made in most of the currently available stellar evolutionary models.

As for the trend with mass, a recent study on the Hyades eclipsing binary vB22 (Lebreton et al. 2001) seems to support it. We show in Fig. 6 that UV Psc and vB22 are clearly not consistent with the solar $`\alpha _{\mathrm{MLT}}`$ value. Moreover, Lebreton et al. (2001) show that the slope of the main-sequence in the Hyades suggests that $`\alpha _{\mathrm{MLT}}`$ could decrease with mass, in particular they obtained $`\alpha _{\mathrm{MLT}}`$$``$1.4 below 1 $`M_{}`$. The results we have obtained for UV Psc give additional support to this trend. However we note that this trend seems to contradict the results of Ludwig et al. (1999) coming from detailed hydrodynamical simulations where, for a fixed gravity, $`\alpha _{\mathrm{MLT}}`$ increases with deacreasing $`T_{\mathrm{eff}}`$.

More results are needed before reaching any definitive conclusion, because a larger sample may well reveal a scatter of $`\alpha _{\mathrm{MLT}}`$ for a given mass. For instance, recent calibrations of $`\alpha _{\mathrm{MLT}}`$ in visual binary stars were performed by P. Morel and collaborators using the CESAM code (with different physical ingredients than those used in this paper). For the stars in the range of mass relevant in this paper, they find $`\alpha `$ Cen B, $`M=`$ 0.97 $`M_{}`$, $`\alpha _{\mathrm{MLT}}`$$`=`$ 1.53 (Morel et al 2000a), $`\iota `$ Peg B, $`M=`$0.81 $`M_{}`$, $`\alpha _{\mathrm{MLT}}`$$`=`$1.36 (Morel et al 2000b) and 85 Peg A, $`M=`$0.84 $`M_{}`$, $`\alpha _{\mathrm{MLT}}`$$`=`$1.80 (Fernandes et al. 2002). Since both $`\iota `$ Peg B and 85 Peg A have very similar masses, this seems to support the idea of a dispersion of $`\alpha _{\mathrm{MLT}}`$ at a given mass. Solar models with and without diffusion yield different values in the internal convective regions (1.7 vs 1.9, Richard et al. 1996). There are no further data for main sequence stars with larger masses. Studies of more massive stars have dealt exclusively with red giant stars. For instance, $`\alpha _{\mathrm{MLT}}`$ seems to increase to values up to 3 to reproduce the red giant branches in open clusters (Stothers & Chin, 1997), while Iwamoto & Saio (1999) present a trend of $`\alpha _{\mathrm{MLT}}`$ increasing with mass and metallicity again for giant stars. In these stars overshooting is playing a major role and so disentangling the relative effects of mass, metallicity and overshooting is far from trivial.

Therefore, at this stage, a relation between mass and $`\alpha _{\mathrm{MLT}}`$ along the main sequence remains highly premature and must be seen with caution.

We also note, in passing, that our results are not in contradiction with the recent analysis by Palmieri et al. (2002), who on the basis of the colours of red giants in globular clusters concluded that there is no dependence of $`\alpha _{\mathrm{MLT}}`$ with metallicity. The extremely narrow range in mass probed by these red giants would be nearly equivalent to a single point on our diagram. Since the detailed modelling of red giants is more complicated than the one of main sequence stars, a larger sample of double-lined, eclipsing binaries with components in the 0.7-1.1 $`M_{}`$would offer a much cleaner and definitive test.

In any case, all these possible variations and their trends reflect a shortcoming of the MLT theory more than a real physical effect. If one insists in using MLT, its parameter $`\alpha _{\mathrm{MLT}}`$ should be varied. This is perhaps the first time that such a result is obtained using observations, since the variablity of the effective $`\alpha _{\mathrm{MLT}}`$ is a well-know result of detailed hydrodynamical simulations of stellar convection (Ludwig et al. 1999).

On the other hand we note that the value of $`\alpha _{\mathrm{MLT}}`$ determined in atmosphere of the stars gives in general different results in relation to those obtained from the stellar interior analysis. For instance, van’t Veer-Menneret & Mégessier (1996) found that the spectral fitting of the first for Balmer lines of the solar spectrum required $`\alpha _{\mathrm{MLT}}`$$`=0.5`$, a result confirmed by Fuhrmann et al. (1993) for other cool stars, while Gardiner, Kupka & Smalley (1999) found that $`\alpha _{\mathrm{MLT}}`$ $`0.5`$ for $`T_{\mathrm{eff}}`$s below 6000 K but $`\alpha _{\mathrm{MLT}}`$ $`1.25`$ in the range 6000–7000 K to possibly decrease again to $`\alpha _{\mathrm{MLT}}`$ $`=0.5`$ above 7000 K.

### 3.4 Implications for the helium to metal ratio

Since we have derived $`Z`$ and $`Y`$ independently for UV Psc, the helium to metal ratio ($`\mathrm{\Delta }`$$`Y`$/$`\mathrm{\Delta }`$$`Z`$) of this system can be inferred assuming a primordial helium abundance. The initial helium content is fixed by taking the currently assumed law $`Y=`$ $`Y_p`$ $`+`$ ($`\mathrm{\Delta }`$$`Y`$/$`\mathrm{\Delta }`$$`Z`$) $`\times `$ $`Z`$, where $`Y_p`$ is the primordial helium abundance. If we consider $`Y_p=`$0.235 (Peimbert et al. 2000) and use the ($`Z`$, $`Y`$) solution that we have found for the UV Psc binary, we derive $`\mathrm{\Delta }`$$`Y`$/$`\mathrm{\Delta }`$$`Z`$ $``$ 6 and 4 by using the Popper (1997) or Alonso et al. (1996) $`T_{\mathrm{eff}}`$s, respectively. In both cases, the ratio is large in comparison to what can be obtained from the solar model described previously: ($`\mathrm{\Delta }`$$`Y`$/$`\mathrm{\Delta }`$$`Z`$) $``$ 1.9. This may suggest that the currently assumed simple linear relationship between $`Z`$ and $`Y`$ is not universal and/or that $`\mathrm{\Delta }`$$`Y`$/$`\mathrm{\Delta }`$$`Z`$ is not unique.
The possible relation between $`Y`$ and $`Z`$ is one of the most interesting and puzzling problems in astrophysics since the pioneering work by Faulkner (1967).

The most recent observations in HII regions and planetary nebulae as well as semi-empical determinations using theoretical stellar models - including the detailed calibration of the Sun - (cf. Ribas et al. 2000, Fernandes 2001 and references therein) support values around 2 or 3 , assuming $`Y_p`$ $``$ 0.23. Nevertheless, some exceptions indicate values larger than 4 (e.g. Pagel et al. 1992 from observations of extragalactic HII regions, and Belikov et al. 1998 from the main sequence of the Pleiades). Our determination clearly supports the idea of a large $`\mathrm{\Delta }`$$`Y`$/$`\mathrm{\Delta }`$$`Z`$ ratio.

## 4 Conclusions

We study the eclipsing binary UV Psc, which is of particular interest since Popper (1997) derived age estimates for both components which differed by some 11 Gyr. Adopting his same data set (masses, radii and $`T_{\mathrm{eff}}`$s) we compute a large grid of CESAM models for each component555These models are available upon request (email: jmfernan@mat.uc.pt).. We obtain a best fit in the $`T_{\mathrm{eff}}`$-radius diagram for a common chemical composition ($`Z`$, $`Y`$)$`=`$(0.012, 0.31), a negligible age difference, and a different MLT parameter ($`\alpha _{\mathrm{MLT},\mathrm{A}}`$$`=`$ 0.95$`\pm `$0.12 and $`\alpha _{\mathrm{MLT},\mathrm{B}}`$$`=`$ 0.65$`\pm `$0.07). The difference $`\alpha _{\mathrm{MLT},\mathrm{A}}`$ $``$ $`\alpha _{\mathrm{MLT},\mathrm{B}}`$ $`=`$ 0.3 is significant and remains robust against changes in temperature calibrations and the inclusion of diffusion.

Within the framework of the MLT theory, this suggests that fixing $`\alpha _{\mathrm{MLT}}`$ to its solar value ($``$1.6), the usual hypothesis made by most stellar evolutionary models, may not be correct. Secondly, since the $`\alpha _{\mathrm{MLT}}`$ is smaller for the less massive component, this suggests that the $`\alpha _{\mathrm{MLT}}`$ parameter may decrease with the stellar mass. This trend needs further confirmation with other binary stars with accurate data. Again, this may reflect a problem of the application of the MLT theory, more than a real physical effect, but it has to be taken into account when using MLT. On the other hand as the main observational errors come from the $`T_{\mathrm{eff}}`$, new and independent determinations are needed for both stars.

So far, few other binaries have data with the precision required to confirm this trend and observational programs are highly needed (e.g. Kurpinska-Winiarska & Oblak, 2000).

###### Acknowledgements.

We warmly thank P. Morel for making available CESAM, a superb tool for stellar evolution, and without which this study could not have been performed. EL is supported by a ”Fundação para a Ciência e Tecnologia” (FCT) postdoctoral fellowship (grant SFRH/BPD/5556/2001). This work was partially supported by the project ”PESO/P/PRO/15128/1999” from the FCT, and by ”Convénio Embaixada de França-ICCTI (59-B0)”. This research has made use of the SIMBAD database operated at CDS, Strasbourg, France, and of NASA’s Astrophysics Data System Abstract Service.


format: text

text:
# 1 Introduction

## 1 Introduction

Magnetic monopole solutions are a generic prediction of grand unified theories. Such solutions possess a topological charge proportional to their magnetic charge. In Yang-Mills-Higgs (YMH) theory with gauge group SU(2) the simplest solution with unit topological charge is the spherically symmetric ’t Hooft-Polyakov monopole . SU(2) monopoles with higher topological charge cannot be spherically symmetric  and possess at most axial symmetry .

In the Bogomol’nyi-Prasad-Sommerfield (BPS) limit of vanishing Higgs potential spherically symmetric monopole and axially symmetric multimonopole solutions, which satisfy the first order Bogomol’nyi equations  as well as the second order field equations, are known analytically . For these solutions all zeros of the Higgs field are superimposed at a single point. Multimonopole solutions of the Bogomol’nyi equations which do not possess any rotational symmetry , have recently been constructed numerically . In these solutions the zeros of the Higgs field are no longer all superimposed at a single point but are located at several isolated points.

As shown by Taubes, in each topological sector there exist in addition smooth, finite energy solutions of the second order field equations, which do not satisfy the Bogomol’nyi equations . Consequently, the energy of these solutions exceeds the Bogomol’nyi bound. The simplest such solution resides in the topologically trivial sector and forms a saddlepoint of the energy functional . It possesses axial symmetry, and the two zeros of its Higgs field are located symmetrically on the positive and negative $`z`$-axis. This solution corresponds to a monopole and antimonopole in static equilibrium .

In this paper we present new axially symmetric saddlepoint solutions, where the Higgs field vanishes at $`m>2`$ isolated points on the symmetry axis. These configurations represent chains of monopoles and antimonopoles, in alternating order. Chains with an equal number of monopoles and antimonopoles, $`m=2k`$, reside in the topologically trivial sector, whereas chains with $`k+1`$ monopoles and $`k`$ antimonopoles reside in the sector with topological charge one.

After briefly reviewing SU(2) Yang-Mills-Higgs theory in Section 2, we discuss the axially symmetric Ansatz in Section 3 together with the boundary conditions. We present the numerical results in Section 4 and our conclusions in Section 5.

## 2 Yang-Mills-Higgs theory

The Lagrangian of SU(2) Yang-Mills-Higgs theory is given by

$$L=\left\{\frac{1}{2}\mathrm{Tr}\left(F_{\mu \nu }F^{\mu \nu }\right)+\frac{1}{4}\mathrm{Tr}\left(D_\mu \mathrm{\Phi }D^\mu \mathrm{\Phi }\right)+\frac{\lambda }{8}\mathrm{Tr}\left[\left(\mathrm{\Phi }^2\eta ^2\right)^2\right]\right\}d^3r,$$
(1)
with su(2) gauge potential $`A_\mu =A_\mu ^a\tau ^a/2`$, field strength tensor $`F_{\mu \nu }=_\mu A_\nu _\nu A_\mu +ie[A_\mu ,A_\nu ]`$, and covariant derivative of the Higgs field $`D_\mu \mathrm{\Phi }=_\mu \mathrm{\Phi }+ie[A_\mu ,\mathrm{\Phi }]`$. $`e`$ denotes the gauge coupling constant, $`\eta `$ the vacuum expectation value of the Higgs field and $`\lambda `$ the strength of the Higgs selfcoupling.

The field equations are derived from the Lagrangian as the variational equations with respect to the gauge potential and the Higgs field,

$`D_\mu F^{\mu \nu }{\displaystyle \frac{1}{4}}ie[\mathrm{\Phi },D^\nu \mathrm{\Phi }]`$ $`=`$ $`0,`$ (2)
$`D_\mu D^\mu \mathrm{\Phi }\lambda \left(\mathrm{\Phi }^2\eta ^2\right)\mathrm{\Phi }`$ $`=`$ $`0.`$ (3)
## 3 Ansatz and Boundary Conditions

We parametrize the gauge potential and the Higgs field by the Ansatz

$`A_\mu dx^\mu `$ $`=`$ $`\left({\displaystyle \frac{K_1}{r}}dr+(1K_2)d\theta \right){\displaystyle \frac{\tau _\phi }{2e}}\mathrm{sin}\theta \left(K_3{\displaystyle \frac{\tau _r^{(m)}}{2e}}+(1K_4){\displaystyle \frac{\tau _\theta ^{(m)}}{2e}}\right)d\phi `$ (4)
$`\mathrm{\Phi }`$ $`=`$ $`\mathrm{\Phi }_1\tau _r^{(m)}+\mathrm{\Phi }_2\tau _\theta ^{(m)}.`$ (5)
which generalizes the axially symmetric Ansatz employed in  for the monopole-antimonopole pair solution. Here the su(2) matrices $`\tau _r^{(m)}`$, $`\tau _\theta ^{(m)}`$, and $`\tau _\phi `$ are defined as

$`\tau _r^{(m)}`$ $`=`$ $`\mathrm{sin}(m\theta )\tau _\rho +\mathrm{cos}(m\theta )\tau _z,`$
$`\tau _\theta ^{(m)}`$ $`=`$ $`\mathrm{cos}(m\theta )\tau _\rho \mathrm{sin}(m\theta )\tau _z,`$
$`\tau _\phi `$ $`=`$ $`\mathrm{sin}\phi \tau _x+\mathrm{cos}\phi \tau _y,`$
where $`\tau _\rho =\mathrm{cos}\phi \tau _x+\mathrm{sin}\phi \tau _y`$, and $`m`$ is an integer to which we refer as $`\theta `$ winding number.

The profile functions $`K_1`$$`K_4`$ and $`\mathrm{\Phi }_1`$, $`\mathrm{\Phi }_2`$ depend on the coordinates $`r`$ and $`\theta `$, only. With this Ansatz the general field equations Eqs. (2) and (3) reduce to six PDEs in the coordinates $`r`$ and $`\theta `$. The ansatz possesses a residual U(1) gauge symmetry. To fix the gauge we impose the condition $`r_rK_1_\theta K_2=0`$. For convenience we change to dimensionless coordinates by rescaling $`rr/(e\eta )`$ and $`\mathrm{\Phi }\eta \mathrm{\Phi }`$.

To obtain regular solutions with finite energy and energy density we have to impose appropriate boundary conditions. Regularity at the origin requires

$$K_1(0,\theta )=0,K_2(0,\theta )=1,K_3(0,\theta )=0,K_4(0,\theta )=1,$$

$$\mathrm{sin}(m\theta )\mathrm{\Phi }_1(0,\theta )+\mathrm{cos}(m\theta )\mathrm{\Phi }_2(0,\theta )=0,$$

$$_r\left[\mathrm{cos}(m\theta )\mathrm{\Phi }_1(r,\theta )\mathrm{sin}(m\theta )\mathrm{\Phi }_2(r,\theta )\right]|_{r=0}=0.$$

To obtain the boundary conditions at infinity we require that solutions in the vacuum sector ($`m=2k`$) tend to a gauge transformed trivial solution,

$$\mathrm{\Phi }U\tau _zU^{},A_\mu i_\mu UU^{},$$

and that solutions in the sector with unit topological charge ($`m=2k+1`$) tend to

$$\mathrm{\Phi }U\mathrm{\Phi }_{\mathrm{}}^{(1)}U^{},A_\mu UA_\mu \mathrm{}^{(1)}U^{}+i_\mu UU^{},$$

where

$$\mathrm{\Phi }_{\mathrm{}}^{(1)}=\tau _r^{(1)},A_\mu \mathrm{}^{(1)}dx^\mu =\frac{\tau _\phi }{2}d\theta \mathrm{sin}\theta \frac{\tau _\theta ^{(1)}}{2}d\phi $$

is the asymptotic solution of a monopole, and $`U=\mathrm{exp}\{ik\theta \tau _\phi \}`$, both for even and odd $`m`$. Consequently, solutions with even $`m`$ have vanishing magnetic charge, whereas solutions with odd $`m`$ possess unit magnetic charge.

In terms of the functions $`K_1K_4`$, $`\mathrm{\Phi }_1`$ and $`\mathrm{\Phi }_2`$ these boundary conditions read

$$K_10,K_21m,$$
(6)
$$K_3\frac{\mathrm{cos}\theta \mathrm{cos}(m\theta )}{\mathrm{sin}\theta }m\mathrm{odd},K_3\frac{1\mathrm{cos}(m\theta )}{\mathrm{sin}\theta }m\mathrm{even},$$
(7)
$$K_41\frac{\mathrm{sin}(m\theta )}{\mathrm{sin}\theta },$$
(8)
$$\mathrm{\Phi }_11,\mathrm{\Phi }_20.$$
(9)
Regularity on the $`z`$-axis requires

$$K_1=K_3=\mathrm{\Phi }_2=0,_\theta K_2=_\theta K_4=_\theta \mathrm{\Phi }_1=0,$$

for $`\theta =0`$ and $`\theta =\pi `$.

## 4 Results

We have constructed numerically axially symmetric solutions of the Yang-Mills-Higgs theory in the BPS limit $`\lambda =0`$ for $`\theta `$ winding number $`1m6`$.

These $`m`$-chains possess $`m`$ nodes of the Higgs field on the $`z`$-axis. Due to reflection symmetry, each node on the negative $`z`$-axis corresponds to a node on the positive $`z`$-axis. The nodes of the Higgs field are associated with the location of the monopoles and antimonopoles. For odd $`m`$ ($`m=2k+1`$) the Higgs field possesses $`k`$ nodes on the positive $`z`$-axis and one node at the origin. The node at the origin corresponds to a monopole if $`k`$ is even and to an antimonopole if $`k`$ is odd. For even $`m`$ ($`m=2k`$) there is no node of the Higgs field at the origin.

The $`m=1`$ solution is the well-known ’t Hooft-Polyakov monopole . The $`m=3`$ (M-A-M) and $`m=5`$ (M-A-M-A-M) chains represent saddlepoints with unit topological charge. The $`m=2`$ chain is identical to the monopole-antimonopole (M-A) pair discussed in . The M-A pair as well as the $`m=4`$ (M-A-M-A) and $`m=6`$ (M-A-M-A-M-A) chains form saddlepoints in the vacuum sector.

In Table 1 we present the dimensionless energy $`E/4\pi \eta `$, the dimensionless magnetic dipole moment $`\mu `$ and the nodes $`z_i`$ of the Higgs field for the solutions with $`\theta `$ winding number $`1m6`$.

We observe that the energy $`E^{(m)}`$ of an $`m`$-chain is always smaller than the energy of $`m`$ single monopoles or antimonopoles (with infinite separation between them), i. e. $`E^{(m)}<E_{\mathrm{}}/4\pi \eta =m`$. On the other hand $`E^{(m)}`$ exceeds the minimal energy bound given by the Bogolmol’ny limit $`E_{\mathrm{min}}/4\pi \eta =0`$ for even $`m`$, and $`E_{\mathrm{min}}/4\pi \eta =1`$ for odd $`m`$. We observe an (almost) linear dependence of the energy $`E^{(m)}`$ on $`m`$. This can be modelled by taking into account only the energy of $`m`$ single (infinitely separated) monopoles and the next-neighbour interaction between monopoles and antimonopoles on the chain. Defining the interaction energy as the binding energy of the monopole-antimonopole pair,

$$\mathrm{\Delta }E=2(4\pi \eta )E^{(2)},$$

we obtain as energy estimate for the $`m`$-chain

$$E_{\mathrm{est}}^{(m)}/4\pi \eta =m+(m1)\mathrm{\Delta }E.$$

In Fig. 1 we show this estimate for the energy, $`E_{\mathrm{est}}^{(m)}/4\pi \eta `$, together with the exact energy for chains with $`\theta `$ winding number $`m=1,\mathrm{},6`$. We note that the deviation of the estimated energy from the exact energy is indeed very small.

In Table 2 we list the possible decompositions of the $`m`$-chains into subchains of smaller length together with the total energy of the subchains (at infinite separation). We observe that the energetically most favourable state is the $`m`$-chain. This supports our interpretation of the $`m`$-chain as an equilibrium state of $`m`$ monopoles and antimonopoles. We note, however, that there are other possible decompositions not included in the table. For example, the M-A-M-A-chain could be decomposed into a pair of a charge $`2`$ multimonopole and a charge $`2`$ anti-multimonopole.

To define the magnetic dipole moment for solutions with even $`m`$, we first transform to a gauge where the Higgs field is constant at infinity, $`\mathrm{\Phi }=\tau _z`$. From the asymptotic expansion , we obtain $`K_3(1\mathrm{cos}(m\theta ))/\mathrm{sin}(\theta )+C_3\mathrm{sin}\theta /r`$. Thus the gauge field assumes the asymptotic form

$$A_\mu dx^\mu =C_3\frac{\mathrm{sin}^2\theta }{2r}\tau _zd\phi ,$$
(10)
from which we read off the (dimensionless) magnetic dipole moment $`\mu =C_3`$. Solutions with odd $`m`$ have vanishing magnetic dipole moment, since in this case the function $`K_3`$ is odd under the transformation $`zz`$. Consequently, the asymptotic form of the gauge potential cannot contain terms like the r.h.s. of Eq. (10).

The magnetic dipole moment also increases with increasing $`m`$, as seen in Table 1. To obtain an estimate for the dipole moment we consider the magnetic charges as point charges located at the nodes of the Higgs field, yielding

$$\mu _{\mathrm{est}}^{(m)}=\underset{i=1,m}{}z_iP_i,$$

with charges $`P_i=1`$ for monopoles and $`P_i=1`$ for antimonopoles, respectively. For comparison the estimated magnetic dipole moments are also given in Table 1. The deviation from the exact values is on the order of $`10`$%.

Concerning the nodes of the Higgs field we observe that the distances between the nodes increase with increasing $`m`$. Remarkably, the distances between the nodes do not vary much within a chain. For example, denoting the location of the nodes by $`z_i`$ in decreasing order we find for the chain with $`\theta `$ winding number $`m=6`$ from Table 1, $`|z_1z_2|5.06`$, $`|z_2z_3|5.11`$ and $`|z_3z_4|4.92`$.

In Fig. 2 we present the dimensionless energy density for the solutions with $`\theta `$ winding number $`m=1,\mathrm{},6`$. The energy density of the $`m`$-chain possesses $`m`$ maxima on the $`z`$-axis, and decreases with increasing $`\rho `$. The locations of the maxima are close to the nodes of the Higgs field, which are indicated by asterisks. For a given $`m`$ the maxima are of similar magnitude, but their height decreases with increasing $`m`$. (Note that the scale for the $`m=1`$ solution is different compared to the $`m2`$ solutions, and that the contour lines are distorted due to the different scaling of the $`\rho `$\- and $`z`$-axis.)

## 5 Conclusion

We have obtained new static axially symmetric solutions of the SU(2) Yang-Mills-Higgs theory which represent monopole-antimonopole chains. They are characterized by the $`\theta `$ winding number $`m`$, which equals the number of nodes of the Higgs field, and the total number of monopoles and antimonopoles. Solutions with even $`m`$ carry no magnetic charge but possess a non-vanishing magnetic dipole moment, whereas solutions with odd $`m`$ carry unit magnetic charge but possess no magnetic dipole moment. The energy of these $`m`$-chains increases (almost) linearly with $`m`$.

We interpret the $`m`$-chains as equilibrium states of $`m`$ monopoles and antimonopoles. As shown long ago , the force between monopoles is given by twice the Coloumb force when the charges are unequal, and vanishes when the charges are equal, provided the monopoles are at large distances. Thus, monopoles and antimonopoles can only be in static equilibrium, if they are close enough to experience a repulsive force that counteracts the attractive Coloumb force. In other words, $`m`$-chains are essentially non-BPS solutions. To see this in another way let us rewrite the energy in the form

$$E=\left\{\frac{1}{4}\mathrm{Tr}\left(\left(\epsilon _{ijk}F_{ij}\pm D_k\mathrm{\Phi }\right)^2\right)\frac{1}{2}\epsilon _{ijk}\mathrm{Tr}\left(F_{ij}D_k\mathrm{\Phi }\right)\right\}d^3r.$$
(11)
The second term is proportional to the topological charge and vanishes when $`m`$ is even. The first term is just the integral of the square of the Bogol’molnyi equations. Thus, for even $`m`$ the energy is a measure for the deviation of the solution from selfduality.

So far we have considered only solutions in the vacuum sector or the sector with unit topological charge. Generalising the ansatz to $`\phi `$ winding number $`n>1`$ leads to new solutions, which carry topological charge $`n`$ if $`m`$ is odd, and to solutions in the vacuum sector if $`m`$ is even . Recently, the $`n=2`$, $`m=2`$ solution has been constructed in an extended model , which includes the YMH theory as a special case. Chains of $`n=2`$ multimonopoles and anti-multimonopoles will be presented elsewhere . For $`n>2`$ a new phenomenon occurs. The zeros of the Higgs field no longer form a set of isolated points. Instead the Higgs field vanishes on rings around the $`z`$-axis . (For odd $`m`$ the node at the origin persists).

Dyonic solutions can be readily obtained from the $`m`$-chains , as outlined in . Interestingly, these solutions carry electric charge even for solutions in the vacuum sector.

When the gravitional interaction is included, we anticipate a different behaviour for $`m`$-chains with finite magnetic charge and those with vanishing magnetic charge. For magnetically charged solutions a degenerate horizon may form at a critical value of the gravitional parameter, as observed for monopoles  and multimonopoles . On the other hand, no formation of a horizon was found for the gravitating monopole-antimonopole pair .

We expect that solutions analogous to the $`m`$-chains exists in Weinberg-Salam theory , generalizing the sphaleron-antisphaleron pair . The axially symmetric Ansatz with $`\phi `$ winding number $`n`$ and $`\theta `$ winding number $`m`$  allows for (multi)sphaleron–anti-(multi)sphaleron chains and solutions with rings of zeros.

| m | $`E[4\pi \eta ]`$ | $`\mu `$ | $`\mu _{\mathrm{est}}`$ | $`z_i`$ |
| --- | --- | --- | --- | --- |
| 1 | $`1.00`$ | $`0.0`$ | $`0.0`$ | $`0.0`$ |
| 2 | $`1.70`$ | $`4.72`$ | $`4.18`$ | $`\pm 2.09`$ |
| 3 | $`2.44`$ | $`0.0`$ | $`0.0`$ | $`0.0`$ $`\pm 4.67`$ |
| 4 | $`3.12`$ | $`9.86`$ | $`9.21`$ | $`\pm 2.39`$ $`\pm 6.99`$ |
| 5 | $`3.78`$ | $`0.0`$ | $`0.0`$ | $`0.0`$ $`\pm 4.79`$ $`\pm 9.61`$ |
| 6 | $`4.40`$ | $`16.06`$ | $`15.40`$ | $`\pm 2.46`$ $`\pm 7.57`$ $`\pm 12.63`$ |

Table 1 The dimensionless energy, the dipole moment $`\mu `$, the estimated dipole moment $`\mu _{\mathrm{est}}^{(m)}`$, and the nodes $`z_i`$ are given for the $`m`$-chains with $`m=1,\mathrm{},6`$.

| m | chain | $`E[4\pi \eta ]`$ | decomposition | $`E[4\pi \eta ]`$ |
| --- | --- | --- | --- | --- |
| 2 | M-A | $`E^{(2)}=1.70`$ | (M)+(A) | $`2E^{(1)}=2.00`$ |
| 3 | M-A-M | $`E^{(3)}=2.44`$ | (M-A)+(M)  (M)+(A)+(M) | $`E^{(2)}+E^{(1)}=2.7`$  $`3E^{(1)}=3.00`$ |
| 4 | M-A-M-A | $`E^{(4)}=3.12`$ | (M-A)+(M-A)  (M-A-M)+(A)  (M-A)+(M)+(A)  (M)+(A)+(M)+(A) | $`2E^{(2)}=3.40`$  $`E^{(3)}+E^{(1)}=3.44`$  $`E^{(2)}+2E^{(1)}=3.70`$  $`4E^{(1)}=4.00`$ |
| 5 | M-A-M-A-M | $`E^{(5)}=3.78`$ | (M-A-M-A)+(M)  (M-A-M)+(A-M)  (M-A)+(M-A)+(M)  (M-A-M)+(A)+(M)  (M-A)+(M)+(A)+(M)  (M)+(A)+(M)+(A)+(M) | $`E^{(4)}+E^{(1)}=4.12`$  $`E^{(3)}+E^{(2)}=4.14`$  $`2E^{(2)}+E^{(1)}=4.40`$  $`E^{(3)}+2E^{(1)}=4.44`$  $`E^{(2)}+3E^{(1)}=4.7`$  $`5E^{(1)}=5.00`$ |
| 6 | M-A-M-A-M-A | $`E^{(6)}=4.40`$ | (M-A-M-A-M)+(A)  (M-A-M-A)+(M-A)  (M-A-M)+(A-M-A)  (M-A)+(M-A)+(M-A)  (M-A-M-A)+(M)+(A)  (M-A-M)+(A-M)+(A)  (M-A)+(M-A)+(M)+(A)  (M-A-M)+(A)+(M)+(A)  (M-A)+(M)+(A)+(M)+(A)  (M)+(A)+(M)+(A)+(M)+(A) | $`E^{(5)}+E^{(1)}=4.78`$  $`E^{(4)}+E^{(2)}=4.82`$  $`2E^{(3)}=4.88`$  $`3E^{(2)}=5.10`$  $`E^{(4)}+2E^{(1)}=5.12`$  $`E^{(3)}+E^{(2)}+E^{(1)}=5.14`$  $`2E^{(2)}+2E^{(1)}=5.40`$  $`E^{(3)}+3E^{(1)}=5.44`$  $`E^{(2)}+4E^{(1)}=5.70`$  $`6E^{(1)}=6.00`$ |

Table 2 The decompositions of the $`m`$-chains into subchains and their energies are given for $`m=1,\mathrm{},6`$.

Fig. 1 The exact energy (asterisks) and the estimate for the energy (crosses) are shown for chains with $`m=1,\mathrm{},6`$. The solid line demonstrates the linear dependence of $`E_{\mathrm{est}}^{(m)}`$ on $`m`$.

Fig. 2 The dimensionless energy density is shown as function of $`\rho `$ and $`z`$ for chains with $`m=1,\mathrm{},6`$. The asterisks indicate the nodes of the Higgs field.


text:
# Motion of a heavy impurity through a Bose-Einstein condensate.

## I Introduction

One of the most important peculiarities of Landau theory superfluidity is the existence of a finite critical velocity. If a body moves in a superfluid at $`T=0`$ with velocity $`V`$ less then $`v_c`$, the motion is dissipationless. At $`V>v_c`$ a drag force arises because of the possibility of emission of elementary excitations. However, both theoretical and experimental investigation in superfluid $`{}_{}{}^{4}He`$ are difficult. The critical velocity in $`{}_{}{}^{4}He`$ is related to creation of rotons, for which one has no simple theoretical description. Further, an important role is played by complicated processes involving vortex rings production.

The situation in low-density weakly-interacting Bose-Einstein condensed (BEC) gases is simpler. The Landau critical velocity in this case is due to Cherenkov emission of phonons which can be described by mean-field theory. Due to the presence in the theory of an intrinsic length parameter - the correlation length $`\xi `$\- the friction force for a small body does not depend on its structure. Vortex rings in the BEC cannot have radius less then $`\xi `$ (see ) and it is reasonable to believe that probability of their creation by a small body is small. Thus quantitative investigation of critical velocities in BEC are very interesting and can be used to probe the superfluidity of a quantum gas.

Recently existence of the critical velocity in a Bose-Einstein Condensed gas was confirmed in a few experiments. At MIT a trapped condensate was stirred by a blue detuned laser beam  and the energy of dissipation was measured. The critical velocity was found to be smaller than the speed of sound due to emission of vortices. The diameter of the laser spot in this experiment was of a macroscopic size and was large compared to the healing length. An improved technique allowed measurement of the drag force acting on the condensate in a subsequent experiment .

The analytical study of flow of the condensate over an impurity is highly nontrivial due to the intrinsic nonlinearity of the problem arising from the interaction of the particles in the condensate. In one dimension the dissipation could occur at velocities smaller than predicted by Landau’s approach due to emission of solitons . The dependence of the critical velocity on the type of the potential was studied both by using a perturbative approach and numerical integration in . The effective two dimensional problem was considered in . In this work generation of excitations in the oscillating condensate in a time dependent parabolic trap in the presence of a static impurity was studied analytically. A three-dimensional flow of a condensate around an obstacle was calculated numerically by integration of the Gross-Pitaevskii (GP) equation and emission of vortices was observed .

In this letter we study a $`\delta `$-function perturbation moving at a constant velocity in a condensate. We find analytically the depletion of the superfluid fraction and the drag force.

## II Three-dimensional system

Let us consider an impurity moving through a three-dimensional condensate at $`T=0`$. One of the possible realizations of this model could be scattering of heavy neutral molecules by the condensate.

### A Effective mass and normal fraction

We start from the three-dimensional energy functional of a homogeneous weakly-interacting Bose gas in the presence of a $`\delta `$-function perturbation (an impurity) moving with a constant velocity $`𝐕`$

$$E=(\frac{\mathrm{}^2}{2m}|\psi |^2+(\mu g_i\delta (𝐫𝐕t))|\psi |^2+\frac{g}{2}|\psi |^4)d^3x,$$
(1)
where $`\psi `$ is the condensate wavefunction, $`\mu `$ is the chemical potential, $`m`$ mass of a particle in the condensate, $`g=4\pi \mathrm{}^2a/m`$ and $`g_i=2\pi \mathrm{}^2b/m`$ are particle-particle and particle-impurity coupling constants, with $`a`$ and $`b`$ being the respective scattering lengths . We will assume that the interaction with impurity is small and we will use perturbation theory. By splitting the wave function into a sum of the unperturbed solution and a small correction $`\psi (𝐫,t)=\varphi _0+\delta \psi (𝐫,t)`$ and linearizing the time-dependent GP equation with respect to $`\delta \psi `$, we obtain an equation describing the time evolution of $`\delta \psi `$

$`i\mathrm{}{\displaystyle \frac{}{t}}\delta \psi `$ $`=`$ $`\left({\displaystyle \frac{\mathrm{}^2}{2m}}\mathrm{}\mu +2g|\varphi _0|^2\right)\delta \psi +g\varphi _0^2\delta \psi ^{}`$ (3)
$`+g_i\delta (𝐫𝐕t)\varphi _0.`$
In a homogeneous system $`\varphi _0`$ is a constant fixed by the particle density $`\varphi _0=\sqrt{n}`$ and $`\mu =gn=mc^2`$.

The perturbation follows the moving impurity, i.e. $`\delta \psi `$ is a function of $`(𝐫𝐕t)`$, so the coordinate derivative is related to the time derivative $`\delta \psi (𝐫𝐕t)/t=𝐕\stackrel{}{}\delta \psi (𝐫𝐕t)`$. We shall work in the frame moving with the impurity $`𝐫^{}=𝐫𝐕t`$ and the subscript over $`𝐫`$ will be dropped.

Eq. (3) for a perturbation in a homogeneous system can be conveniently solved in momentum space. In order to do this we introduce the Fourier transform of the wave function $`\delta \psi _𝐤=e^{i𝐤𝐫}\delta \psi (𝐫)d^3x`$. Eq. (3) becomes

$$\left(\mathrm{}𝐤𝐕\frac{\mathrm{}^2k^2}{2m}mc^2\right)\delta \psi _𝐤mc^2(\delta \psi _𝐤)^{}g_i\varphi _0=0.$$
(4)
Substitution of $`kk`$ and complex conjugation of (4) give the second equation. The obtained system of linear equations can be easily solved

$$\delta \psi _𝐤=g_i\varphi _0\frac{\mathrm{}𝐤𝐕+\frac{\mathrm{}^2k^2}{2m}}{(\mathrm{}𝐤𝐕)^2\frac{\mathrm{}^2k^2}{2m}\left(\frac{\mathrm{}^2k^2}{2m}+2mc^2\right)}.$$
(5)
Let us calculate the energy of the perturbation. By neglecting terms of the order of $`g_i|\delta \psi |^2`$ the energy functional (1) becomes

$$E=E_0+g_i\varphi _0^2+\mathrm{}𝐤𝐕|\delta \psi _𝐤|^2\frac{d^3k}{(2\pi )^3}+\frac{g_i\varphi _0}{2}(\delta \psi +\delta \psi ^{})_{𝐫=0}.$$
(6)
Here $`E_0=Ngn/2`$ is the energy of the system in absence of the perturbation. We expand the $`3^{rd}`$ and $`4^{th}`$ terms in powers of $`V`$. To avoid the large-$`k`$ divergency, one must also introduce a renormalization of the scattering amplitude. It is sufficient to express the coupling constant in the $`2^{nd}`$ term of eq. (6) in terms of the scattering amplitude $`b`$ using the second order Born approximation:

$$g_i=\frac{2\pi \mathrm{}^2b}{m}\left(1+\frac{2\pi \mathrm{}^2b}{m}\left(\frac{\mathrm{}^2k^2}{2m}\right)^1\frac{d^3k}{(2\pi )^3}\right).$$
(7)
Finally, by carrying out the integration over momentum space and considering $`N_{imp}`$ impurities with a concentration given by $`\chi =N_{imp}/N`$ we obtain the energy per particle

$`{\displaystyle \frac{E}{N}}`$ $`=`$ $`\left\{2\pi na^3\left(1+\chi {\displaystyle \frac{b}{a}}\right)+8\pi ^{3/2}(na^3)^{3/2}\chi \left({\displaystyle \frac{b}{a}}\right)^2\right\}{\displaystyle \frac{\mathrm{}^2}{ma^2}}`$ (9)
$`+{\displaystyle \frac{2\sqrt{\pi }}{3}}(na^3)^{1/2}\chi \left({\displaystyle \frac{b}{a}}\right)^2{\displaystyle \frac{mV^2}{2}}.`$
If we set $`V=0`$ we recover Bogoliubov’s corrections to the energy in the presence of quenched impurities . Note that even if the “mean-field” energy obtained from the GP equation in the absence of impurities ($`\chi =0`$) leaves out terms of the order of $`(na^3)^{3/2}`$, the equations we obtain in the presence of impurities in a perturbative manner still correctly describe the effect of the disorder up to the terms of the order of $`(na^3)^{3/2}`$.

If $`V0`$ a quadratic term in the impurity contribution to the energy is present. It can be denoted as $`\chi m^{}V^2/2`$ with

$$m^{}=\frac{2\sqrt{\pi }}{3}(na^3)^{1/2}\left(\frac{b}{a}\right)^2m$$
(10)
being the induced mass, i. e. the mass of particles dragged by an impurity . Applicability of the perturbation theory demands $`m^{}`$ to be small compared to $`m.`$ This gives the condition $`(na^3)^{1/2}\left(\frac{b}{a}\right)^21`$. At zero temperature the interaction between particles does not lead to depletion of the superfluid density and the suppression of the superfluidity comes only from the interaction of particles with impurities. Thus $`\left(\text{10}\right)`$ defines the normal density $`\rho _n=\chi \rho m^{}/m`$.

This result is in agreement with the one obtained by the means of Bogoliubov transformation starting from the Hamiltonian written in the second-quantized form in the presence of disorder .

The normal density of a superfluid is an observable quantity. It was evaluated in liquid $`{}_{}{}^{4}He`$ by measuring of the moment of inertia of a rotating liquid or by measuring of the second sound velocity. Both methods can be, in principle, developed for BEC gases.

### B Drag force and energy dissipation

The force with which the impurity acts on the system is

$$𝐅=|\psi \left(𝐫\right)|^2\stackrel{}{}(g_i\delta \left(𝐫\right))d^3x=g_i(\stackrel{}{}|\psi (𝐫)|^2)_{𝐫=0}.$$
(11)
Expanding the wave function into the sum of $`\varphi _0`$ and $`\delta \psi `$ and neglecting terms of order $`\delta \psi ^2`$ we obtain

$`𝐅`$ $`=`$ $`g_i\varphi _0{\displaystyle i𝐤\left[\delta \psi _𝐤+\left(\delta \psi _𝐤\right)^{}\right]\frac{d^3k}{(2\pi )^3}}`$ (12)
$`=`$ $`{\displaystyle \frac{2\left(g_i\varphi _0\right)^2i𝐤\left(\mathrm{}^2k^2/2m\right)}{(\mathrm{}𝐤𝐕+i0)^2\frac{\mathrm{}^2k^2}{2m}\left(\frac{\mathrm{}^2k^2}{2m}+2mc^2\right)}\frac{d^3k}{\left(2\pi \right)^3}},`$ (13)
where we added an infinitesimal positive imaginary part $`+i0`$ to the frequency $`𝐤𝐕`$ according to the usual Landau causality rule. The drag force is obviously directed along to the velocity $`𝐕`$.

We can carry out the integration with respect to $`\mathrm{cos}\vartheta ,`$ where $`\vartheta `$ is the angle between the momentum $`𝐤`$ and velocity $`𝐕`$, using the formula $`\frac{1}{x+i0}=𝒫\frac{1}{x}i\pi \delta (x)`$. Due to the integration between symmetric limits, only the imaginary part contributes to the integral with respect to $`\mathrm{cos}\vartheta `$. The poles in the integration over $`\mathrm{cos}\vartheta `$ appear if the square root in the denominator is smaller than one, which leads to the restriction on the values of momentum which contribute

$$|k|k_{\mathrm{max}}=2m(V^2c^2)^{1/2}/\mathrm{}.$$
(14)
Thus the energy dissipation takes place only if the impurity moves with a speed larger than the speed of sound. Integration with respect to $`k`$, taking into account restriction (14), finally gives

$$F_V=4\pi nb^2mV^2(1c^2/V^2)^2.$$
(15)
The energy dissipation, $`\dot{E}=F_VV`$, can be evaluated by measuring the heating of the gas.

For large $`V`$ the force is proportional to $`V^2`$. The energy dissipation per unit time can then be presented as $`\dot{E}=\gamma E`$ with the damping rate $`\gamma nb^2V.`$

Note in conclusion that our perturbative calculations can not describe processes involving dissipation of energy due to creation of quantized vortex rings. Such a creation is possible at $`V<c`$ but has a small probability for low velocity and for a weak point-like impurity.

## III Low dimensional systems

In this type of experiment the role of the impurity can also be played by a laser beam with small enough size and intensity. The Fourier components of the perturbed wavefunction $`\delta \psi _𝐤`$ are given by the formula (5), which is derived in an arbitrary number of dimensions. The only difference is in the substitution of $`d^3k/(2\pi )^3`$ with $`d^Dk/(2\pi )^D`$ in the integrals.

### A Two-dimensional system

There are different possible geometries of the experiment. One can create a two-dimensional perturbation in the three-dimensional condensate. Such a two-dimensional impurity can be created, analogously to the MIT experiment , by a thin laser beam. Such a beam creates a cylindrical hole in the condensate, which is stirred by moving the position of the laser beam. Another possibility is to fix the position of the laser beam along the long axis of an elongated condensate, so that the dissipation can be studied by shaking the trap and exciting the breathing modes. The problem is to create a beam with a diameter which is small with compared to the correlation length. The theory can be easily generalized for beams of finite diameter. The intensity of the beam can be tuned to satisfy the condition of a weak perturbation.

The more interesting possibility is the investigation of true two-dimensional condensates, which can be created in plane optical traps, produced by a standing light wave. If the light intensity is large enough, tunneling between planes is small and the condensates behave as independent two dimensional systems. The impurity can again be created by a laser beam perpendicular to the condensate plane. Another possibility is to use impurity atoms, which can be drive by a laser beam, with a frequency close to the atomic resonance of the impurity.

In two dimensions one has $`d^2k=kdkd\vartheta `$. The drag force $`𝐅`$ is different from zero only if the denominator has poles, which means that the velocity $`V`$ must be larger than the speed of sound $`c`$. Only momenta smaller than $`k_{max}`$ (see eq. (14)) contribute to the integral.

Integration for the velocities $`V>c`$ gives us the drag force in the two-dimensional case:

$$F_V^{2D}=g_i^2n_2m^2(V^2c^2)/(\mathrm{}^3V),$$
(16)
where $`n_2`$ is two-dimensional density.

In a quasi two-dimensional system, i.e. when the gas is confined in the $`z`$-direction by the harmonic potential $`\left(1/2\right)m\omega _z^2r^2`$, the two-dimensional coupling constant equals $`g_i=\sqrt{2\pi }\frac{\mathrm{}^2b}{ma_z}`$, where $`a_z=\sqrt{\mathrm{}/m\omega _z}`$ is the oscillator length and $`b`$ is the three-dimensional scattering length. (We consider here only the mean-field $`2D`$ situation. See , § 17 for a more detail discussion.)

The calculation of the effective mass gives

$$m^{}=g_i^2n_2m/(4\pi ^2\mathrm{}^2c^2).$$
(17)
Notice again that our calculations do not take into account creation of vortex pairs which is possible at $`V<c`$.

### B One-dimensional system. Mean-field theory

In one dimension the integration is straightforward. The integration $`(\text{13})`$ over $`k`$ gives $`2\pi i`$ if $`V>c`$ and zero otherwise. So, the force is $`F_V^{1D}=2g_i^2n_1m/\mathrm{}^2`$, where $`n_1`$ is the linear density. In a quasi one dimensional system (i.e. a very elongated trap or a waveguide) there are no excitations in the radial harmonic confinement and the coupling constant is given by $`g_i=\frac{\mathrm{}^2}{mb_{1D}}`$ with $`b_{1D}=a_{}^2/b`$, where $`a_{}=\sqrt{\mathrm{}/m\omega _{}}`$ and one has $`F^{1D}=2n_1\mathrm{}^2/mb_{1D}^2`$.

An interesting peculiarity is that the result does not depend on the velocity $`V`$ (where, of course, the velocity must be larger than the speed of sound). The same result for the $`\delta `$-potential was found in . Calculation of the effective mass gives $`m^{}=g_i^2n_1/2\mathrm{}c^3`$.

In a 1D system energy dissipation is possible at $`V<c`$ due to creation of the “gray solitons” first considered in . Non-linear calculations  show that the critical velocity for this process decreases with increasing coupling constant $`g_i`$.

This theory can be checked in an experiment in a three-dimensional condensate. The two-dimensional impurity can be presented by a moving light sheet.

### C One-dimensional system. Bethe-ansatz theory

We saw in the previous subsection that for a weakly interacting impurity the drag force appears only when the impurity velocity $`V`$ is larger that the Landau critical velocity, which is equal to the velocity of sound $`c`$. The situation is, however, different in the Bethe-ansatz Lieb-Liniger theory of a 1D Bose gas . According to this theory excitations in the system actually have a fermionic nature. Even a low frequency perturbation can create a particle-hole pair with a total momentum near $`2p_F2\mathrm{}k_F=\mathrm{}2\pi n_1`$. To calculate the drag force for this case we will use the dynamic form factor of the system $`\sigma (\omega ,k)`$ (we follow notation of , § 87). The dissipated energy at $`T=0`$ can be calculated as

$$\dot{E}=\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{dk}{2\pi }\underset{0}{\overset{\mathrm{}}{}}\frac{d\omega }{\pi }\omega \frac{n_1}{2\mathrm{}}\sigma (\omega ,k)\left|U(\omega ,k)\right|^2,$$
(18)
where $`U(\omega ,k)=2\pi g_i\delta (\omega kV)`$ is the Fourier transform of the impurity potential $`U(t,z)=g_i\delta (zVt)`$. One has $`\left|U(\omega ,k)\right|^2=2\pi g_i^2t\delta (\omega kV)`$, where $`t`$ is ”time of observation”. Thus the energy dissipation per unit of time is

$$\dot{E}=F_VV=\frac{g_i^2n_1V}{\mathrm{}}\underset{0}{\overset{\mathrm{}}{}}\frac{dk}{2\pi }k\sigma (kV,k),$$
(19)
where $`F_V`$ is the drag force. We will try to estimate the velocity dependence of $`F_V`$.

For low frequency dissipation the important values of $`k`$ are near $`2k_F`$. According to

$$\sigma (\omega ,2k_F)\omega ^{\left(\eta 2\right)},\omega 0,$$
(20)
where $`\eta =\frac{2\mathrm{}k_F}{mc}=\frac{2\pi \mathrm{}n_1}{mc}2`$ is the characteristic parameter of a 1D Bose gas. In the mean-field limit when $`n_1\mathrm{}`$ the parameter $`\eta \mathrm{}`$. In the opposite case of a small density bosons behave as impenetrable particles (Girardeau limit ) and the dynamic form-factor coincides with the one of an ideal Fermi gas. In this limit $`\eta =2.`$

In the general case one can calculate $`\sigma (\omega ,k)`$ at small $`\omega `$ and $`k2k_F`$generalizing the method of Haldane  for the case of time-dependent correlation functions. Calculations give

$$\sigma (\omega ,k)=\frac{n_1c}{\omega ^2}\left(\frac{\mathrm{}\omega }{mc^2}\right)^\eta f\left(\frac{c\mathrm{\Delta }k}{\omega }\right),\omega >0,k>0,$$
(21)
where $`k=2k_F+\mathrm{\Delta }k`$ and the function $`f(x)`$ is

$$f\left(x\right)=A\left(\eta \right)\left(1x^2\right)^{\eta /21}$$
(22)
in the interval $`\left|x\right|<1`$ and is equal to zero at $`\left|x\right|1`$ (see also ). The constant $`A\left(\eta \right)`$ can be calculated in two limiting cases: $`A\left(\eta =2\right)=\pi /4`$ (see  § 17.3) and $`A\left(\eta \right)8\pi ^2/\left[\left(8C\right)^\eta \mathrm{\Gamma }^2\left(\frac{\eta }{2}\right)\right]`$ , where $`C=1.78`$… is the Euler ’s constant, for $`\eta 1`$ (details of the calculation will be published elsewhere).

Substituting $`\left(\text{21}\right)`$ into $`\left(\text{19}\right)`$ we finally find velocity dependence of the drag force:

$$F_V=\frac{\mathrm{\Gamma }\left(\frac{\eta }{2}\right)}{2\sqrt{\pi }\mathrm{\Gamma }\left(\frac{\eta +1}{2}\right)}A\left(\eta \right)\frac{g_i^2n_1^2}{\mathrm{}V}\left(\eta \frac{V}{c}\right)^\eta .$$
(23)
Equation $`\left(\text{23}\right)`$ is valid for the condition $`Vc`$.

Thus in the Girardeau strong-interaction limit $`F_VV`$ and Bose gas behaves, from the point of view of friction, as a normal system, where the drag force is proportional to the velocity. On the contrary, in the mean-field limit the force is very small and the behavior of the system is analogous to a 3D superfluid. However, even in this limit the presence of the small force makes a great difference. Let us imagine that our system is twisted into a ring, and that the impurity rotates around the ring with a small angular velocity. If the system is superfluid in the usual sense of the word, the superfluid part must stay at rest. Presence of the drag force means that equilibrium will be reached only when the gas as a whole rotates with the same angular velocity. From this point of view the superfluid part of the 1D Bose gas is equal to zero even at $`T=0.`$ Notice that in an earlier paper  the author concluded that $`\rho _s=\rho `$ at $`T=0`$ for arbitrary $`\eta `$. We believe that this difference results from different definitions of $`\rho _s`$ and reflects the non-standard nature of the system.

Equation $`\left(\text{23}\right)`$ is equivalent to a result which was obtained by a different method in , with a model consisting of an impurity considered as a Josephson junction. Notice that the process of dissipation, which in the language of fermionic excitations can be described as creation of a particle-hole pair, corresponds in the mean-field limit to creation of a phonon and a small-energy soliton. It seems that such a process cannot be described in the mean-field approach in the linear approximation.

Experimental confirmation of these quite non-trivial predictions demands a true one-dimensional condensate, where non mean-field effects can be sufficiently large. Such condensates have been investigated for the first time in experiments . In experiments  condensates have been created in the form of elongated independent “needles” in optical traps, consisting of two perpendicular standing laser waves. The role of an impurity in this case must be played by a light sheet, perpendicular to the axis of condensates and moving along them.

Notice also, that application of the additional light waves in this experiments of this type allows one to create a harmonic perturbation of the form

$$U(t,z)=U_0\mathrm{cos}\left(\omega tkz\right),k=2k_F+\mathrm{\Delta }k$$
(24)
with small $`\omega `$ and $`\mathrm{\Delta }k`$. Such potential with was used in  for experiments with 1D condensate in a periodic lattice. However, for a small amplitude $`U_0`$, measurement of the dissipation energy $`Q`$ gives, according to $`\left(\text{18}\right),`$ the dynamic form-factor $`S(\omega ,k)`$ directly.

## IV Conclusions

We have studied motion of an impurity through the condensate at zero temperature by considering the perturbation of a stationary solution of the GP equation. We calculated the induced mass which contributes to the mass of normal component. We find that the motion at small velocities is dissipationless in one-, two-, and three- dimensional systems, although movement with velocities larger than the speed of sound leads to a non-zero drag force due to Cherenkov radiation of phonons. The expressions for the drag force are calculated. We used results for the dynamic form factor of exact Lieb-Liniger theory to investigate the velocity dependence of the drag force in a 1D system. The form factor was calculated with the help of the Haldane method of calculations of correlation functions. The drag force exists at an arbitrarily small velocity of motion, but is very small in the mean-field limit. The dynamic form-factor can be also directly measured by applying a harmonic time-dependent perturbation on one-dimensional condensates .

Acknowledgments: We thank T. Esslinger, F. D. M. Haldane, B. Jackson, Yu. Kagan, V. Korepin, E. Sonin and S. Stringari for useful discussions.


text:
# Probing dark energy perturbations: the dark energy equation of state and speed of sound as measured by WMAP

## I Introduction

With the recent unveiling of the Wilkinson Microwave Anisotropy Probe (WMAP) results, measuring the Cosmic Microwave Background (CMB) anisotropy Bennett03 , the on-going supernovae searches super  and the upcoming completion of the Sloan Digital Sky Survey, amongst others, we are seeing a wealth of precision observational data being made available. To a great extent the standard $`\mathrm{\Lambda }`$-CDM scenario fits the data well Spergel03 .

However, the WMAP data might suggest that some modifications to the standard scenario are needed. One possible hint at required modifications is the deficit of large scale power in the temperature map, and in particular, the low CMB quadrupole whose posterior probability is less than a few hundredth (see e.g. for possible interpretations Efstathiou03 ; Contaldi03  and Efstathiou03b  for a discussion of this number). One possibility is that this lack of large scale power might point to some particular properties of the dark energy. The dominant contribution to fluctuations on these scales is the Integrated Sachs Wolfe (ISW) effect which describes the fluctuations induced by the passage of CMB photons through the time evolving gravitational potential associated to nearby ($`z<5`$) large scale structures (LSS). One property we expect of dark energy is that it suppresses the gravitational collapse of matter at relatively recent times, which in turn suppresses the gravitational potential felt by the photons, thereby leaving a signature in the ISW correlations. Since this signature is created by the time evolving potential associated with relatively close LSS, it should be intimately correlated with any tracer of LSS Kofman85 ; Crittenden96 . A positive detection of such a cross correlation using WMAP data, assuming a cosmological constant as the dark energy, has recently been measured Boughn03 ; Nolta03 ; Fosalba03 .

However the underlying cause of the dark energy is still unknown; and such observational inferences offer rich prospects for guiding and leading the theoretical effort. A wide variety of models have been proposed to explain observations, from the unperturbed cosmological constant, to a multitude of scalar field quintessence and exotic particle theories (see Peebles02  for a review).

Much effort has been put into determining the equation of state of dark energy, in an attempt to constrain and direct theories. Since the equation of state affects both the background expansion and the evolution of matter perturbations there are a wealth of complementary observations available (again see Peebles02  and references therein).

An equally insightful, but less investigated, characteristic of dark energy is the speed of sound within it. This does not affect the background evolution but is fundamental in determining a dark energy s clustering properties, through the Jeans scale. It will, therefore, have an effect on the evolution of fluctuations in the matter distribution.

Following the papers laying the foundations for cosmological perturbation theory Bardeen80 ; Kodama84 , the effect of the speed of sound on observables was considered in more detail: for the CMB and large scale structure Hu98 ; Hu99  and in the context of lensing Hu02  . Observational implications of the speed of sound in a variety of dark energy models have also recently been discussed: for example for k-essence Erikson01 ; Dedeo03 , condensation of dark matter Bassett02  and the Chaplygin gas, in terms of the matter power spectrum Sandvik02 ; Beca03 ; Reiss03  and combined full CMB and large scale structure measurements Bean03 ; Amendola03 .

Minimally coupled scalar field, quintessence models, commonly have a non-adiabatic speed of sound close or equal to unity (in units of c, the speed of light),(see for example Ferreira97 ). By contrast however, the adiabatic Chaplygin gas model,( e.g. motivated by a rolling tachyon Gibbons02 ) has a speed of sound directly proportional to the equation of state, both of which are approximately zero up until late times when the dark energy starts to dominate. It is conceivable therefore that distinctions between such models might be able to be made through the detection of a signature of the dark energy speed of sound: in the large scale ISW correlations, and in the cross correlation of the CMB with the distribution of large scale structure Crittenden96 .

In section II we briefly review parameters describing a general fluid and the issues that arise when establishing a fluid’s speed of sound. In section III we describe the implications of equation of state and speed of sound on perturbation evolution in the fluid and CDM. We consider a toy model with slowly varying equation of state and sound speed applicable in a wide variety of minimally coupled scalar field theories. In section IV we discuss the potential for detecting the speed of sound using late time perturbation evolution, in the ISW effect, through the autocorrelation of the CMB temperature power spectrum. In section V we extend the analysis to the cross correlation of the WMAP CMB data with distribution. Finally, in section VI we summarize our findings.

## II The speed of sound within general matter

For a perfect fluid the speed of sound purely arises from adiabatic perturbations in the pressure, $`p`$, and energy density $`\rho `$ and the adiabatic speed of sound, $`c_a^2`$ is purely determined by the equation of state $`w`$,

$`w_i`$ $``$ $`{\displaystyle \frac{p_i}{\rho _i}}`$ (1)
$`c_{ai}^2`$ $``$ $`{\displaystyle \frac{\dot{p}_i}{\dot{\rho }_i}}=w_i{\displaystyle \frac{\dot{w}_i}{3(1+w_i)}}`$ (2)
where the subscript $`i`$ denotes a general specie of matter, where dots represent derivatives with respect to conformal time and where $``$ is the Hubble constant with respect to conformal time.

In imperfect fluids, for example most scalar field or quintessence models, however, dissipative processes generate entropic perturbations in the fluid and this simple relation between background and the speed of sound breaks down and we have the more general relation

$`c_{si}^2`$ $``$ $`{\displaystyle \frac{\delta p_i}{\delta \rho _i}}.`$ (3)
In order to establish the speed of sound in these cases we must look to the full action for the fluid described often through the form of an effective potential. In this case, the speed of sound can be written in terms of the contribution of the adiabatic component and an additional entropy perturbation $`\mathrm{\Gamma }_i`$ and the density fluctuation in the given frame $`\delta _i`$ Kodama84 ,

$`w_i\mathrm{\Gamma }_i`$ $``$ $`(c_{si}^2c_{ai}^2)\delta _i`$ (4)
$`=`$ $`{\displaystyle \frac{\dot{p}_i}{\rho _i}}\left({\displaystyle \frac{\delta p_i}{\dot{p}_i}}{\displaystyle \frac{\delta \rho _i}{\dot{\rho }_i}}\right)`$ (5)
$`\mathrm{\Gamma }_i`$ is the intrinsic entropy perturbation of the matter component, representing the displacement between hypersurfaces of uniform pressure and uniform energy density. In this paper, we are solely interested in probing the intrinsic entropy of the dark energy component. It is worth noting that in a multi-fluid scenario, in addition to the intrinsic entropy perturbations denoted by $`\mathrm{\Gamma }_i`$, further contributions to the total entropy perturbation of the system can arise from the relative evolution of two or more fluids with different adiabatic sound speeds, and through non-minimal coupling (see for example Balakin03 ).

Whereas the adiabatic speed of sound, $`c_{ai}`$, and $`\mathrm{\Gamma }_i`$ are scale independent, gauge invariant quantities, $`c_{si}`$ can be neither. As such the general speed of sound is gauge and scale dependent and issues of preferred frame arise. Looking at equation (5), since the fluid $`i`$ rest frame is the only frame in which $`\delta _i`$ is a gauge invariant quantity, this is the only frame in which a matter component’s speed of sound is also gauge-invariant.

A useful transformation Kodama84  relates the gauge-invariant, rest frame density perturbation, $`\widehat{\delta }_i`$, to the density and velocity perturbations in a random frame, $`\delta _i`$ and $`\theta _i`$,

$`\widehat{\delta }_i=\delta _i+3(1+w_i){\displaystyle \frac{\theta _i}{k^2}}`$ (6)
where we assume that the component is minimally coupled to other matter species and henceforth dark energy rest frame quantities are denoted using a circumflex (^) .

Using equations (5) and (6) we can rewrite the pressure perturbation in a general frame, $`\delta p_i`$, in terms of the rest frame speed of sound,

$`\delta p_i`$ $`=`$ $`\widehat{c}_{si}^2\delta \rho _i+3(1+w_i)(\widehat{c}_{si}^2c_{ai}^2)\rho _i{\displaystyle \frac{\theta _i}{k^2}}.`$ (7)
## III The speed of sound and perturbation evolution

Herein, we use the synchronous gauge and follow the notation of Ma95 . CDM rest frame quantities are denoted $`\delta _i`$ and $`\theta _i`$, while dark energy rest frame quantities use the circumflex (^). The two are related by equation (6).

The energy density and velocity perturbation evolution of a general matter component in the CDM rest frame is given by

$`\dot{\delta }`$ $`=`$ $`(1+w)\left\{\left[k^2+9^2(\widehat{c}_s^2c_a^2)\right]{\displaystyle \frac{\theta }{k^2}}+{\displaystyle \frac{\dot{h}}{2}}\right\}`$ (8)
$`3(\widehat{c}_s^2w)\delta `$
$`{\displaystyle \frac{\dot{\theta }}{k^2}}`$ $`=`$ $`(13\widehat{c}_s^2){\displaystyle \frac{\theta }{k^2}}+{\displaystyle \frac{\widehat{c}_s^2}{1+w}}\delta .`$ (9)
This set of equations illustrate clearly that linear perturbations can be fully characterised by two numbers (and their potential time evolution): the equation of state and the rest frame speed of sound.

Let us now consider a toy model with a general fluid in which the time variation in $`w`$ and $`\widehat{c}_s^2`$ is small in comparison to the expansion rate of the universe so that we can model it with constant $`w`$ (i.e. $`c_a^2w`$) and $`\widehat{c}_s^2`$. Such models are not impractical and can be used as the basis for comparison with scalar field theories such as those with scaling potentials and Chaplygin gases during the radiation and matter dominated eras.

In the matter dominated era, ignoring baryons for simplicity, CDM density perturbations are affected by the speed of sound of dark energy (denoted ‘$`x`$’) through the relation,

$`\ddot{\delta }_c+\dot{\delta }_c{\displaystyle \frac{3^2}{2}}\mathrm{\Omega }_c\delta _c=`$ (10)
$`{\displaystyle \frac{3^2\mathrm{\Omega }_x}{2}}\left[(1+3\widehat{c}_s^2)\delta _x+9(1+w)(\widehat{c}_s^2w){\displaystyle \frac{\theta _x}{k^2}}\right]`$
In the radiation and matter dominated eras the expansion rate obeys an effective power law $`a\tau ^m`$ and we find that the evolution equations admit a solution of the form $`\delta _c\delta _x/(1+w)\tau ^2`$, and $`\theta _x\tau ^3`$.

Equation (8) shows how $`w`$ and $`\widehat{c}_s^2`$ affect the relative size of dark energy and CDM perturbations. In Fig. 1 we see that, as one expects, increasing the speed of sound accentuates this suppression produced by reducing $`w`$, lowering $`\delta _x/\delta _c`$. Note that the dark energy density perturbation is well defined (remaining positive) in the rest frame, while the transformation into the CDM rest frame can make $`\delta _x`$ negative; this is just a foible of the frame one chooses, however.

The presence of dark energy perturbations leaves a $`w`$ and $`c_s^2`$ dependent signature in the ISW source term. This can be written in terms of the time variation of the anisotropic stress and the rest frame density perturbations of each matter component.

$`\dot{\mathrm{\Phi }}\dot{\mathrm{\Psi }}`$ $`\stackrel{z0}{}`$ $`{\displaystyle \frac{1}{k^2}}{\displaystyle \frac{d}{d\tau }}\left[^2\left(\mathrm{\Omega }_c\delta _c+\mathrm{\Omega }_x\widehat{\delta }_x\right)\right]`$ (11)
Since $`\widehat{\delta }_x`$ is suppressed in comparison to $`\delta _c`$, the dominant contribution to the ISW will come from the CDM perturbations. Subsequently it will be suppression of these (in comparison to a $`\mathrm{\Omega }_c=1`$ scenario) through the effect of the dark energy speed of sound and equation of state that will leave a signature in the ISW. Fig 2. shows how as one increases $`\widehat{c}_s^2`$ the ISW effect increases.

## IV Constraints using WMAP temperature fluctuations at large scales

In this section, we investigate the joint constraints on the equation of state and the speed of sound that can be infered using the CMB temperature power spectrum. As was discussed earlier, the main effect of a sound speed smaller than the speed of light will be felt at late times and large scales and will thus only affect the very large scale CMB temperature fluctuations arising from the late ISW effect, for which WMAP already provide us will full sky cosmic variance limited measurement Bennett03 . The Fourier component of the fluctuations arising from the ISW effect are given by (for this equation only, we ignore the well known issues related to the sphericity of the observed sky).

$`{\displaystyle \frac{\delta T}{T_0}}(z,k)`$ $`=`$ $`{\displaystyle _{\eta (z)}^{\eta _0}}2\dot{\mathrm{\Phi }}(\eta ,k)𝑑\eta `$ (12)
$`=`$ $`{\displaystyle \frac{3H_0^2\mathrm{\Omega }_m^0}{c^2k^2}}\delta _c^0(k){\displaystyle _0^z}{\displaystyle \frac{dg}{dz}}𝑑z`$ (13)
where $`c^2`$ is the square of the speed of light, $`H_0`$ is the Hubble constant today, $`\delta _c^0(k)=\delta _c(z=0,k)`$ and $`\mathrm{\Omega }_m^0`$ is the fractional energy density in matter (CDM+baryons) today and $`g(z,k)=(1+z)D(z,k)`$ where $`D(z,k)`$ is the linear growth factor given by $`\delta _c(z,k)=D(z,k)\delta _c^0(k)`$. In the limit that $`w`$ tends to -1 $`D(z,k)`$ is scale independent and can be approximated by

$`D(z)`$ $`=`$ $`{\displaystyle \frac{5\mathrm{\Omega }_c^0E(z)}{2}}{\displaystyle _z^{\mathrm{}}}{\displaystyle \frac{(1+z)dz}{E(z)^3}}`$ (14)
$`E(z)`$ $`=`$ $`{\displaystyle \frac{H(z)}{H_0}}=\left[{\displaystyle \underset{i}{}}\mathrm{\Omega }_i^0(1+z)^{3(1+w_i)}\right]^{\frac{1}{2}}.`$ (15)
however this approximation does not provide the degree of precision that is required for $`w>1`$, even in the absence of dark energy perturbations Linder03 . Because of this and in order to factor in the late time effect of dark energy perturbations and their scale dependence, for $`\widehat{c}_s^20`$, we explicitly calculate the linear growth function, $`D(z,k;\widehat{c}_s^2,w)`$ for each model.

Note that the effect of the speed of sound comes in solely through the value of $`\delta _c^0`$ while the equation of state affects both $`\delta _c^0`$ and the linear growth factor.

The associated auto power spectrum is given subsequently by

$`C_{\mathrm{}}^{TT}`$ $`=`$ $`4\pi {\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dk}{k}}{\displaystyle \frac{\delta T}{T_0}}(k){\displaystyle \frac{\delta T}{T_0}}(k^{})j_{\mathrm{}}^2(k\eta ).`$ (16)
In order to probe solely the effect of $`\widehat{c}_s^2`$ on $`C_{\mathrm{}}^{TT}`$, we will compare to WMAP observations a family of models lying along the angular diameter degeneracy surface present in the CMB spectrum. To do so, we keep $`\mathrm{\Omega }_ch^2=0.135`$ and $`\mathrm{\Omega }_bh^2=0.0224`$ to be consistent with the WMAP best fit Spergel03  and choose $`w`$ and $`h`$ such that the angular diameter distance to last scattering is the same. Other parameters correspond to the best fit model of Spergel03  (table 7). Doing so, we ensure that only the large scale correlations vary with each model and that in all other respects they fit the WMAP data well. Note that we also have vary slightly the overall amplitude due to the change in the first peak height ISW plateau ratio. We consider $`w`$ values between 0 and -1 and $`\widehat{c}_s^2`$ between 0 and 1. Given this grid of model, we can then deduce easily the likelihood of the data using the publically available code provided by the WMAP team WMAPanalysis03 , from which we can deduce some joint constraints on $`w`$ and $`\widehat{c}_s^2`$.

In Fig. 3 we show the variation of the CMB TT power spectrum as one varies $`\widehat{c}_s^2`$ from 0 to 1 for a model with $`w=0.3`$ and -0.9. One can see that increasing $`\widehat{c}_s^2`$ increases the suppression of the CDM perturbations and therefore increases the power on large scales. The effect decreases though as one decreases $`w`$; at low $`w`$, the suppression due to the equation of state itself will generate a dominant ISW effect on top of which a subdominant contribution from $`\widehat{c}_s^2`$ is then superimposed. Those results agree with the one obtained in Hu02 .

In Fig. 4 we show the likelihood plot from the WMAP data in the $`w\widehat{c}_s^2`$ plane. The low quadrupole, and other low $`\mathrm{}`$ $`C_{\mathrm{}}`$’s lead to a value of $`\widehat{c}_s^2<0.04`$ being preferred by the data, at the 1$`\sigma `$ level, although as one moves to lower $`w`$ the ability to distinguish between different values of the speed of sound disappears, because of cosmic variance.

The cosmic variance thus limit our ability to constraint the dark energy speed of sound using temperature $`C_{\mathrm{}}`$ only. However, given the fact that all the constraints comes from the ISW effect, it is natural to consider the cross correlation the CMB with the large scale distribution of matter near us, correlation that is a direct probe of the late ISW. In theory then, this might give us a better and different probe into $`\widehat{c}_s^2`$, so that both should be combine eventually. We consider this in the next section.

## V Constraints using CMB and large scale structure cross correlation

As stated earlier, the dark energy affects very large scale modes of dark matter density perturbations. As shown on Fig. 5 , those modes are outside the range of current wide field galaxy surveys. For example, the SDSS measured the galaxy power spectrum down to $``$ 0.01 $`h`$Mpc<sup>-1</sup> ”only” (see e.g. Dodelson01 ). Full sky survey exists though, but their particular properties and intrinsic limitations restricted their use as direct probes of the matter power spectrum at those scales. For example, the NRAO VLA Sky Survey (NVSS) Condon98  encompass such a wide variety of objects that the difficulties in modeling the biases at stake prevented its usage to directly measure (dark) matter density fluctuations at those scales and infer this way any precise cosmological constraint. However, their use in conjunction with large scale CMB fluctuation measurements allows us to circumvent somehow this difficulty.

Indeed, within a given cosmological model, one can look at the surveyed objects as a simple linearly bias tracer of dark matter perturbations, a reasonnable approximation on those very large scales. By measuring the auto-correlation function (ACF) of those objects on those scales, one can infer the model dependant effective bias for this composite population. Since this population trace the large scale gravitational potential, it should correlate with the CMB fluctuations induced by the same potential through the ISW effect. The angular dependance of this cross-correlation function (CCF) and its amplitude both depend on the tracer properties (bias and redshift distribution) and on the particular cosmological model considered. In particular, we would expect an important dependacy on the dark energy properties which drive the evolution of the universe at those late times and large scales. Modeling the tracer properties, we can thus in principle constrain the cosmology. This has been advocated first in Crittenden96 , studied in details in Peiris00 ; Cooray02  and performed effectively using as a tracer the NVSS sources Boughn02b ; Boughn03 ; Nolta03 , HEAO-1X-ray sources Boughn98 ; Boughn02  or APM galaxies Fosalba03 .

So far, this correlation has been probed to prove the very existence of dark energy and to constrain its overall density. We here extend this approach and try to investigate the potential constrains on its very perturbative properties, i.e. jointlly its equation of state and its sound speed. We will use as a data-set the ACF and CCF measurements of Nolta03  performed using the NVSS catalog and the WMAP 1-year maps. For the sake of simplicity we will follow the same notations that we recall brievly.

The measurements of fluctuations in the nearby matter distribution from measuring the radio sources distribution can be expressed in terms of the fractional source count perturbation given by,

$`{\displaystyle \frac{\delta N}{N_0}}(\widehat{n})`$ $`=`$ $`b_r\delta _c^0(\widehat{n}){\displaystyle _0^z}{\displaystyle \frac{d\stackrel{~}{N}}{dz}}D(z,\widehat{n})𝑑z`$ (17)
where $`b_r`$ is the linear bias in the matter distribution, $`N_0`$ is the mean source count per pixel (147.9 for 1.8$`\mathrm{deg}`$ square pixels used in Nolta03 ), and $`d\stackrel{~}{N}/dz`$ is the normalised redshift distribution of galaxies, such that $`(d\stackrel{~}{N}/dz)𝑑z=1`$. For the latter we adopt the model of Dunlop90 .

The dimensionless two point correlation function between two quantities $`X`$ and $`Y`$ with background values $`X_0`$ and $`Y_0`$ in positions $`\widehat{n}`$ and $`\widehat{n}^{}`$ in the sky is given by

$`C^{XY}(\theta )`$ $`=`$ $`{\displaystyle \frac{\delta X}{X_0}}(\widehat{n}){\displaystyle \frac{\delta Y}{Y_0}}(\widehat{n}^{})_{\widehat{n}\widehat{n}^{}=cos\theta }`$ (18)
$`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle \underset{\mathrm{}=2}{\overset{\mathrm{}}{}}}(2l+1)C_{\mathrm{}}^{XY}P_{\mathrm{}}(cos\theta )`$
For the fractional source count and CMB temperature cross and auto correlations

$`C_{\mathrm{}}^{NT}`$ $`=`$ $`4\pi {\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dk}{k}}{\displaystyle \frac{\delta N}{N_0}}(k){\displaystyle \frac{\delta T}{T_0}}(k^{})j_{\mathrm{}}^2(k\eta )`$ (19)
$`=`$ $`4\pi {\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dk}{k}}\mathrm{\Delta }^2(k)f_{\mathrm{}}^N(k)f_{\mathrm{}}^T(k),`$
$`C_{\mathrm{}}^{NN}`$ $`=`$ $`4\pi {\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dk}{k}}{\displaystyle \frac{\delta N}{N_0}}(k){\displaystyle \frac{\delta N}{N_0}}(k^{})j_{\mathrm{}}^2(k\eta )`$ (20)
$`=`$ $`4\pi {\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dk}{k}}\mathrm{\Delta }^2(k)f_{\mathrm{}}^N(k)f_{\mathrm{}}^N(k),`$
where

$`\mathrm{\Delta }^2(k)=\delta _c^0(k)\delta _c^0(k^{})=\delta (kk^{})k^3P(k)/2\pi ^2`$ (21)
and where the filter functions for the temperature and number count fluctuations respectively are given by

$`f_{\mathrm{}}^T(k)`$ $`=`$ $`{\displaystyle \frac{3H_0^2\mathrm{\Omega }_m^0}{c^2k^2}}{\displaystyle _0^z}{\displaystyle \frac{dg(z,k)}{dz}}j_{\mathrm{}}(k\eta (z))𝑑z`$ (22)
$`f_{\mathrm{}}^N(k)`$ $`=`$ $`b_r{\displaystyle _0^z}{\displaystyle \frac{d\stackrel{~}{N}}{dz}}D(z,k)j_{\mathrm{}}(k\eta (z))𝑑z,`$ (23)
where $`D(z,k)`$ has been previously defined in section IV.

The results of this approach are summarised in figures 6 and 7. As was used in section IV, we consider here a family of models lying along the CMB angular diameter degeneracy surface, and move along it by varying simultaneoulsy $`w`$ and $`h`$. For each of these background, we then consider various $`c_s^2`$ and compute the linear transfer function using a modified version of the CMBfast cmbfast  or CAMB camb  softwares . For a sample of those models, we plotted both the predicted ACF, $`C^{NN}(\theta )`$, from which we infer the bias (figure 6) and the subsequent predicted CCF, $`C^{NT}(\theta )`$ which one can compare with the measurements of Nolta03  (figure 7). Obviously, even if this signal does indeed have some dependance with regards to the dark energy pertubations parameters, $`w`$ and $`c_s^2`$, most of the effect is absorbed in the bias measurement illustrated in figure 6. The uncertainty in the bias is the main hindrance in placing constraints on both the equation of state and speed of sound, so that current data do not allow this correlation is not measured well enough to constrain our models. We obtain a bias range and a bias evolution similar as the one obtained by Nolta03 , i.e. $`b_r1.4`$ for our fiducial $`\mathrm{\Lambda }`$CDM model, and $`b_r`$ tends to decrease when $`w`$ increases (in range comprised between 1.4 and 2.2). Note that the plotted error bars are heavily correlated. The knowledge of the correlation matrix computed in Nolta03  allows us however to compute a well defined $`\chi ^2`$ statistic. Note also that those error bars include cosmic variance only but do not take into account the shot noise in the radio source catalog. It would be important to consider the uncertainties in the bias measurement that propagates directly in the signal prediction but we leave this issue for future work. Nevertheless, given the fact that all our models lie within the one sigma error bars, we will not illustrate here by a contour plot those weak joint constraints on $`w`$ and $`c_s^2`$.

A natural and important question that arises at this level concerns the future prospect for the measurements of this correlation, depending on the used LSS tracers as well as the survey considered. Although some studies have already been done Hu99 ; Peiris00 ; Cooray02 , more specific investigations are necessary. In particular an independent measurement of the bias, along with improved scale and depth of survey will all contribute to vastly improving on the current observational uncertainties.

## VI Conclusions

We have reviewed the effect of the speed of sound of dark energy on CDM and dark energy perturbations. While a positive dark energy speed of sound suppresses the CDM perturbations, it is the deviation from adiabaticity, in combination with the equation of state that determine the degree of suppression of the amplitude dark energy perturbation in comparison to those of the CDM.

We have found the CMB large scale temperature fluctuations, dominated by the ISW effect, are a promising tool to measure the speed of sound. The suppression of CDM matter perturbations drive the late time ISW effect.

From the auto correlation of the WMAP data with itself we obtain a $`1\sigma `$ ”constraint” on the speed of sound $`\widehat{c}_s^2<0.04`$, using scenarios that minimise contributions to the likelihood on small scales (from the peaks) as much as possible by using well known degeneracies to follow the WMAP best fit model as closely as possible. The main limitation in obtaining constraints from the auto correlation is the cosmic variance.

We have also investigated the cross correlation of the large scale CMB with fluctuations in the nearby mass distribution using the NVSS radio source catalogue. We here again find that cosmic variance is a strong limitation and prevent us from placing any strong constraint in the $`wc_s^2`$ plane.

However, since the potential of such an analysis might be unique in unveiling the mysteries of the dark energy, it is important to explore further out the prospect of future potential large scale probe of the gravitational portential and so of the ISW (LSST, PLANCK, CMBPOL). We have presented some estimates of prospective constraints that one might obtain from cross correlation of large scale probes with CMB however we leave this exploration for future work.

Acknowledgements We would like to thank Robert Caldwell, Sean Carroll, Anthony Challinor, Rob Crittenden, Joe Hennawi, Mike Nolta, Hiranya Peiris, Martin White and, especially, David Spergel for very helpful discussions and questions in the course of this work. O.D. acknowledges the Aspen Center for Physics where part of this work was pursued. R.B. and O.D. are supported by WMAP and NASA ATP grant NAG5-7154 respectively.