text:
# Revisiting the Cosmic Cooling Crisis

## 1 Introduction

Gas cooling in the expanding Universe is an intrinsically unstable process because cooling acts to increase the density of the gas, which in turn increases the cooling rate. The consequence of this is that as soon as the gas within a system is able to cool at all it tends to do so catastrophically, only regulated by the speed at which the gas can respond to the new configuration. Systems which collapse at low redshift have a low mean density and, thus, long cooling times which generally exceed their dynamical times. Systems collapsing at higher redshifts have a higher mean density and, because the cooling rate depends strongly upon the gas density, they cool much more rapidly. In a Universe dominated by cold dark matter (CDM), the hierarchical growth of structure consequently results in almost all of the baryons cooling by the present day, unless additional physics is considered \[Cole 1991, White & Frenk 1991, Blanchard et al.  1992\].

Unfortunately, the naive theoretical prediction that all the baryonic material in the Universe should have cooled into galaxies and formed stars conflicts with the observation of X-ray emission from galaxy clusters. This emission demonstrates that large amounts of hot gas persist in the Universe. This material can arise from two possible sources: either it was never contained within a collapsed halo in which cooling was efficient, or energy injection (for instance feedback of energy due to supernovae) has either reheated or prevented the cooling of the gas. The commonly used theory of Press & Schecter \[Press & Schechter 1974\] (see also Bower 1991) implicitly assumes that all the matter in the Universe is contained within haloes. The best numerical models \[Jenkins et al.  2000\] indicate that the fraction of material contained within haloes is certainly high, and only leave room for a small uncollapsed fraction. Therefore, we appeal to feedback mechanisms \[Larson 1974, White & Rees 1978, Cole 1991, White & Frenk 1991\] to reduce the amount of cold gas within collapsed haloes.

In this paper we review the observational constraints on the density of cold and hot baryons, both globally and in rich clusters (§2). These data establish the global cooled gas fraction at $`f_{\mathrm{c},\mathrm{global}}`$$`0.073h`$, and there is no convincing evidence that it is much higher in clusters. We then examine cosmological simulation results in §3, and show that they can easily cool far too much gas to be consistent with these data. Finally, in §4 we discuss the necessity of combining simulations with analytic, sub-grid physics models, and the merits of various initial attempts to do this.

Throughout this paper, we parametrise the Hubble constant as $`H_{}=100h`$ km s-1 Mpc-1. We will need to consider two cosmologies, $`\mathrm{\Lambda }`$ CDM ($`\mathrm{\Omega }_{}=0.3`$, $`\mathrm{\Lambda }=0.7`$) and standard CDM ($`\mathrm{\Omega }_{}=1`$, $`\mathrm{\Lambda }=0`$). Matter densities are given relative to the critical density, and mass-to-light ratios ($`M/L`$) are in solar units.

## 2 Observational Constraints

### 2.1 The Global Fraction of Cold Baryons

The mass density of cooled baryons in the form of stars, $`\mathrm{\Omega }_{\mathrm{stars}}`$, can be obtained from the observed luminosity density (e.g. Blanchard, Valls-Gabaud & Mamon 1992) as long as 1) we do not miss a significant contribution from low surface brightness galaxies; and 2) the mass-to-light ratio is well understood. Fortunately, surface brightness is tightly correlated with total luminosity, and there does not appear to be a significant contribution to the baryon density from low surface brightness galaxies (Driver 1999; Cole et al. 2000). Although there is an uncomfortably large range in the Schechter parameters of the local optical luminosity function (due to the large degeneracy between $`\varphi ^{}`$ and $`L^{}`$), Fukugita, Hogan & Peebles (1998, hereafter FHP) demonstrate that the total integrated $`B`$band luminosity is constrained to within about 15%. However, the use of the optical luminosity function requires large and uncertain corrections for the $`M/L_B`$ ratio, which is sensitive to stellar populations and, thus, must be determined separately for different types of galaxy.

A more robust estimate is obtained from the $`K`$band luminosity function, since $`M/L_K`$ is less sensitive to star formation history, varying by only a factor of about two due to stellar population age differences \[Bell & de Jong 2000\]. Recently, Cole et al. \[Cole et al.  2000b\] have computed the mass function using over 17000 galaxies with redshifts from the 2dF galaxy redshift survey and magnitudes from the 2-micron all sky survey (2MASS). From the galaxy colours and stellar population synthesis models they determine $`M/L_K`$ for each galaxy (independent of $`h`$). Their ($`\mathrm{\Lambda }CDM`$) result is most sensitive to the assumed initial mass function (IMF): for a Kennicutt \[Kennicutt 1983\] IMF, $`\mathrm{\Omega }_{\mathrm{stars}}=0.0014h^1`$, while for the Salpeter \[Salpeter 1955\] IMF, $`\mathrm{\Omega }_{\mathrm{stars}}=0.0026h^1`$.

The average stellar $`M/L_K`$ from which Cole et al.  \[Cole et al.  2000b\] derive $`\mathrm{\Omega }_{\mathrm{stars}}`$ is 0.73 and 1.32 for the Kennicutt and Salpeter IMFs, respectively (independent of $`h`$). A substantially higher $`M/L_K`$ would require a large contribution from low mass stars and brown dwarfs. Although these objects are likely to be a significant component of the mass budget, there is no strong evidence that they are more abundant than expected from standard IMFs \[Fuchs et al.  1998, Gizis et al.  2000, Lucas & Roche 2000\]. In particular, Reid et al. \[Reid et al.  1999\] claim that brown dwarfs contribute no more than 15% of the mass of the Galactic disk. Earlier microlensing results which suggested that there might be a large population of dark objects in the Galaxy halo are no longer compelling, following a reanalysis of the MACHO project data (Alcock et al. 2000).

Dynamical measurements of early type galaxies \[van der Marel 1991, Bell & de Jong 2000\] favour $`M/L_K<1.1h`$; this is an upper limit to the global average because later type galaxies have lower $`M/L_K`$. We therefore prefer the estimate of $`\mathrm{\Omega }_{\mathrm{stars}}=0.0014h^1`$ based on the Kennicutt \[Kennicutt 1983\] IMF with a mean $`M/L_K=0.73`$.

Cooled gas is also present in the form of neutral and molecular gas, but these make only small contributions. The density of neutral gas is $`\mathrm{\Omega }_{\mathrm{atomic}}0.000188h^1`$, as recently calculated by the HI Parkes All Sky Survey (HIPASS) team \[Kilborn et al.  1999\]. For the amount of gas in molecular form, we adopt the relation $`\rho _{H_2}/\rho _{HI}=0.81`$ used by FHP, gleaned from the CO-based compilation of Young & Scoville \[Young & Scoville 1991\]. Together, the neutral and molecular gas only represent $`10`$% of the stellar mass, which is the fraction we will adopt, so that $`\mathrm{\Omega }_{\mathrm{cold}}=1.1\mathrm{\Omega }_{\mathrm{stars}}`$.

The total baryon content of the universe is unknown observationally, because baryons in the warm plasma phase believed to occupy normal galaxy haloes are very difficult to see \[Benson et al.  2000\]. However, if nucleosynthesis calculations are correct, the baryon density $`\mathrm{\Omega }_b`$ can be determined from deuterium abundances. The best current estimate, $`\mathrm{\Omega }_b=0.019h^2`$ \[Burles & Tytler 1998, Burles et al.  2000\], therefore implies $`f_{\mathrm{c},\mathrm{global}}=\mathrm{\Omega }_{\mathrm{cold}}/\mathrm{\Omega }_b=0.073h`$. There are indications, from recent analysis of the combined BOOMERANG and MAXIMA data, that $`\mathrm{\Omega }_b`$ is much larger than this, $`\mathrm{\Omega }_b=0.039h^2`$ (Jaffe et al. 2000), which would reduce $`f_{\mathrm{c},\mathrm{global}}`$ by a factor of two.

In summary, only about 5% of the available baryons in the universe have cooled (for $`h0.7`$ and a Kennicutt IMF). This result is slightly lower than the minimum of the range $`0.062<f_{\mathrm{c},\mathrm{global}}<0.167`$ found by FHP from optical data, but much more certain as it does not rely on knowing the relative abundances of different galaxies (disks, spheroids and irregulars) with very different (and uncertain) $`M/L_B`$ ratios. As we will show in §3, this low value is in stark contrast to the results obtained in numerical simulations which do not employ an adequate feedback model.

### 2.2 The Abundance of Hot Baryons

In the previous section we claimed that most of the baryons in the Universe are in a warm or hot phase, which is difficult to observe because it is not hot enough to be observable in X-ray radiation \[Blanchard et al.  1992, Davé et al.  2000, Benson et al.  2000\]. The existence of such a warm component is compatible with constraints on the anisotropy of the microwave background, a long as the gas temperature satisfies the constraint $`T<4h\times 10^7`$ K (at $`z<1`$) \[Wright et al.  1994\].

In clusters, however, the same gas is much hotter, and is directly observable. As originally shown by White et al. \[White et al.  1993\], the total baryon content of rich clusters, including this plasma, is fully consistent with $`\mathrm{\Omega }_b`$ from element abundance determinations, if $`\mathrm{\Omega }_{}0.3`$. FHP compile data on the baryon contributions in clusters due to stars, cold gas and plasma, and demonstrate that $`\mathrm{\Omega }_b/\mathrm{\Omega }_{}=0.112\pm 0.05`$ ($`h=0.7`$). In this case, $`\mathrm{\Omega }_b`$ is a sum of all observed baryons in clusters. Hence, for $`\mathrm{\Omega }_{}=0.3\pm 0.1`$, we have that $`\mathrm{\Omega }_b=0.034\pm 0.019`$ ($`h=0.7`$) which is consistent with the baryon fraction deduced from deuterium abundances \[Burles & Tytler 1998\]. Thus, the fact that $`\mathrm{\Omega }_{\mathrm{stars}}\mathrm{\Omega }_b`$ in the Universe is almost certainly a reflection of the fact that the warm gas has not been directly observed.

### 2.3 The Dependence on Halo Mass

Although we have demonstrated that $`f_{\mathrm{c},\mathrm{global}}`$$`0.073h`$, it is possible that the cooled fraction in some environments, rich clusters for example, could be very different. In particular, Bryan \[Bryan 2000\] claims that the efficiency of galaxy formation is dependent on halo mass, with more gas cooling in small haloes, relative to large ones. Most of the cluster data in the literature has been obtained at optical wavelengths, so the correction for $`M/L`$ variations is not as robust as for the global value computed from the $`K`$-band field luminosity function in §2.1. We will use the optimal values from FHP, who compile data from various sources, including dynamical measurements and population synthesis model estimates. The $`M/L_B`$ depends strongly on morphological type, and FHP find $`M/L_B=6.5_{2.0}^{+1.8}`$ for E/S0 galaxies; $`M/L_B=1.5\pm 0.4`$ for spiral galaxies; and $`M/L_B=1.1\pm 0.25`$ for irregular galaxies. For the morphological composition of a “typical” cluster, they determine $`M/L_B=4.5\pm 1`$, and we will use this number. Assuming a pure elliptical population would increase stellar mass estimates by $`45`$%.

For low mass clusters and groups, the X-ray data only extend out to a small fraction of the virial radius. This can lead to a strong bias, since the gas distribution is less concentrated than the stars \[Roussel et al.  2000\]. We will consider the cluster and group data recently compiled by Roussell et al. \[Roussel et al.  2000\] specifically for the purpose of addressing this issue. In that paper, care was taken to ensure that stellar and gas masses are computed to the same radius, and only those clusters for which the X-ray data extend to at least 25% of the virial radius are considered here. Nonetheless, we warn that the data for groups with $`kT<5`$ keV must be treated with caution, as their gas masses still require a substantial extrapolation to the virial radius. In addition to the standard analysis of the X-ray data, Roussel et al.  also consider hot gas mass estimates based on simulation-motivated scaling laws, which do not require an assumption of isothermal, $`\beta `$model profiles. We use the measurements made based on the scaling laws of Bryan & Norman \[Bryan & Norman 1998\]. In addition, we consider the 15 rich clusters with reliable velocity dispersions from the CNOC sample \[Carlberg et al.  1996, Carlberg et al.  1997\], with X-ray data taken from Lewis et al.  \[Lewis et al.  1999\].

The stellar masses are computed from the integrated stellar luminosity function, assuming $`M/L_B=4.5`$, and corrected for undetected galaxies by extrapolating the luminosity function with an assumed faint end slope of $`\alpha 1.2`$. Finally, we assume that the mass of cold baryons is 10% larger than the stellar mass, to account for the presence of neutral and molecular gas (see § 2.1).

We consider two approaches to compute $`f_{\mathrm{c},\mathrm{cluster}}`$ from $`M_{\mathrm{cold}}`$ in the clusters. In the first case, we compute

$$f_{\mathrm{c},\mathrm{cluster}}=M_{\mathrm{cold}}/\left(M_{\mathrm{cold}}+M_{\mathrm{hot}}\right),$$
(1)
where $`M_{\mathrm{hot}}`$ is determined from the observable X-ray emission from the hot plasma, extrapolated to the virial radius. This will be valid if the extrapolation is accurate (less likely for lower mass clusters), and will be an overestimate of $`f_{\mathrm{c},\mathrm{cluster}}`$ if there is an unobserved, warm component which does not contribute to the X-ray emission (as is the case in galactic mass haloes). We show $`f_{\mathrm{c},\mathrm{cluster}}`$ for the two samples in the bottom panel of Figure 1, for a $`\mathrm{\Lambda }CDM`$ cosmology with $`h=0.7`$. Both $`M_{\mathrm{cold}}`$ and $`M_{\mathrm{hot}}`$ are distance dependent measurements, with $`M_{\mathrm{cold}}h^1`$ and $`M_{\mathrm{hot}}h^{1.5}`$. Since, usually, $`M_{\mathrm{hot}}M_{\mathrm{cold}}`$, $`f_{\mathrm{cool}}h^{0.5}`$ when calculated from Equation 1.

Most of the $`f_{\mathrm{c},\mathrm{cluster}}`$ values in the bottom panel of Figure 1 are higher than the global fraction $`f_{\mathrm{c},\mathrm{global}}`$ computed in §2.1, shown as the horizontal line. There is only weak evidence for a trend, in the sense that $`f_{\mathrm{c},\mathrm{global}}`$ is lower in higher temperature systems. However, this trend is much weaker than suggested by Bryan \[Bryan 2000\] and, we caution that the $`kT<5`$ keV groups may still be biased toward high fractions because of the limited extent of the X-ray data; none of these groups have X-ray emission detected beyond half the virial radius. Also, both the random and systematic uncertainties in observationally determining $`M_{\mathrm{hot}}`$, especially for the low temperature groups, are large. This is emphasized in the top panel of Figure 1, in which the same data are plotted, but $`f_{\mathrm{c},\mathrm{cluster}}`$ is computed as

$$f_{\mathrm{c},\mathrm{cluster}}=\frac{M_{\mathrm{cold}}}{M_{\mathrm{total}}}\frac{\mathrm{\Omega }_{}}{\mathrm{\Omega }_b},$$
(2)
where $`M_{\mathrm{total}}`$ is the total dynamical mass of the cluster, and we take $`\mathrm{\Omega }_b`$ from Burles & Tytler \[Burles & Tytler 1998\]. This does not depend on the difficult measurement of $`M_{\mathrm{hot}}`$, but does require an estimate of $`\mathrm{\Omega }_{}`$, which is 0.3 in our $`\mathrm{\Lambda }CDM`$ cosmology. Note that, in this case, $`f_{\mathrm{c},\mathrm{cluster}}`$$`\mathrm{\Omega }_{}h`$. Not only has the apparent trend in $`f_{\mathrm{c},\mathrm{cluster}}`$ with temperature now virtually disappeared, but the scatter has been reduced considerably. In this case the $`f_{\mathrm{c},\mathrm{cluster}}`$ values are in much closer agreement with $`f_{\mathrm{c},\mathrm{global}}`$; the remaining difference is not very compelling, since a slightly lower $`\mathrm{\Omega }_{}0.2`$ would result in concordance between the cluster and global fractions.

These results differ from those of Bryan \[Bryan 2000\], who claims that $`f_{\mathrm{c},\mathrm{cluster}}`$ increases from cluster to group scales. However, this trend could well be the result of a bias due to the limited radial extent of the X-ray observations, $`R_x`$, in the group samples \[Roussel et al.  2000\]. For example, from the Mulchaey et al.  \[Mulchaey et al.  1996\] sample considered by Bryan \[Bryan 2000\], the groups with the smallest $`R_x`$ have among the highest $`f_{\mathrm{c},\mathrm{cluster}}`$ values. On the other hand, the two groups for which $`R_x>0.5`$ Mpc have the lowest $`f_{\mathrm{c},\mathrm{cluster}}`$. This emphasizes the importance of ensuring the X-ray data extend out to at least a sizable fraction of the virial radius, as we have done when considering the Roussel et al.  compilation. We have shown that there is no discernable trend of $`f_{\mathrm{c},\mathrm{cluster}}`$ with temperature in this sample, particularly if the baryon fraction is assumed from nucleosynthesis arguments, rather than relying on the measurement of $`M_{\mathrm{hot}}`$.

## 3 Numerical Simulations

Hydrodynamical numerical simulations which attempt to follow the effects of gas cooling have proved very difficult to perform across the full range of the gravitational mass hierarchy. Allowing gas to dissipate its energy via radiative cooling generates enormous density contrasts and a huge range of interesting spatial scales within the model, even before the intrinsically sub-resolution processes of star formation and the feedback of energy due to supernovae and stellar winds are considered.

The earliest attempt at such a simulation was made by Thomas & Couchman \[Thomas & Couchman 1992\], though the model which included radiative cooling was only briefly discussed. This model had a very high mass resolution threshold, allowing only the largest galaxies to form, and only a small fraction of the gas to cool. This was quickly followed by the work of Katz & White \[Katz & White 1993\] who produced a model that was several years ahead of its time. They demonstrated that the formation of a Virgo-like cluster, including the effects of cooling, could be followed by a computer simulation. However, the model had some problems: in particular, as the authors make clear, much more gas cools than is observed. These results are shown as the large, filled circles in Figure 2. The lower point shows the actual fraction of cold baryons in the simulation; the higher point to which it is connected shows this fraction corrected by extrapolating the simulated luminosity function to include the additional cooled gas within unresolved objects, as described in Katz & White. The simulation of Katz & White was stopped at $`z=0.13`$, so the redshift zero fraction will be even higher.

Recently there has been renewed interest in this problem. Suginohara & Ostriker \[Suginohara & Ostriker 1998\] also produced a cluster model, and noted that additional physical processes (such as star formation and feedback) would be required to prevent the cooling of an excessively large fraction of the gas; unfortunately, the fraction of cooled baryons was not published in their paper so we cannot consider their results in Figure 2. Their work was followed by that of Pearce et al. \[Pearce et al.  1999, Pearce et al.  2000\], Lewis et al. \[Lewis et al.  2000\] and Davé et al. (2000). Lewis et al.  used a more advanced code to improve upon the earlier simulation of Katz & White \[Katz & White 1993\], and they also found that around 40% of the gas within the cluster cools. This is shown in Figure 2 as filled triangles; the upper triangle includes a correction for galaxies below the resolution limit, by extrapolating the resultant luminosity function assuming a faint end slope $`\alpha =0.96`$. Davé et al.  used both a parallel treecode and a Eulerian grid based code in several cosmologies and, although the cluster-by-cluster fractions are not quoted, find a universal cooling fraction of around 30% to 40% for their high resolution simulations, in both cases. This is shown as a lower limit in Figure 2, since it does not include a correction for galaxies formed in haloes below the resolution limit.

Simulations like those discussed above which cool too much gas cannot be expected to model galaxy formation in a physically correct way. An alternative, heuristic approach can be taken to force the simulations to satisfy the constraint on $`f_{\mathrm{c},\mathrm{global}}`$, by limiting the resolution (see § 3.1). This is the method used by Pearce et al. \[Pearce et al.  2000\], who deliberately set their mass resolution to cool $``$15% of the gas within their simulation volume; the cooled fractions for the twenty most massive of these clusters are shown in Figure 2 as the filled squares (both before and after correction for sub-resolution galaxies). By construction, the cooled gas fractions of these simulated haloes are closer to the observed values, though they are still probably too high. The obvious problem with this simulation is that it does not allow the formation of small galaxies.

### 3.1 Resolution

It has long been theorised that progessively increasing the resolution within a simulation leads to a continual increase in the fraction of material which cools \[Cole 1991, White & Frenk 1991, Suginohara & Ostriker 1998\], a situation commonly referred to as the cooling catastrophe. We demonstrate this in Figure 3, by presenting the results of standard cold dark matter simulations, at different resolutions (the mass of the smallest object that can effectively cool). This figure shows how the fraction of gas particles in the cold phase increases as resolution is increased. At the low resolution end, the relation is too steep, due to the artificial heating problems described by Steinmetz & White \[Steinmetz & White 1997\]. However, the trends discussed below are qualitatively retained, even when this problem is overcome \[Kay 2000\]. Apart from this effect, the shape of the curve will be approximately independent of the simulations. However, the scale along the abscissa will vary; that is, the curves may be translated laterally, for example, by changing $`\mathrm{\Omega }_b`$. More subtly, in a model with a positive cosmological constant, structure forms earlier and so more gas cools at the same threshold mass. This causes the curves plotted on Figure 3 to move to the right by about a factor of two. Changes in metallicity can affect the result in the same way, though even primoridal abundances cannot prevent the cooling catastrophe \[Davé et al.  2000\]. In addition, differences in the numerical implementation can also shift the curve by factors of order two.

Figure 3 makes it clear why some groups find cooled gas fractions around 40% whilst others find much lower fractions. Values around 40% will be achieved for a wide range of simulation resolutions, as this is the value obtained when the resolution mass is less than the characteristic mass of the halo mass function. More gas will cool as the resolution is increased but a large increase is required to see a significant change, since the slope of the mass function is shallow in this regime. Any required cold gas fraction below 40% can be obtained by choosing the resolution appropriately, although care must be taken to account for artificial heating \[Steinmetz & White 1997\]. However, in order to get a good match to the observed gas fraction within a standard cold dark matter simulation the resolution limit must be set to $`10^{11}\mathrm{M}_{}`$ or higher, which is clearly unphysically high. This is the situation we referred to as the “cooling crisis”; it is a crisis because its resolution will result in a significant improvement in the utility of simulation results.

## 4 Discussion

We have demonstrated that it is easy to exceed the observed cooled gas fraction with numerical simulations of sufficient resolution. The failure to match this constraint reflects the well-known need for some method of reheating the gas in order to prevent excessive star formation at high redshift \[Larson 1974, White & Rees 1978, White & Frenk 1991\]. Until an effective feedback model is implemented in simulations, it will not be possible to obtain physically illuminating results.

Present attempts to incorporate feedback within cosmological simulations all suffer from the same drawback, that the physical processes responsible for the redistribution of energy occur on scales well below those that can presently be resolved. Hence, an accurate treatment of the propagation of energy from the local ISM to super-galactic scales is currently an intractable problem. However, several groups have attempted to include feedback by injecting the energy into regions at the resolution limit of the simulation. For example, an approach that has been taken in smoothed particle hydrodynamics (SPH) simulations is simply to add available energy from newly-formed ”star” particles to the thermal energy of the surrounding gas \[Katz et al.  1996, Lewis et al.  2000, Davé et al.  2000\]. However, the majority of this gas is cold and dense, and so any excess energy is rapidly re-radiated. It is plausible that this problem is due to the failure of the SPH simulations to resolve a multiphase medium, in which the feedback would be able to create hot bubbles of diffuse gas, with significantly longer cooling times than the cold dense gas from which they were spawned. Several groups have attempted to ”correct” for this by preventing reheated gas from radiatively cooling over a finite timescale (e.g. Gerritsen 1997; Springel 2000; Thacker & Couchman 2000). Another method is to supply energy in kinetic form (e.g. Navarro & White 1993; Kay 2000), which is efficient at limiting the fraction of cooled gas. However, the effectiveness of this implementation may be artificial since the SPH simulations fail to resolve shocks that would efficiently thermalize the reheated material \[Katz et al.  1996\].

An alternative approach to feedback is that of Cen & Ostriker \[Cen & Ostriker 1999\], who use a grid-based cosmological code, and implement analytic rules, similar to those of semi-analytic models \[Somerville & Primack 1999, Cole et al.  2000a, Kauffmann et al.  1999, Somerville et al.  2001\] to regulate the formation of galaxies within a single cell. Unfortunately, the resolution of grid-based methods such as this is still too coarse ($`100h^1`$ kpc in the highest resolution simulation of Cen & Ostriker) to provide useful results on galaxy scales. It is not clear that the analytic rules implemented to model the sub-grid physics in this code are robust to changes in resolution, as changes in the cell size can result in large changes in the gas densities, for which the analytic formalism may not appropriately compensate. A more promising route may be the use of adaptive-mesh refinement codes \[Bryan & Norman 1998, Bryan & Voit 2001\], which will allow the implementation of an analytic treatment of gas cooling within cells of sufficient resolution.

The goal now is to find a physically motivated feedback model which is able to sufficiently reduce $`f_\mathrm{c}`$. It is therefore useful to consider, qualitatively, what is required of such a model. The amount of gas which can cool in a given halo cannot be computed by considering a halo structure at a single redshift and treating it as if it had existed unchanged over a Hubble time. For example, the cooling time in clusters today is much longer than the Hubble time, and this has been used to explain why not all of the gas in these environments has cooled \[Binney 1977, Silk 1977, Rees & Ostriker 1977\]. However, this argument does not hold in hierarchical models, where cluster progenitors at earlier epochs were low mass haloes in which the cooling time was short, and simulations of hierarchical cluster formation demonstrate that a large fraction of the baryons cool by the present day. To treat the problem analytically, it is necessary to consider the hierarchical growth of a halo through the plane of Figure 4. The smooth curved lines represent the evolution of structure, for a standard CDM model ($`\sigma _8=0.8`$, $`\mathrm{\Gamma }=\mathrm{\Omega }_{}h=0.5`$); the central line shows the median mass halo which virializes as a function of redshift, and the two others bracket the region in which 50% of the mass in the Universe collapses into haloes. In a manner similar to that of Rees & Ostriker \[Rees & Ostriker 1977\], we can define an efficient cooling mass, $`M_{\mathrm{cool}}`$, as the mass which can cool all the gas out to its virial radius, within a Hubble time at redshift $`z`$. The cooling rate, which depends on halo density and temperature, is computed using the models of Raymond et al. \[Raymond et al.  1976\], with one-third solar metallicity. This cooling mass, for the standard CDM model described above, is shown in Figure 4, and is given by $`M_{\mathrm{cool}}10^{12}`$$`\mathrm{M}_{}`$, almost completely independent of redshift for $`z<10`$. This non-trivial result is a consequence of the fact that, when cooling is dominated by line emission, the increase in the cooling rate with increasing redshift is closely balanced by the decrease in the Hubble time. This forms an approximate upper mass limit to efficient cooling; haloes with masses smaller than this limit can cool all of their gas, out to the virial radius, while to the right of this line, only a rapidly decreasing fraction of the gas is able to cool.

At low redshifts, most of the haloes which are collapsing are more massive than $`M_{\mathrm{cool}}`$; therefore, cooling plays a minor role today. By $`z1`$, however, the characteristic mass crosses the $`M_{\mathrm{cool}}`$ threshold, and cooling becomes a dominant process, leading to the cooling catastrophe. In simulations, the resolution threshold imposes a mass limit, $`M_{\mathrm{th}}`$, which prevents smaller mass haloes from cooling. Efficient cooling is therefore restricted to the shaded region in Figure 4. As $`M_{\mathrm{th}}`$ is lowered, more and more mass cools, and will quickly exceed the tight observational constraints. If the threshold is high enough to prevent efficient cooling at high redshift, no low mass galaxies are formed.

The final fraction of cooled gas then depends on the trajectory a particular halo takes through this plane, as it grows. This growth will occur via both smooth accretion, and discrete jumps due to mergers. The resolution threshold imprints a scale on the accumulation of cooled baryons, as the final fraction $`f_\mathrm{c}`$ in a given halo will depend on its trajectory only after it becomes sufficiently resolved. This could have an effect on the dependence of $`f_\mathrm{c}`$ on halo mass in simulations (but see Kay 2000).

It is clear that any analytic feedback model which behaves in a manner similar to a sharp, fixed mass threshold will fail to simultaneously match the observed $`f_{\mathrm{c},\mathrm{global}}`$ and still produce a reasonable galaxy luminosity function. A physical feedback scheme, however, does not provide a hard mass threshold, but gradually alters the cooling efficiency of haloes as a function of mass and redshift; thus the shaded region in Figure 4 becomes a “fuzzy” region without a well defined lower boundary. There are many models for such feedback, including supernovae or AGN reheating \[Larson 1974, White & Rees 1978, Lacey & Silk 1991, Cole 1991, Wu et al.  1998, Bower et al.  2000, Mac Low & Ferrara 1999, Efstathiou 2000\], ionizing radiation from the ultraviolet background \[Efstathiou 1992, Weinberg et al.  1997\], and an entropy threshold established at high redshift \[Kaiser 1991, Blanchard et al.  1992, Ponman et al.  1999, Balogh et al.  1999\]. A very successful implementation of a feedback scaling law is that used by semi-analytic models of galaxy formation \[Cole et al.  2000a, Somerville & Primack 1999, Kauffmann et al.  1999, Somerville et al.  2001\], in which the mass of reheated gas is assumed to scale in proportion with the halo circular velocity raised to some power $`\alpha `$. Though this scaling is constructed to allow the models to match the observed luminosity function, it also succeeds in limiting $`f_\mathrm{c}`$ to more reasonable values (though generally still high compared with the constraints described here). Furthermore, there may be some theoretical motivation for this scaling relation from supernova-driven wind models \[Efstathiou 2000\]. The next step in this work is to combine these models with realistic halo merger trees to determine, analytically, how they affect the cooling rate as a function of halo mass and redshift, for implementation into numerical simulations.

## 5 Conclusions

In this paper, we have presented an up-to-date review of the cooling catastrophe. Although the issues we have discussed have been known for a long time, recent advances in observations and simulation techniques warrant a reevaluation of the situation, which can be summarized in two key points:

* The global fraction of baryons in the form of stars and cold gas is now well constrained observationally from the $`K`$band luminosity function to be $`f_{\mathrm{c},\mathrm{global}}`$$`0.073h`$. Although there is some uncertainty about how this may vary with halo mass, there is no convincing evidence for a significant difference from the global value.
* The amount of gas which cools in simulations is directly related to the resolution, unless a working reheating scheme is implemented. No current simulation is able to produce both (1) the correct $`f_{\mathrm{c},\mathrm{global}}`$ and (2) a population of low mass galaxies, without implementing analytic models of the sub-grid physics on uncomfortably large scales.

The use of simulations to aid our understanding of galaxy properties and of cluster scaling laws therefore awaits the implementation of a sub-grid, physical feedback scheme which allows the models to satisfy observational constraints on $`f_{\mathrm{c},\mathrm{global}}`$ in a resolution-independent way.

## Acknowledgements

MLB acknowledges support from a PPARC rolling grant for extragalactic astronomy and cosmology at Durham. We thank Simon White, Eric Bell, Greg Bryan and Cedric Lacey for several suggestions and corrections which greatly improved this paper. Finally, we thank the referee for providing a thorough report which allowed us to improve the focus and clarity of this work.


format: text

text:
# 1 Candidate Halo White Dwarfs. WD names are new discoveries; otherwise, Luyten LHS or LP names are given; J0014–3937 was discovered by Scholz et al. (20), who classify this object as “cool no H𝛼”, and F351–50 was discovered by Ibata et al. (3). The tangential velocity, 𝑣ₜₐₙ(est) in km s-1, is an estimate based on photometric parallax. Distances are estimates based on a color-magnitude relation (see text). Star positions are given in right ascension (RA) and declination (Dec.) at the equinox and epoch of J2000. The proper motion, 𝜇, and position angle, PA, are given in ′′ year-1 and degrees, respectively. The photometric measurements in 𝐵_J, 𝑅59𝐹, and 𝐼_N are given in magnitudes, with uncertainties of ±0.1 mag. on average. The distance is in parsecs from the Sun.

Direct Detection of Galactic Halo Dark Matter

B.R. Oppenheimer, N.C. Hambly, A.P. Digby, S.T. Hodgkin, D. Saumon111B.R. Oppenheimer (bro@astron.berkeley.edu) is at Astronomy Department, University of California–Berkeley, Berkeley, CA 94720-3411, USA. N.C. Hambly (nch@roe.ac.uk) and A. P. Digby (apd@roe.ac.uk) are at Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, UK. S.T. Hodgkin (sth@ast.cam.ac.uk) is at Institute of Astronomy, Cambridge University, Madingley Road, Cambridge, CB3 0HA, UK. D. Saumon (dsaumon@cactus.phy.vanderbilt.edu) is at Department of Physics and Astronomy, Vanderbilt University, Nashville, TN 37235 USA.

The Milky Way Galaxy contains a large, spherical component which is believed to harbor a substantial amount of unseen matter. Recent observations indirectly suggest that as much as half of this “dark matter” may be in the form of old, very cool white dwarfs, the remnants of an ancient population of stars as old as the Galaxy itself. We conducted a survey to find faint, cool white dwarfs with large space velocities, indicative of their membership in the Galaxy’s spherical halo component. The survey reveals a substantial, directly observed population of old white dwarfs, too faint to be seen in previous surveys. This newly discovered population accounts for at least 2% of the halo dark matter. It provides a natural explanation for the indirect observations, and represents a direct detection of Galactic halo dark matter.

Dark matter in the spherical halo of the Milky Way Galaxy has been inferred because the gravitational field due to the known distribution of luminous matter, primarily stars, cannot explain the observed rotational characteristics of the Galaxy’s spiral disk. A substantial portion of this unseen matter may be old, very cool white dwarfs (1-4). A white dwarf is the extremely dense end-state in the evolution of stars with masses less than about eight times the mass of the Sun ($`M_{}`$). Once a star becomes a white dwarf, it no longer produces energy through nuclear fusion, and therefore cools and fades. The first four examples of ultracool white dwarfs, whose temperatures are below 4000 K, were identified only in the past 2 years (2, 3, 5-7), principally because they were too faint to detect in previous surveys. They have spectral energy distributions that differ dramatically from those of the hotter white dwarfs, consistent with white dwarf atmosphere models (8-10). The difference is due to the formation of H2 molecules in white dwarf atmospheres with effective temperatures below $`4500`$ K. Although H2 molecules are symmetric and thus have no dipole moment, in the high densities of white dwarf atmospheres, collisions between the molecules are common. These collisions induce momentary dipole moments, which produce opacity at wavelengths longer than 0.6 $`\mu `$m (11). We conducted a survey to search for nearby, high-proper-motion white dwarfs that could be halo members and might exhibit this opacity. Regardless of their spectral appearance, nearby halo stars can be distinguished from disk stars by their space motion, because the two populations possess different kinematic properties.

Observations

Our survey used digitized, photographic plates in the $`R59F`$ and $`B_\mathrm{J}`$ passbands from the SuperCOSMOS Sky Survey (12-14). $`R59F`$ and $`B_\mathrm{J}`$ roughly correspond to optical wavelengths of 0.59 and 0.45 $`\mu `$m respectively. We searched for objects with proper motions, $`\mu `$, between 0.33 and 10.00 ′′ year-1 as faint as $`R59F=19.8`$, using plates near the South Galactic Cap (SGC) at three epochs in each field. The SGC plates comprise 196 fields (15). The non overlapping area of each field is 25 square degrees. Hence, the area of our survey is 4900 square degrees, or about 12% of the sky. Because of image blending problems and large halos of scattered light around bright stars, the actual area surveyed is 4165 square degrees, or 10% of the sky. This is several times larger than the areas searched in two previous surveys (3, 16). At the very faint end of our sample, we found many objects not included in the Luyten catalogs of proper motion stars (17), primarily because Luyten’s survey in the Southern Hemisphere did not attempt to find objects as faint as those we have found. These previously undetected stars are the principal members of our sample of halo white dwarfs.

We used the technique of reduced proper motion (RPM) (18, 19) to identify white dwarfs by their subluminosity in comparison to main-sequence stars and their high intrinsic space motions. RPM is defined as $`H_\mathrm{R}=R59F+`$ 5 log $`\mu +5`$, and is an estimate of the absolute magnitude based on the proper motion, $`\mu `$, in arcseconds per year. We selected all 126 objects in the RPM diagram (Fig. 1) having either $`H_\mathrm{R}>22.5`$ or lying in the subluminous regions as potential halo white dwarfs. Sixty-three of these were previously cataloged (12, 20). Of these 63 known, high-proper-motion stars, 29 had no published spectroscopic follow-up. This yielded an observing list of 92 stars for follow-up, including 63 new discoveries.

Over four nights (11 to 14 October 2000) at the Cerro Tololo Interamerican Observatory’s V. M. Blanco 4-m Telescope, we obtained spectra in the wavelength range 0.48 to 0.98 $`\mu `$m for 69 of the 92 candidates (filled symbols in Fig. 1). We inspected the spectra within minutes of collecting the data, so we were able to discern by the middle of the second night which types of objects populated which parts of the RPM diagram (Fig. 1). For example, as expected, we found only white dwarfs with H$`\alpha `$ features on the far left of the subluminous sequence. Dwarf and subdwarf M–type stars are restricted to the right, which is also the bottom of the main sequence. We concentrated our observing resources on the cool objects, which occupy the lower part of the subluminous sequence. Of the 69 candidate objects observed, 16 are M-dwarfs or M-subdwarfs, 2 are hot He white dwarfs, 13 are white dwarfs with H$`\alpha `$ features, and 38 are new cool white dwarfs (Fig. 2a), a few of which are probably cooler than WD0346$`+`$246, the prototypical ultracool white dwarf, which exhibits the H2 opacity (2, 7, 21, 22).

Two of the 38 new cool white dwarfs have unusual spectra (Fig. 2b). One, LHS 1402, has a spectrum very similar to those of the peculiar stars LHS 3250 (5, 7) and SDSS 1337+00 (6), but with a steeper slope toward longer wavelengths, suggesting a cooler temperature. The other object, WD2356$``$209, possesses a bizarre spectrum, incomparable to any other known object. We reanalyzed the data on this object several times and found no evidence for residual instrumental effects.

White Dwarfs in the Galactic Halo

The 38 cool white dwarfs in our sample show no spectral lines (e.g., Fig. 2a), preventing determination of their radial velocities. Nevertheless, we can estimate their space motions from the tangential velocity, $`v_{\mathrm{tan}}=\mu d`$, where $`\mu `$ is the proper motion and $`d`$ is the distance between each star and Earth. To estimate distances, we used a photometric parallax relation derived from a linear least–squares fit to the cool white dwarf sample of Bergeron, Ruiz, and Leggett (23), supplemented by the ultracool WD0346+246, yielding the absolute magnitude in the $`B_\mathrm{J}`$ filter, $`M_{\mathrm{B}_\mathrm{J}}=12.73+2.58(B_\mathrm{J}R59F)`$. Our calculations of the spectra of very cool white dwarfs also show a linear relation between $`M_{\mathrm{B}_\mathrm{J}}`$ and $`B_\mathrm{J}R59F`$. The color turnoff in $`B_\mathrm{J}R59F`$ due to collision-induced absorption by H2 molecules in the pure hydrogen atmospheres becomes apparent in $`B_\mathrm{J}R59F`$ only at effective temperatures below 2500 K. The distances derived by this method are listed in Table 1. The scatter of the data around the least-squares fit results in a 20% uncertainty in our distance estimates. For LHS 147 and LHS 542, the only stars in Table 1 with measured distances (16, 24), our estimates and the measured distances agree within 1 and 20%, respectively.

In order to use the velocities of the stars as an indicator of membership in the Galactic Halo, we transformed $`v_{\mathrm{tan}}`$ into the components $`U`$, $`V`$, and $`W`$, in galactic coordinates. $`U`$ is radial away from the Galactic center, $`V`$ is in the direction of rotation, whereas $`W`$ is perpendicular to the Galactic disk. The calculations used to derive these velocities take into account the deviation of the velocity of the Sun with respect to the average velocity of nearby stars in the Galactic disk. Because our survey was centered around the SGC, our observations are most sensitive to the $`U`$ and $`V`$ velocities (Fig. 3). Therefore, the most reasonable assumption is that $`W`$ is zero for all of our stars. Zero is also the average value of $`W`$ for all stars in the Milky Way. The point where $`U`$, $`V`$, and $`W`$ are all equal to zero defines the so-called “local standard of rest,” a frame of reference which is rotating in the disk about the Galactic center at the average rotation velocity of the disk near the Sun. Halo objects, because they have an isotropic velocity distribution with little or no overall rotation, should lag behind the local standard of rest with a distribution centered near $`V=220`$ km s-1, because the local standard of rest orbits the Galaxy at 220 km s-1 (25).

Our survey contains 38 stars that unambiguously possess halo kinematics (Fig. 3). We selected a sample of halo stars by accepting all objects whose velocities exceed the 2$`\sigma =94`$ km s-1 velocity dispersion for the old disk stars (25). Our velocity cut excludes some halo stars with lower velocities, but permits only minimal contamination of the sample by disk stars. By choosing stars above 2$`\sigma =94`$ km s-1, we excluded 95% of the disk stars.

Furthermore, if a portion of our sample contains binary stars \[the fraction of white dwarfs that are in white dwarf-white dwarf pairs is between 2 and 20% (26)\], we expect that the velocity distribution of halo stars (Fig. 3) will be biased toward the point $`(V,U)=(0,0)`$. This is because we assumed that all of the stars are solitary in our estimate of distance. For those that are actually binary, we underestimated the intrinsic luminosity by, at most, a factor of 2 and the distance by a factor of $`\sqrt{2}`$. Therefore the values of $`U`$ and $`V`$ are underestimated by the same factor of $`\sqrt{2}`$. This bias is apparent in our sample (Fig. 3), implying that while we may have included a few disk stars in our sample, we may also have excluded a few binary halo stars.

Because we searched for stars having $`\mu >0.33^{\prime \prime }`$ year-1 there is a lack of stars with very low space motions. Furthermore, because of large differences in the epochs of the plates, an object with a proper motion above 3′′ year-1 will have a displacement of at least 30′′ between the plates. This compromises the detection of objects with extremely high space motions: even though we searched for objects having annual motions of up to 10′′, our survey probably had less than a 10% chance of finding very faint stars moving by 3′′ year-1 or more. These biases against stars with very low and very high space motions result in the dearth of points within 50 km s-1 of the point $`(U,V)=(0,0)`$ and the paucity of points at the extreme left of the $`UV`$ velocity diagram (Fig. 3). Therefore, we underestimated the numbers of halo white dwarfs, and we measured a lower limit to the space density of these objects.

The subset of objects with velocities outside the 2$`\sigma `$ velocity dispersion of the old disk stars (Fig. 3) contains 24 cool white dwarfs and 14 white dwarfs with H$`\alpha `$ features (Table 1). An initial estimate of the space density of these objects is obtained using only the 12 coolest white dwarfs, whose average absolute magnitude, $`M_{\mathrm{R59F}}`$, is 16. With the survey depth of 19.8 and the survey area, a rough volume is computed yielding a space density of $`2\times 10^4`$ pc-3, or $`1.2\times 10^4`$ $`M_{}`$ pc-3.

To obtain a more precise estimate of the space density of all 38 of the halo white dwarfs in our survey, we used the $`1/V_{\mathrm{max}}`$ technique (27), which determines the maximum volume, $`V_{\mathrm{max}}`$, in which the survey could have found each star, given the brightness detection limit. The sum of $`1/V_{\mathrm{max}}`$ for all of the stars in the sample is the number density. Because our survey was limited to stars brighter than $`R59F_{\mathrm{lim}}=19.8`$, we have computed $`V_{\mathrm{max}}`$ using the absolute $`R59F`$ magnitudes, $`M_{\mathrm{R59F}}`$, implied by the distances inferred above.

$`V_{\mathrm{max}}=(4/3)\pi d_{\mathrm{max}}^3\mathrm{\Omega }/(4\pi )`$ pc3 (1)

Here, $`\mathrm{\Omega }`$ is the surveyed area in steradians, and $`d_{\mathrm{max}}`$, the maximum distance in parsecs (pc) that determines $`V_{\mathrm{max}}`$, is the minimum of the two following relations.

log $`d_{\mathrm{max}}(\mathrm{pc})=0.2(R59F_{\mathrm{lim}}M_{\mathrm{R59F}})+1`$ (2)

$`d_{\mathrm{max}}`$ (pc) $`=d\mu /0.33`$ (3)

The number density derived in this manner is $`1.8\times 10^4`$ pc-3. Using 0.6 $`M_{}`$ as the average white dwarf mass in the halo (28), the local space density of halo white dwarfs is $`1.1\times 10^4`$ $`M_{}`$ pc-3. However, our survey is not complete to the $`R59F=19.8`$ level. A measure of the uniformity of our sample is given by the average value of $`V/V_{\mathrm{max}}`$, where $`V`$ is the volume of space given by the distance to the given star and the survey area. In the case of this survey, the average of $`V/V_{\mathrm{max}}`$ is 0.46. For a sample distributed uniformly in space, this value should be 0.5. An average of $`V/V_{\mathrm{max}}=0.5`$ is found for $`R59F=19.7`$, indicating that all stars brighter than 19.7, with the proper motion constraint mentioned above, were found in the survey. Recalculating the space density using $`R59F_{\mathrm{lim}}=19.7`$ yields $`1.3\times 10^4`$ $`M_{}`$ pc-3.

The density of white dwarfs in the halo predicted from the subdwarf star counts and a standard initial mass function is 1.3$`\times 10^5`$ $`M_{}`$ pc-3 (29), which is ten times smaller than the value we calculate. This means that star formation in the early Galaxy must have favored higher-mass stars that would evolve into the white dwarfs we are detecting now. At the same time, this early star formation must have under produced the low-mass subdwarfs relative to the more recent star formation processes observed in the disk (4).

The estimated density of halo dark matter near the Sun is approximately 8$`\times 10^3`$ $`M_{}`$ pc-3 (30). Thus, the population we have detected accounts for about 2% of the local dark matter. We treat this number as a lower limit for the reasons discussed above and because it appears that we have not probed the entire range of intrinsic luminosities of these objects. Indeed, the number density as a function of $`M_{\mathrm{R59F}}`$ suggests that we are only detecting the rising power law of the luminosity function and have not seen any indication of the expected turnover due to the finite life time of the halo white dwarf population (31). This means that there may be an undetected, larger population of even fainter and cooler white dwarfs in the halo.

Discussion

Microlensing experiments have indirectly revealed a population of massive, compact objects in the line of sight to the Large Magellanic Cloud. The Massive Compact Halo Object project (MACHO) (1) estimates that these compact objects contribute between 8 and 50% of the local halo dark matter at the 95% confidence level, and the inferred lens mass suggests that they are white dwarfs. The similar Expérience pour la Recherche d’Objets Sombres project (EROS) places a strong upper limit of 35% for such objects (32). Because the halo white dwarf space density that we derive from our survey represents a lower limit and is consistent with the two measurements mentioned above, the population of white dwarfs with halo kinematics that we have discovered provides a natural explanation for the microlensing results.

A number of the white dwarfs revealed in this survey may be cooler than any other white dwarf previously known. There are at least three objects (Table 1) that are comparable in temperature to WD0346$`+`$246 (7), representing a doubling of the number of ultracool white dwarfs known. The nature of cool white dwarfs can be discussed on the basis of their positions in color-color diagrams. Although a combination of optical and infrared colors is necessary at these low effective temperatures (7, 23), for now, we must contend with the limited spectral coverage of photographic plates for this sample. The comparison of $`B_\mathrm{J}`$, $`R59F`$, and $`I_\mathrm{N}`$ (a bandpass centered on 0.8 $`\mu `$m) colors is a poor diagnostic of atmospheric composition for effective temperatures, $`T_{\mathrm{eff}}`$, greater than 3500 K (Fig. 4). Given the large photometric error bars, there is little that distinguishes the colors of these stars from a sample of cool white dwarfs with disk kinematics (33).

Stars with $`B_\mathrm{J}R59F<1.2`$ are relatively warm, with $`T_{\mathrm{eff}}>5000`$K (except for the strange stars LHS 3250, SDSS 1337$`+`$00 and LHS 1402). They are also the most luminous and most distant stars in the sample (Table 1). Although it is not possible to guess the atmospheric composition of any of these stars from the present data, we point out that the sample may contain a few unusual objects, in addition to the peculiar LHS 1402 and WD2356$``$209 (Fig. 2b). Other inferences, however, can be made from the photographic color-color diagram (Fig. 4), using stars previously studied. LHS 542 has a pure He atmosphere and $`T_{\mathrm{eff}}=4747`$K (24). WD0346$`+`$246 and F351$``$50 are similar to each other with predominantly He atmospheres, a very small admixture of H, and $`T_{\mathrm{eff}}`$ of $`3750`$K and $`3500`$K, respectively (7). The most extreme star of the sample, and also the faintest, is WD0351$``$564 with $`B_\mathrm{J}R59F=1.98`$. The stars with $`B_\mathrm{J}R59F>1.6`$ (WD0351$``$564, WD0205$``$053 and WD0345$``$362) may have either pure He atmospheres with $`T_{\mathrm{eff}}>4000`$K or, more interestingly, be other examples of very cool white dwarfs with He atmospheres polluted by a small amount of H, like WD0346$`+`$246. The stars most likely to have pure H atmospheres residing below the color turnoff due to strong H2 collision-induced opacity are WD0135$``$546, WD0340$``$330, and, perhaps, WD0351$``$564.

We still do not understand the nature of the three strange objects LHS 3250, LHS 1402, and SDSS 1337$`+`$00 (5-7), all of which lie at $`R59FI_\mathrm{N}`$ colors near $``$0.5 (well outside of Fig. 4). The kinematics of these stars, with the exception of LHS 1402, suggest that they may be disk objects. Although they are certainly cool white dwarfs, we know neither their effective temperatures nor their atmospheric compositions. They may represent another, perhaps unusual, stage in the spectral evolution of cool white dwarfs.

The direct detection of a significant white dwarf halo population opens promising avenues of investigation. Besides the opportunity to study individual white dwarfs that are older and cooler than any known so far, dedicated studies of this population of stellar remnants will bear directly on the ancient history of the Galaxy. For instance, a determination of the luminosity function of halo white dwarfs and a characterization of the low luminosity cutoff will reveal the age of the first generation of halo stars, lead to a determination of its initial mass function, and provide strong constraints on formation scenarios of the Galaxy. The interpretation of observations of young galaxies at high redshifts will benefit from a better understanding of the early history of our Galaxy. Finally, the complete characterization of this new stellar population will reveal whether it is responsible for the microlensing.

References

1. C. Alcock et al., Astrophys. J. 542, 281 (2000).

2. S. T. Hodgkin et al., Nature 403, 57 (2000).

3. R. Ibata, M. J. Irwin, O. Bienaymé, R. Scholz, J. Guibert, Astrophys. J. Lett. 532, L41 (2000).

4. G. Chabrier, Astrophys. J. Lett. 513, L103 (1999).

5. H. C. Harris et al., Astrophys. J. 524, 1000 (1999).

6. H. C. Harris et al., Astrophys. J. Lett., 549, L109 (2001).

7. B. R. Oppenheimer et al., Astrophys. J., 550, 448 (2001).

8. B. Hansen, Nature 394, 860 (1998).

9. P. Bergeron, D. Saumon, F. Wesemael, Astrophys. J. 443, 764 (1995), and D. Saumon, S. Jacobson, Astrophys. J. Lett. 511, L107 (1999).

10. B. Hansen, Astrophys. J. 520, 680 (1999).

11. A. Borysow, L. Frommhold, Astrophys. J. Lett. 348, L41 (1990).

12. N. C. Hambly et al., in preparation (2001).

13. N. C. Hambly, M. J. Irwin, H. T. MacGillivray, in preparation (2001).

14. N. C. Hambly, A. C. Davenhall, M. J. Irwin, H. T. MacGillivray, in preparation (2001).

15. N. C. Hambly, in The 12th European Conference on White Dwarfs, Newark, Delaware, 12 to 16 July 2000, H. Shipman and J. Provencal, Eds. (Astronomical Society of the Pacific, San Francisco, 2001).

16. D. G. Monet et al., Astron. J. 120, 1541 (2000).

17. W. J. Luyten, NLTT Catalog (Univ. Minnesota Press, Minneapolis, 1979).

18. R. A. Knox, M. R. S. Hawkins, N. C. Hambly, Mon. Not. Royal Astron. Soc. 306, 736 (1999).

19. C. Flynn, J. Sommer–Larsen, B. Fuchs, D. S. Graff, Mon. Not. Royal Astron. Soc., 322, 553 (2001).

20. R. D. Scholz, M. J. Irwin, R. Ibata, H. Jahreiß, O. Yu. Malkov, Astron. Astrophys. 353, 958 (2000).

21. N. C. Hambly, S. J. Smartt, S. T. Hodgkin, Astrophys. J. Lett. 489, L157 (1997).

22. N. C. Hambly et al., Mon. Not. Royal Astron. Soc. 309, L33 (1999).

23. P. Bergeron, M.-T. Ruiz, S. K. Leggett, Astrophys. J. Suppl. Ser. 108, 339 (1997).

24. S. K. Leggett, M.-T. Ruiz, P. Bergeron, Astrophys. J. Lett. 497, L294 (1998).

25. M. Chiba, T. C. Beers, Astron. J. 119, 2843 (2000).

26. P. F. L. Maxted, P. R. Marsh, Mon. Not. Royal Astron. Soc. 307, 122 (1999).

27. M. A. Wood, T. D. Oswalt, Astrophys. J. 497, 870 (1998).

28. P. Bergeron, R. A. Saffer, J. Liebert, Astrophys. J. 394, 228 (1992).

29. A. Gould, C. Flynn, J. N. Bahcall, Astrophys. J. 503, 798 (1998).

30. E. I. Gates, G. Gyuk, M. S. Turner, Astrophys. J. Lett. 449, L123 (1995).

31. G. Fontaine, P. Brassard, P. Bergeron, Publ. Astron. Soc. Pac., 113, 409 (2001). See also C. M. Tamanaha, J. Silk, M. A. Wood, D. E. Winget, Astrophys. J. 358, 164 (1990).

32. T. Lasserre et al., Astron. Astrophys. 355, L39 (2000).

33. P. Bergeron, S. K. Leggett, M.-T. Ruiz, Astrophys. J., in press.

34. We thank M. Irwin, I. King, G. Chabrier and T. Nakajima for helpful discussions. B.R.O. is supported by a Hubble Postdoctoral Research Fellowship, grant number HF-01122.01-99A. D.S. acknowledges support from NSF Grant AST97-31438.

Note added after Science Express online publication: This text reflects the correction of a few typographical errors in the online version of the table. It also includes the new constraint on the calculation of $`d_{\mathrm{max}}`$ which accounts for the fact that the survey could not have detected stars with proper motions below 0.33 arcseconds per year.

16 February 2001; accepted 13 March 2001

Published online 22 March 2001;

10.1126/science.1059954

Include this information when citing this paper.


text:
# Quasiparticle Liquid in the Highly Overdoped Bi2Sr2CaCu2O8+δ

## Abstract

We present results from the study of a highly overdoped (OD) Bi2Sr2CaCu2O8+δ with a $`T_c=51`$K using high resolution angle-resolved photoemission spectroscopy. The temperature dependent spectra near the ($`\pi ,0`$) point show the presence of the sharp peak well above $`T_c`$. From the nodal direction, we make comparison of the self-energy with the optimally doped and underdoped cuprates, and the Mo(110) surface state. We show that this OD cuprate appears to have properties that approach that of the Mo. Further analysis shows that the OD has a more $`k`$-independent lineshape at the Fermi surface than the lower-doped cuprates. This allows for a realistic comparison of the nodal lifetime values to the experimental resistivity measurements via Boltzmann transport formulation. All these observations point to the validity of the quasiparticle picture for the OD even in the normal state within a certain energy and momentum range.

The question of whether the Fermi Liquid (FL) model is a valid description for the high-$`T_c`$ superconductors has been one of the central issues in condensed matter. For the optimally doped (OP), and certainly the underdoped (UD) cuprates, there are many results which indicate that the FL picture may not be valid especially in the normal state. It is speculated that as doping increases beyond the optimal value, the cuprate may become a FL. However, experimental evidence from the highly overdoped (OD) regime showing FL properties is scant. Moreover, the FL term has become imprecise in recent literature. In the strictest definition of a FL, transport properties and the self-energy, $`\mathrm{\Sigma }`$, of the single-particle excitations or quasiparticles (QP) are proportional to the square of binding energy, $`\omega `$, and temperature, $`T`$, reflecting the electron-electron interactions. For regular metals, this is not the case for most $`T`$ and $`\omega `$ because the electron-phonon interaction dominates the scattering rate and the resistivity. At the other extreme, the term FL has been used to describe any material that has well defined single-particle excitation regardless of the nature of the QP scattering or transport. In this Letter, we describe in detail the nature of $`\mathrm{\Sigma }`$ of the highly OD cuprate Bi2Sr2CaCu2O8+δ (Bi2212). We find a state where there are well-defined single-particle excitations but the details do not match either a strict definition of a FL nor does $`\mathrm{\Sigma }`$ appear to be dominated by the electron-phonon interaction as found in normal metals. We use the term “quasiparticle liquid” (QPL) as a way to describe such a system that has: (i) well-defined single-particle excitations around the Fermi surface, and (ii) excitations that govern the transport properties even if these properties do not have the classic $`\omega ^2`$ and $`T^2`$ behavior or electron-phonon-like behavior.

Angle-resolved photoemission spectroscopy (ARPES) is, in principle, the best method to determine the nature of the QP in a 2D solid. The presence of a sharp peak in the spectra potentially indicates the presence of a long-lived QP. Indeed, this has been used to argue for the presence of QP for OP and UD cuprates below $`T_c`$ However, the existence of a sharp peak is not in itself proof of the existence of a FL. In the normal state of the OP and UD, there may be a sharp peak but only in the nodal direction. In the superconducting state, the well-defined sharp peak appears around the Fermi surface but only for a very short range of $`\omega `$ and $`k`$ and appears to be unrelated to transport measurements. In fact, for most $`\omega `$ and $`k`$ values, the spectral function blows up into a broad, ill-defined excitation. Analysis of the peak in the nodal direction for OP cuprates reveals that $`\mathrm{\Sigma }`$ is unlike that for a typical metal. In model 2D metals, the extracted $`\mathrm{\Sigma }`$ values are amenable to analyses in terms of traditional scattering mechanisms. More exotic phenomenologies have been invoked to understand similar data for the OP. Our approach is to not only present experimental data and analysis for highly OD Bi2212 but also to compare the results to the lower-doped compounds and to a metallic 2D surface state of Mo(110). This allows us to determine crucial differences amongst cuprates of different dopings, and between the cuprates and other real materials.

Single crystals of Bi2212 were synthesized using the floating-zone technique that yielded optimal $`T_c`$ onset of 95 K. Overdoping is achieved by annealing the samples in oxygen at 400 C for 5 days, yielding $`T_c`$ 51K. The samples were kept in the cells under oxygen pressure until they were mounted in the vacuum chamber. The crystals were cleaved in situ at $`T`$ below 150 K under vacuum with a base pressure of $`8\times 10^{11}`$ Torr. Crystal orientation was determined using low-energy electron diffraction (LEED). All data shown here were obtained within 12 hours after cleaving to minimize variation in doping content of the exposed surface with time. $`T_c`$ of the samples were determined to be $``$51K $`\pm `$ 5K using a SQUID magnetometer. All ARPES measurments were performed at beamline U13UB at the NSLS with photon energy of 22 eV and a Scienta SES 200 hemispherical analyzer that simultaneously collects a large energy (0.5 to 1 eV) and angular (120) window. This allows for the simultaneous collection of the photoemission intensity as a function of binding energy (energy distribution curve or EDC) and momentum (momentum distribution curve or MDC) along a particular direction in the BZ. The resulting energy resolution is $``$10 meV with an angular resolution of better than 0.20.

Fig. 1a shows the $`T`$ dependence of the EDCs near the M-point ($`\mathrm{\Gamma }`$M direction) of the BZ. The spectra show a well-defined sharp peak along with a broader hump that persist well into the normal state, a behavior that was also seen in a similar OD material. In OP and UD cuprates the sharp peak disappears around $`T_c`$. A recent report shows the sharp peak disappearing at slightly above $`T_c`$ on lightly OD crystals and claims the peak intensity is a measure of the superfluid density. This cannot be the case for our more OD samples since the sharp peak persists to temperatures well above $`T_c`$ where there is no superfluid. The only change in lineshape as $`T`$ goes above $`T_c`$ is a shift in the leading edge to signify the closing of the gap, as shown in Fig 1b. The gap size estimated from this shift is 10 meV, which is consistent with that observed in tunneling measurements on similar crystals..

A visual demonstration of the different energy ranges over which the ARPES excitation is well defined has been shown in Ref.  for the Bi2212 cuprates and Ref.  for Mo(110). When compared to the lower-doped cuprates, the OD ARPES spectra show a larger range of $`\omega `$ and $`k`$ over which the spectral peak is well defined. $`k_F`$ for the OD in this direction is measured to be 0.39A-1 or 0.34($`\pi ,\pi `$), compared to 0.446A-1 or 0.391($`\pi ,\pi `$) for OP Bi2212, indicating a shift of the Fermi surface towards the $`\mathrm{\Gamma }`$ point. In contrast to these samples, the Mo shows well defined peaks over the entire band. Clearly the $`\omega `$ range over which the ARPES spectra is well defined is much greater for the OD sample versus either the UD or OP, but still less than that for the Mo.

The width $`\mathrm{\Delta }k`$ of the MDC spectra in the $`\mathrm{\Gamma }`$Y region is used to obtain Im$`\mathrm{\Sigma }`$ using the relation $`\mathrm{}v_k\mathrm{\Delta }k2`$Im$`\mathrm{\Sigma }`$, where $`v_k`$ is the non-interacting band velocity. Fig. 2 shows Im$`\mathrm{\Sigma }`$ along with similar data for the OP and the Mo. None of these data sets show Im$`\mathrm{\Sigma }`$ varying as $`\omega ^2`$ for small energies. As noted previously, the Mo data can be broken into contributions from three scattering mechanisms: electron-electron, electron-phonon, and defect scattering. The phonon mechanism dominates in the vacinity of $`E_F`$. Unlike Mo, 2Im$`\mathrm{\Sigma }`$ for the OD and the OP Bi2212 have a predominantly linear dependence on $`\omega `$. For the OP sample, the slope of Im$`\mathrm{\Sigma }`$ versus $`\omega `$ is much greater and follows a characteristic Marginal Fermi Liquid (MFL) behavior where Im$`\mathrm{\Sigma }`$ max($`\omega ,T`$). Such an identification is less clear in the OD sample, but appears possible. Certainly the OD does not seem to be easily analyzed in terms of the typical scattering processes in metals. Overall, the peak broadening that gives Im$`\mathrm{\Sigma }`$ for the OD sample lies in between the behavior of the OP and the Mo.

In both the normal and superconducting state, $`\mathrm{\Sigma }`$ of the OD cuprate is much less $`k`$ dependent than for the lesser-doped compounds. This is manifested in the ARPES lineshape at a given $`\omega `$ at different parts of the BZ. Fig. 3a and b compare the OD EDC’s from the two in-plane symmetry directions at several binding energies, both below and above $`T_c`$. There is a remarkable similarity in the lineshape along $`\mathrm{\Gamma }`$M and $`\mathrm{\Gamma }`$Y. The main difference is that the remnant broad hump seen for $`\omega _p`$ 12 meV is only present in the $`\mathrm{\Gamma }`$M spectra. This hump becomes less pronounced at higher $`\omega `$ where the spectra look more similar to each other. We compare this to the spectra of OP Bi2212 (Fig. 3d) where there is a dramatic difference between the spectra from the two directions. This is especially true above $`T_c`$ where there is no well defined peak in the $`\mathrm{\Gamma }`$M spectrum. The difference in lineshapes is even stronger in UD Bi2212. We conclude that well-defined ARPES peaks occur over a much greater $`k`$ range in the normal state of the OD compound versus the OP or UD material. In particular $`\mathrm{\Sigma }`$ seems to be $`k`$ dependent in the UD and OP samples but largely $`k`$ independent in the OD material. It could be argued that with increasing hole doping, the influence of the antiferromagnetic “environment” in the system becomes weaker. This is consistent with the weaker coupling strength with increasing doping observed in the nodal direction.

In conventional metals, transport is governed by QP that scatter due to interactions with other excitations. The resistivity $`\rho `$ in these materials are related to the lifetime or the mean free path via Boltzmann transport theory that involves an integral over the Fermi surface. Although the lifetimes of ARPES QP, and thus ARPES peak widths, are not identical to that obtained in transport measurements, they should be roughly governed by the same scattering processes. This peak width is directly measured in ARPES via $`\mathrm{\Delta }k`$, the width of the MDC peak at $`E_F`$. For the OD samples, we have QP that are well defined over the whole Fermi surface and scattering rates that are largely $`k`$-independent as seen in Fig. 3. Thus, it is realistic to consider that in the OD case, the Boltzmann transport formulation would produce a meaningful comparison between the lifetime and the measured in-plane resistivity $`\rho _{ab}`$, and that the lifetime or mean free path obtained in the nodal direction can be taken as the average value. Such consideration results in the general relationship where $`\rho _{ab}`$ is directly proportional to $`\mathrm{\Delta }k`$. Figure 4 shows a comparison of the ARPES $`\mathrm{\Delta }k`$ data with $`\rho _{ab}`$. The resistivity data comes from a differently-prepared crystal that is overdoped to have roughly the same $`T_c`$, the best comparison available. We find that both $`\rho _{ab}`$ and $`\mathrm{\Delta }k`$ have the same functional form. Note that $`\mathrm{\Delta }k`$ appears to be insensitive to the superconducting transition. The similarity in the temperature dependence of the normal state ARPES $`\mathrm{\Delta }k(E=E_F)`$ and $`\rho _{ab}`$ is a strong evidence that both transport and ARPES are dominated by quasiparticles with similar scattering interactions.

Recent ARPES studies of OD cuprates have shown the possiblity of two separate bands arising out of the bilayer splitting. While our data (not shown) may indicate evidence for a two-component Fermi surface, the specific labeling of the hump feature in the EDC spectra near ($`\pi ,0`$) as the bonding band is inconsistent with our data in terms of the $`T`$ dependence and trends of relative peak intensities with doping as shown in Fig 1. We note that the mere presence of separate bilayer split bands has little impact on the conclusions of this paper. In fact, the idea of interplane coupling appearing in the OD samples is understandable in terms of the results reported here. Since the interplane coupling is dominated by states near the M-points of the BZ, the presence of more well-defined QP in those regions for the OD would increase the interplane coupling that possibly leads to the bilayer splitting. Furthermore, the well-defined QP would enhance the coherent interlayer tunneling, resulting in a greater $`c`$-axis conductivity. This has been shown to be the case where the $`c`$-axis resistivity decreases and appears to be metallic with increasing doping.

Taken as a whole, there are several elements of the ARPES data which indicate that the OD Bi2212 may be considered a QPL in the normal state, though the primary interactions are considerably different than those in a prototypical metal. The key elements arguing for the QPL nature of the normal state of the OD sample are: (i) the presence of well-defined peaks in the ARPES spectra over a wide range of $`T`$, $`\omega `$ and $`k`$, along with well defined peaks for the entire Fermi surface, (ii) $`\mathrm{\Sigma }`$ that is largely $`k`$ independent, unlike the UD and OP samples, and (iii) $`\rho _{ab}`$ that is proportional to the mean free path measured in ARPES as is expected from simple Boltzman transport. The lower-doped cuprates do not seem to be well described as a QPL in the normal state. The fundamental interactions between the electrons and other excitations is similar for all of the cuprates, and not like that for a typical metal. The strongest evidence for this is the linear $`\omega `$ dependence of Im$`\mathrm{\Sigma }`$ as shown in Fig. 3 which do not reflect contributions from several different scattering processes and no sign of contributions from the electron-phonon interaction. One alternative scenario involves some form of coupling to the magnetic modes observed in inelastic neutron scattering. The overall changes in the ARPES $`\mathrm{\Sigma }`$ as a function of $`T`$ and doping can be matched to the magnetic behavior. However, regardless of the details of the interactions, there is a clear difference between the OD and OP cuprates in both the presence of QP and their relation to transport properties.

Work supported in part by Dept. of Energy under contract number DE-AC02-98CH10886, DE-FG02-00ER45801, DOE-BES W-31-109-ENG-38, and in part by the New Energy and Industrial Technology Development Organization, Australia. The authors acknowledge discussions with John Zasadzinski, Liam Coffey, Phil Allen, Donglai Feng, and Bob Gooding.

FIG. 1. (color) EDCs for OD Bi2212 ($`T_c=51`$K) near the M-point (see inset) at various temperatures. The data in Fig. 1a have been shifted vertically for clarity. Six representative spectra are shown unshifted in Fig. 1b.

FIG. 2. 2Im$`\mathrm{\Sigma }`$ in the $`\mathrm{\Gamma }`$Y direction (see inset) for OD ($`T_C=51`$K), OP ($`T_C=91`$K), and Mo(110) surface state. The values shown in the legend are the temperatures where the ARPES spectra were obtained. 2Im$`\mathrm{\Sigma }`$ for the OD was obtained exclusively from the MDC spectra, whereas the OP values were obtained from both MDC and EDC. The scattering due to impurities has been subtracted from the Mo data. The straight lines are guides to the eye.

FIG. 3. (color) EDC lineshapes for the OD (a,b, and c) and OP (d) taken at the M-point and the nodal direction. Each data have been normalized to the intensity of the sharp peak, except for the normal state of OP where it has not been normalized. $`\omega _p`$ is the binding energy of the peak intensity.

FIG. 4. A comparison of the ARPES $`\mathrm{\Delta }k`$ (open circles) and experimental $`\rho _{ab}`$ (solid line) from OD Bi2212 with $`T_c=58`$K..