text:
# The impact of unresolved binaries on searches for white dwarfs in open clusters

## 1. Introduction

White dwarfs (WDs) are the final endpoint of stellar evolution for the vast majority of stars. WDs are gaining importance in tracing the history of stellar populations in the galaxy, including the galactic disk (e.g. Winget et al., 1987), open star clusters (e.g. Richer et al., 1998), and globular star clusters (Hansen et al., 2002). WDs are also useful in studies of supernova physics, as the upper mass of WD progenitors represents the critical progenitor mass ($`M_{crit}`$)111Also referred to as $`M_w`$, $`m_w`$, and $`M_{up}`$ for core-collapse supernovae. The value of $`M_{crit}`$ is relatively uncertain, with the best estimates of $`5.5M_{\mathrm{}}M_{crit}10M_{\mathrm{}}`$.

The first significant population of white dwarfs (WDs) in open clusters was identified by Eggen & Greenstein (1965) in the Hyades (catalog ), Praesepe (catalog ), and the Pleiades (catalog ). Tinsley (1974) compared the number of Hyades WDs to the estimated number of Hyades stars having completed their lifetimes. Her work found that the WD numbers are compatible with $`M_{crit}4M_{\mathrm{}}`$, though this determination depends sensitively on the assumed initial-mass function (IMF), the turnoff mass of the Hyades, and the Hyades distance modulus. This value of $`M_{crit}`$ is lower than current lower limits ($`M_{crit}6M_{\mathrm{}}`$), set in large part by the existence of LB 1497 (catalog ), the lone WD in the Pleiades (catalog ) (turnoff mass $`5.4M_{\mathrm{}}`$).

Weidemann et al. (1992) discuss the Hyades WD population and determine a deficit of 21 WDs (28 WDs predicted versus 7 observed WDs) for $`M_{crit}=8M_{\mathrm{}}`$ and ascribe this white dwarf deficit to dynamical evaporation of WDs from the Hyades. Similar deficits of WDs have since been observed in other open clusters, including M67 (catalog ) (Richer et al., 1998), Praesepe (catalog ) (Claver et al. 2001, hereafter CLBK ), and NGC 2099 (catalog ) (Kalirai et al., 2001c). The existence of the white dwarf deficit is still debatable, as discussed by von Hippel (1998).

Three explanations for the white dwarf deficit in open clusters readily come to mind. First, WDs may evaporate from open clusters due to dynamical evolution. Mass segregation and galactic tidal fields alone are not sufficient to remove WDs from open clusters, as a 0.6$`M_{\mathrm{}}`$ WD is still more massive than typical cluster stars and thus unlikely to suffer preferential evaporation. This is borne out by modern $`N`$-body simulations of open clusters, such as those of Portegies Zwart et al. (2001); Baumgardt & Makino (2003) and Hurley & Shara (2003), which find that WDs remain bound in open clusters. However, if a WD receives a sufficient velocity kick from asymmetric mass loss during its post-main sequence evolution, the WD may become unbound from the open cluster. This scenario was first suggested by Weidemann et al. (1992) to explain the WD deficit in the Hyades (catalog ). Recent, simple $`N`$-body simulations by Fellhauer et al. (2003) confirm that this mechanism can preferentially remove WDs from an open cluster, though scant observational evidence exists.

A second explanation for the white dwarf deficit, though not exclusive of the dynamical evolutionary argument, is that the WDs may be hidden in binary systems in which the intrinsically faint WDs are not detected due to the overwhelming light of brighter companion (e.g. Kalirai et al., 2001c). Searches for WDs in binaries have been undertaken in the Hyades. Two of the known Hyades WDs, EGGR 38 (catalog ) and V471 Tau (catalog ), are in binary systems. More recently, additional Hyades WDs have been discovered hidden in unresolved binaries, including HD 27483 (catalog ) (Böhm-Vitense, 1993), VA351 (catalog ) (Franz et al., 1998), and four potential WDs in Am binary star systems (Debernardi et al., 2000). Clearly some WDs lie hidden in unresolved binary systems, but the exact numbers are not known.

The third explanation for the small number of white dwarfs seen in the Hyades and other well studied young clusters could be that our expectations of WD numbers are incorrect. The deficit discussed for the Hyades would go away if $`M_{crit}`$ is low ($`4M_{\mathrm{}}`$), if the IMF for masses above the cluster turnoff mass is steeper than the present-day mass function around the turnoff mass, or a combination of these two effects. The low value of $`M_{crit}`$ seems unlikely, given that WDs have been observed in clusters with turnoff masses higher than 4$`M_{\mathrm{}}`$, including the Pleiades, NGC 2516 (catalog ) (Koester & Reimers, 1996) and NGC 2168 (catalog ) (Reimers & Koester, 1988).

We have developed a Monte Carlo method of calculating the number of WDs detectable in observations of specific open clusters, a calculation designed to aid in the interpretation of data from the Lick-Arizona White Dwarf Survey (LAWDS, Williams & Bolte, 2003) and could be used in conjunction with other ongoing cluster WD surveys, such as the CFHT Open Cluster Survey (Kalirai et al., 2001a) and the WIYN Open Cluster Study (von Hippel & Sarajedini, 1998). This calculation utilizes the observed characteristics of specific open clusters, including the cluster age, distance and reddening, and determines how many WDs would be detected given photometric WD selection criteria. The structure of this calculation, its limitations, and tests of the calculation are presented in §2. In §3 we present the results of the calculations for the Hyades, Praesepe, and the Pleiades and compare these results to the observed WD populations. In §4 we discuss the results and discuss the usefulness of these calculations in regard to current searches for WDs in open clusters.

## 2. Calculating the WD populations in open clusters

The number of WDs detectable by photometric observations of open clusters can be calculated by appropriate modeling of the cluster, including the effects of unresolved binary stars, choice of IMF, WD detection criteria, $`M_{crit}`$, metallicity and cluster age. Assuming that such an open cluster model is sufficiently realistic, any significant difference between the calculated numbers of observable WDs and the actual numbers of observed WDs represents the “true” white dwarf deficit that must be explained by other means, such as evaporation or other dynamical evolution.

Monte Carlo techniques are used to calculate the observable WD population. Stars are randomly drawn from an input IMF. For this work, four input IMFs $`\xi (M)M^{(1+\mathrm{\Gamma })}`$ were considered: a Salpeter IMF with $`\mathrm{\Gamma }=1.35`$ (Salpeter, 1955), a steeper power law that accounts for unresolved binary systems with $`\mathrm{\Gamma }=2`$ (Naylor et al., 2002), and a broken power law with a flat low-mass slope: $`\mathrm{\Gamma }=0.2(M1M_{\mathrm{}});\mathrm{\Gamma }=1.8(M>1M_{\mathrm{}})`$ (Naylor et al., 2002), hereafter called the “Naylor IMF,” and the Kroupa IMF with $`\mathrm{\Gamma }=0.3(M<0.5M_{\mathrm{}});\mathrm{\Gamma }=1.3(M0.5M_{\mathrm{}})`$ Kroupa (2001).

Stars are assigned a binary companion with a probability based on an input binary fraction; the binary mass fraction is assumed to be random, i.e. the primary and secondary stars are drawn from the same IMF (see §2.1 below). The stars are then evolved to the input cluster age using the $`Z=0.008,Z=0.019`$ or $`Z=0.040`$ solar-scaled metallicity stellar evolutionary models of Salasnich et al. (2000). These models provide evolutionary data on stars with masses up to 20$`M_{\mathrm{}}`$, sufficiently massive to explore all likely values of $`M_{crit}`$.

Two outcomes are assumed for stars having completed their evolution. Stars with $`M_{}>M_{crit}`$ undergo supernovae. Stars with $`M_{}M_{crit}`$ form a WD. The WD mass is determined from a linear approximation of the initial-final mass relation in Weidemann (2000),

$$M_{WD}=\frac{3}{35}M_{}+0.414M_{\mathrm{}}.$$
(1)
Each WD is assigned a cooling age equivalent to the progenitor lifetime subtracted from the input cluster age. The radii, effective temperature, and surface gravities of each WD are interpolated from the carbon/oxygen models of Wood (1995, hereafter W (95)).

Photometric indices for each system are determined by summation of the flux from each component. Supernova remnants are assumed to emit no flux. WD fluxes are determined from the updated model atmospheres of Finley, Koester, & Basri (1997) graciously provided by D. Koester. Stellar fluxes are taken from the isochrones of Girardi et al. (2002) for the Salasnich et al. (2000) evolutionary models. These indices are corrected for the input distance modulus and reddening, assuming the interstellar reddening law of Rieke & Lebofsky (1985) and $`R_V=3.1`$.

Once the photometric indices for each system have been calculated, these indices are compared to the input WD photometric-selection criteria. These criteria include color indices and single-band magnitude limits, with the intention that the photometric criteria used to select WD candidates in actual observations of a given cluster are used in the simulation of that cluster. If a system meets these input photometric selection criteria, the number of detected WDs is incremented by one. If it fails to meet the selection criteria, the WD is considered undetected.

The Monte Carlo calculation continues to select stellar systems until the number of stars brighter than a limiting magnitude $`V_{lim}`$ reaches the input value. For example, if the cumulative luminosity function of an actual open cluster has $`N_{}`$ stars brighter than $`V_{lim}`$, then each simulated cluster will also contain $`N_{}`$ stars brighter than $`V_{lim}`$. $`10^4`$ realization of each cluster are used to ensure high precision in the output.

### 2.1. Limitations of the calculations

These calculations represent a simplified version of actual open clusters. The simulation assumes that all binary stars are unresolved, an assumption whose validity depends strongly on the cluster distance and the distribution of physical separations among binary stars. Binary interactions during stellar evolution, such as mass transfer systems, common envelope stages, cataclysmic variables, and the potential disruption of binaries due to mass loss are ignored. Open cluster dynamical evolution, such as evaporation of stars and evolution of the binary fraction, is also ignored. The binary mass ratios are assumed to be random, which is likely not be the case for close binaries, but is reasonable for wider binaries (Mason et al., 1998; Abt, Gomez, & Levy, 1990). Finally, stellar systems with more than two components are not considered.

Ignoring binary-star evolution and the different mass ratios of close binaries than a random IMF results in far fewer double-degenerate systems being created than as seen in realistic $`N`$-body simulations, such as the significant double-degenerate population seen in Hurley & Shara (2003). The presence of any observed double-degenerate sequence in an open cluster therefore can provide a constraint on the relative importance of close-binary evolution in open clusters.

The lack of consideration of dynamical evolution is not a deficiency in the calculations, but rather is one purpose of these calculations – to determine the actual number of missing WDs that dynamical evolution must be invoked to explain. The strong, observed evaporation of low-mass stars has no great impact on this calculation. Low-mass stars in binaries with a (comparatively) high-mass WD should not be affected greatly by evaporation processes, and low-mass stars not in a binary system with a WD are not used in the calculation beyond the initial selection of cluster stars.

In short, these calculations are not a simulation of WDs in the larger sense of open cluster evolution. Rather, this calculation provides a more sophisticated means of comparing the observed WD populations of specific open clusters to the expected observable WD population. Any differences between the actual observed WD population and the calculated observable population provide a measure of the strength of the dynamical evolution of the WD population.

### 2.2. Trials of the calculations

We ran numerous tests to verify that these calculations return reasonable results for certain limiting cases. Fig. 1 shows the output IMF of a single simulated cluster with 5000 stars more massive than 1$`M_{\mathrm{}}`$ for each of the input IMFs. As one would hope, the output IMFs agree with the input IMFs.

Simulated clusters with a $`\mathrm{\Gamma }=2`$ power-law or the Naylor IMF are found to contain fewer WDs than those with a Salpeter or Kroupa IMF, merely due to the fact that fewer high-mass stars formed due to the steeper slope. The Naylor and Kroupa IMFs result in a higher fraction of WDs being hidden in binary systems than the single power-law IMFs due to the relatively flat low-mass IMF slope, as a companion to a WD is more likely to be higher mass and therefore more luminous than the WD.

Fig. 2 shows output color-magnitude diagrams (CMDs) for calculations of clusters at four different ages containing no binary stars. Occasional minor numerical issues arise during the interpolation of photometric indices for evolved stars due to the large changes in magnitudes with minute changes in stellar mass. These issues arise primarily at inflection points in the isochrones, such at the helium flash and at the tip of the AGB. These errors do not affect the simulated detection of WDs, as both the correct and erroneous magnitudes for these stars are bright enough to hide any WD companions, and the number of stars affected is very small.

The effect of the addition of binary stars is shown in Fig. 3. The sequence of binaries consisting of two main-sequence (MS) stars is readily visible in the $`\mathrm{B V},M_V`$ CMD. The fan of points leading from the MS to the WD cooling sequence is composed of MS-WD binaries. The $`\mathrm{B V},\mathrm{U B}`$ color-color diagram illustrates the advantage of three-band photometry over two-band photometry. The WD-MS binaries, which would blend in with background disk MS stars in a $`\mathrm{B V},V`$ CMD, are well-separated from the locus of main-sequence in color-color space. In addition, the color-color space information clearly separates the WD cooling sequence from massive MS stars.

A WD-WD binary sequence is visible only when the number of stars in a cluster is very high, as seen in the lower-right panel of CMD in Fig. 3. Any significant observed WD-WD binary population in open clusters would therefore confirm that the binary component mass ratios are not random and/or that dynamical processes such as binary interactions are not negligible. Hurley & Shara (2003) find evidence of a WD-WD sequence in the photometry of NGC 6819 of Kalirai et al. (2001b).

Multiple ($`10^4`$) realizations of each cluster are used to calculate precise means and standard deviations for various categories of WDs, including the total number of WDs, the total observed number of WDs, the number of WDs in binary systems, etc.

## 3. Simulations of specific clusters

We have calculated the predicted observed WD populations for three open clusters: the Hyades, the Praesepe, and the Pleiades. The output populations are compared to the number of observed WDs in each system. The calculated numbers of WDs in each cluster fulfilling various criteria, including whether the WD met the photometric selection criteria (these are given individually for each cluster below), whether a WD was in a binary system, and WD progenitor masses were calculated, and the mean and standard deviation of each over the ensemble of trials were calculated.

The standard deviations are found to be nearly equal to those expected from Poisson statistics. Therefore, we use Poisson statistics to compare the calculated populations to observations of the actual cluster WDs. The probability $`P`$ of obtaining the actual number of observed WDs in a category or fewer is calculated from Poisson statistics given the mean value obtained from the Monte Carlo calculation. If $`P0.05`$, then the derived WD population is considered to be inconsistent with the observed population. For probabilities $`0.15P>0.05`$, the calculated is termed to be mildly inconsistent with the data.

*The Hyades* — The Hyades (catalog ) is the open cluster with the best-studied WD system. Nine WDs are traditionally considered cluster members, though some questions as to the membership of some of these have been raised (Weidemann et al., 1992). Of these, the seven lone WDs and the WD-MS binary EGGR 38 (catalog ) are detectable as blue excess objects in $`UBV`$photometry.

For the Monte Carlo realizations the Hyades, the following physical parameters are used. The intrinsic distance modulus is assumed to be 3.33 (Perryman et al., 1998), the color excess $`E(\mathrm{B V})=0.01`$ (Mermilliod, 1995), and the logarithmic age is assumed to be 8.80 (Perryman et al., 1998). The binary star fraction is chosen to be 0.4, approximately the lower limit of $`0.41\pm 0.05`$ determined from the combination of speckle, spectroscopic, and direct imaging surveys (Patience et al., 1998).

The Hyades WDs have been detected and confirmed through a variety of photometric, spectroscopic, and astrometric surveys, and a direct comparison of these selection criteria with the calculated population is not performed. In order to determine if a cluster WD is to be considered observable, we apply the photometric criteria of $`\mathrm{U B}0.0`$ and $`\mathrm{B V}0.6`$ to select WDs. These values encompass the range of indices observed in the WD sample of Eggen & Greenstein (1965), including the seven single Hyades WDs and the binary EGGR 38 (catalog ). The MS-WD binaries not selected via these criteria represent the contribution of binaries to the WD deficit.

The metallicity of the Hyades is super-solar, with $`[\mathrm{Fe}/\mathrm{H}]=0.14\pm 0.05`$ (Perryman et al., 1998), corresponding to $`Z=0.034`$ for a solar $`Z=0.0245`$ (Christensen-Dalsgaard et al., 1996). We therefore calculate WD populations for each of two metallicities, $`Z=0.019`$ and $`Z=0.040`$, for comparison. The calculated population also explores three values of $`M_{crit}`$, 10$`M_{\mathrm{}}`$, 8$`M_{\mathrm{}}`$, and 6$`M_{\mathrm{}}`$, bracketing the range of likely values. Finally, the calculations were performed using each of the four IMFs discussed above.

Each realization of the cluster contains 138 stellar systems in the magnitude range $`3.0B10.0`$, consistent with the Hyades luminosity function of Oort (1979), a complete luminosity function of the Hyades upper main sequence. The resulting WD populations of the Hyades are given in Table 1.

Our results indicate that the eight photometrically-selected Hyades WDs are inconsistent ($`P0.03`$) with what would be expected from a Salpeter or Kroupa IMF, independent of the value of $`M_{crit}`$. However, the calculated number of WD detections for the $`\mathrm{\Gamma }=2`$ power-law IMF and the Naylor IMF are fully consistent with the observed number of WDs for $`Z=0.019`$, and consistent ($`P0.2`$) for $`Z=0.040`$. The lack of any observed Hyades WDs with progenitor masses $`4M_{\mathrm{}}`$ (Weidemann et al., 1992; CLBK, ) is mildly inconsistent for the $`\mathrm{\Gamma }=2`$ power law and Naylor IMFs at $`M_{crit}=6M_{\mathrm{}}`$ ($`P=0.13`$ and $`P=0.07`$, respectively) and for the $`\mathrm{\Gamma }=2`$ power-law IMF for $`M_{crit}=8M_{\mathrm{}}`$ ($`P=0.07`$), but is inconsistent with the survey results for the remaining explored parameter space. This result will be explored in the discussion (§4) below.

Of the eight photometrically-detectable Hyades WDs, one is in a binary system. This is inconsistent with the number of binaries detectable from a Salpeter IMF ($`P0.05`$) and mildly inconsistent with the Kroupa IMF and $`\mathrm{\Gamma }=2`$ power-law IMF for the binary fraction of 0.5, but is otherwise consistent with the calculations for the $`Z=0.019`$ systems. For $`Z=0.04`$, the single observed binary is inconsistent with one observed binary for each IMF except the Naylor IMF, which is fully consistent with the observations.

Böhm-Vitense (1995) reports on an IUE survey for WD companions to F stars in the Hyades, with only one WD-F star system (actually a triple system with two F stars and a WD) being detected. In our calculations, we have identified all stars meeting the B V selection criteria in that study. The average number of WD companions to these F stars are presented in the final column of Table 1. On average, $`1`$ F star-WD binary is present in a simulated cluster, fully consistent with the Böhm-Vitense (1995) results.

*Praesepe*CLBK  list ten candidate Praesepe (catalog ) WDs. Based on their analysis of six spectra, five are cluster WDs, and the other is a confirmed non-member. We obtained spectroscopy of two of the remaining candidates, WD 0834+209 (catalog ) and LB 1839 (catalog ), with the upgraded blue camera of LRIS (Oke et al., 1995) on the Keck 10-m telescope in December 2002. Both of these objects are found to be QSOs, with LB 1839 (catalog ) at $`z=0.86`$ and WD 0834+209 (catalog ) at $`z=1.09`$. Therefore, at most seven of the Praesepe candidates are WDs, and this is the number of Praesepe WDs used in the comparison to the calculated populations.

The realizations of the Praesepe calculations were normalized by the luminosity function of Jones & Stauffer (1991), who find 179 stars (completeness-corrected) with $`2.5M_V7.5`$. While the Jones & Stauffer (1991) data extend to fainter magnitudes, CLBK  find that the Padova isochrones do not match the Praesepe MS below this magnitude, and this mismatch would systematically bias our calculations.

The age and metallicity of the Praesepe are nearly indistinguishable from that of the Hyades (CLBK, ), so we adopt the Hyades logarithmic cluster age of 8.80 for Praesepe. As the metallicity was found to have only a small effect on the calculated observable number of WDs in the Hyades, we simulate Praesepe using only the Z=0.019 isochrones. The binary star fraction of the Praesepe is $`0.4`$ (Mermilliod & Mayor, 1999), so this value is used in the simulation.

For WD detection criteria, we require $`\mathrm{U B}0.0`$ and $`\mathrm{B V}0.6`$, limits encompassing the CLBK  WD sample. As in the Hyades, the form of the IMF and $`M_{crit}`$ were varied. The calculated WD populations are given in Table 2.

As with the Hyades, the observed WD numbers are inconsistent with the calculated output for a Salpeter and Kroupa IMFs. The numbers of observable WDs calculated using the Naylor law IMF are consistent with the observed number for $`M_{crit}=6M_{\mathrm{}}`$ ($`P=0.24`$) but mildly inconsistent for higher values of $`M_{crit}`$ ($`P=0.15`$ and $`P=0.11`$ for $`M_{crit}=8M_{\mathrm{}}`$ and $`M_{crit}=10M_{\mathrm{}}`$, respectively). The $`\mathrm{\Gamma }=2`$ power-law IMF simulation results are fully consistent with the observed numbers.

CLBK  detected one WD with a progenitor mass $`4M_{\mathrm{}}`$. This number is inconsistent with a Salpeter and Kroupa IMF for all values of $`M_{crit}`$ ($`P0.04`$). The simulated number of WDs with massive progenitors is inconsistent with the Naylor IMF for $`M_{crit}=10M_{\mathrm{}}`$ ($`P=0.05`$), but consistent with the one observed WD for $`M_{crit}=8M_{\mathrm{}}`$ and $`M_{crit}=6M_{\mathrm{}}`$ ($`P=0.17`$ and $`P=0.20`$, respectively). The $`\mathrm{\Gamma }=2`$ power law WD numbers are reasonably consistent with the observations for all values of $`M_{crit}`$.

None of the five confirmed Praesepe WDs are known to be in binary systems. If neither of the two remaining candidates is in a binary system, then the lack of binaries is inconsistent with the simulated results for a Salpeter or Kroupa IMF. The lack of binaries is mildly inconsistent ($`P=0.08`$) with the $`\mathrm{\Gamma }=2`$ IMF, and consistent ($`P0.19`$) with the Naylor IMF. It is emphasized that this inconsistency is between the lack of observed WD binaries in Praesepe and the number of simulated WDs in binaries which would be detected given the stated photometric selection criteria; this discrepancy does *not* include the WDs in binary systems with photometric indices outside the selection criteria.

*Pleiades* — The Pleiades (catalog ) has one known member WD, LB 1497 (catalog ). For this cluster, we have normalized the IMF such that each simulated cluster contains 107 stars in the magnitude range $`2.0V10.0`$, consistent with the luminosity function of the Pleiades from Jones (1970). We assume the Pinsonneault et al. (1998) distance modulus of 5.60 and reddening of $`E(\mathrm{B V})=0.04`$. The metallicity of the Pleiades is near solar, with $`[\mathrm{Fe}/\mathrm{H}]=0.034\pm 0.024`$ (Boesgaard & Friel, 1990), so we use the Z=0.019 isochrones. We adopt the Mermilliod et al. (1992) multiple-star frequency of 0.36 as the binary frequency. Age estimates of the Pleiades have shown scatter, with recent estimates between 100 Myr (Pinsonneault et al., 1998) and 125 Myr (Stauffer, Schultz, & Kirkpatrick, 1998). We therefore produce calculations for logarithmic cluster ages of 8.0 and 8.1.

The results are given in Table 3. For $`M_{crit}=8M_{\mathrm{}}`$ and $`M_{crit}=10M_{\mathrm{}}`$, the Salpeter and Kroupa IMFs are inconsistent with the one observed WD, regardless of the adopted cluster age. For $`M_{crit}=6M_{\mathrm{}}`$, the Salpeter and Kroupa IMFs are mildly inconsistent with the observed WD if the Pleiades have a logarithmic age of 8.1 ($`P=0.12`$), but are consistent for a logarithmic age of 8.0. The $`\mathrm{\Gamma }=2`$ and Naylor IMFs produce WD populations mildly inconsistent with the Pleiades for $`M_{crit}=10M_{\mathrm{}}`$ and for the $`M_{crit}=8M_{\mathrm{}},\mathrm{log}(\mathrm{age})=8.1`$ simulations. Otherwise, these two IMFs result in WD populations consistent with the single observed WD.

## 4. Discussion

A few observations can be made from the comparisons between the calculated WD populations and the observed populations presented above. The calculations consistently over-predict the number of observed, color-selected WDs and the number of color-selected WDs in binary systems. However, the discrepancy between the expected number of WDs and the observed number is lower than in previous observational studies due to the inclusion of binaries in the simulation. Therefore, preferential evaporation of WDs from these clusters due to dynamical evolution need not be as severe as suggested by Weidemann et al. (1992) but cannot be ruled out. It should also be noted that the binary fractions used in this work represent the high end of the range of commonly-quoted cluster binary populations. A lower actual binary fraction would result in a larger WD deficit.

The number of observed WDs can be brought into agreement with the calculated numbers via at least four different methods. First, the slope of the IMF could be steeper than the slopes used here. Second, the binary fractions may be higher, resulting in more WDs being hidden in unresolved binaries. Third, the binary mass ratio could be closer to unity than the random parings discussed here, which would result in more binary WDs being hidden in systems. Fourth, some sort of dynamical evolution may be removing WDs from the open clusters. From the simulation alone, it is not possible to differentiate between these scenarios.

As mentioned above, the lack of observed WDs in the Hyades and Praesepe with progenitor masses $`4M_{\mathrm{}}`$ is inconsistent with the calculated populations. A steeper high-mass IMF would explain the lack of these WDs, as would dynamical evolution. However, there is some evidence that binary stars in this mass range tend to have a fairly flat distribution of mass ratios, especially for orbital periods $`0.1\mathrm{yr}`$ (Abt, Gomez, & Levy, 1990; Mason et al., 1998). If this is true, WDs with massive progenitors may be more likely to be in binary systems with massive, unevolved stars than in binary systems with fainter, low-mass stars, and therefore more likely to go undetected.

There are several avenues of research which could shed some light on these issues. First, systematic, thorough searches for WD-mass companions to stars in the Hyades and Praesepe are needed to complete the census of WDs in each cluster, thereby determining whether or not binary systems can account for the majority of “missing” WDs in open clusters. Such a survey could include a combination of spectroscopic and high-resolution imaging (e.g. speckle and adaptive optics surveys) of these clusters, such as those published in Patience et al. (1998) and similar surveys. These surveys could also determine if the binary mass ratios for any close binaries are non-random, as opposed to the random ratios assumed in this work.

Also, WD searches in open clusters of a wide range of ages, such as the LAWDS survey we are currently undertaking, can provide evidence as to whether the deficit of massive WDs is increasing with time, as would be expected for dynamical evolution. In addition, detailed comparisons of clusters of similar ages but differing compactness and richness should also show differences in the white dwarf deficit if dynamical evolution of the WD population is significant.

The calculations presented here will be useful in interpreting data from the WD observational programs. Most importantly, the calculations permit us to estimate how strong the white dwarf deficit is in a cluster, given basic assumptions about the IMF and binary fraction. This provides a better estimate on the significance of any apparent deficit and the cluster-to-cluster variations in the deficit than could be gleaned from simpler WD population estimators, such as simple integration of the IMF.

The author wishes to thank Michael Bolte for many helpful comments regarding this paper and for his skillful advising throughout this project. Thanks are also due to Matthew Wood, William Mathews, Graeme Smith, James Liebert, Brad Hansen and Jasonjot Kalirai for helpful discussions regarding this work. Thanks to Chuck Claver for finding charts for the Praesepe WD candidates. White dwarf spectral models were graciously and freely provided by Detlev Koester, and white dwarf evolutionary models were graciously and freely provided by Matt Wood; these are greatly appreciated. The author would also like to thank the anonymous referee for insightful comments which led to improvements in this manuscript. Some of the data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. The author wishes to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.


format: text

text:
# A thought experiment on consciousness

## 1 Introduction

One of the most fundamental problem in dealing with conscious experiences and consciousness is the following: if I am able to completely describe the physical state of my brain (conceding that all the physics necessary to such description is already known), may I safely say to have completely described my mental state, my subjective experience too? The point is that my subjective experiences, like for example those of pain, joy or smell (generally referred to as qualia), seem not to get exhausted in a physical-functional description of my cerebral states, even in the most complete description we are able to imagine to. Actually, the description of the physical processes which take place in my brain, when I experience pain for example, seems not to be a complete description of my subjective experience of pain; at most, it seems to be only a complete description of the cerebral states of my brain during my pain experience.

In other words, it seems that a barrier, impassable to every physical theory, forbids any complete objective description of subjective experience, or, at least, every complete objective description of a subjective experience simply does not include the subjective experience itself. The objective description and the subjective experience seem to belong to different and ‘orthogonal’ dimensions, the outside and the inside.

What I have described above briefly summarizes the well-known Mind-Body Problem, the main ingredient of the philosophical investigations of the mind and a thorn in the side of physicalism, namely of those who believe in a complete reduction of consciousness to peculiar physical processes of the brain (for accessible and exhaustive reviews of the Mind-Body Problem see, for example, Nagel  and Chalmers ).

In this paper I provide an equivalent formulation of the Mind–Body Problem, which I will call the ‘brain duplication’ argument, and I will show that if we assume that consciousness is in any case a physical process which takes place in the physical world, in the most general sense of these terms, namely that it does not transcend physical reality, then the human brain, in giving rise to the mind, might be characterized by the astonishing property of an extensive interdependence with the whole surrounding physical world.

## 2 The ‘brain duplication’ argument

For the sake of thought experiment, let us suppose that we manage to create an exact, physically identical duplicate of my brain, as it is in this precise moment. Actually, it does not matter whether we do not known operatively how to do it. The point is that since at least a brain exists physically, then nothing forbids us to imagine an identical physical duplicate of it, as well as nothing forbids us to imagine an identical physical duplicate, atom by atom, of the sheet of paper on which these words are written, just for the fact that such sheet of paper exists, even if probably we will never be able to do the actual duplicate.

Now, just after the creation of this duplicate, how would my own consciousness react? As everyone of you can experience directly, one of the leading characteristics of consciousness is the perception of its own uniqueness, the uniqueness of oneself and of its own conscious experience, and, in addition, the perception of the persistence of such uniqueness (I feel to be myself also in different periods and different places). Therefore, if I were able to exactly duplicate my brain, would my own consciousness change? And if it changes, how does it change? Would I feel to be here, where I was before the duplication, but at the same time would I feel to be there, where my brain duplicate is now?

I believe that the most natural answer for everyone is that I will continue to be myself as I was before the duplication. But then, who or what is my exact physical duplicate? If the two brains are physically identical and their consciousnesses are different, in what do they differ? Here two possible approaches to the problem are presented.

## 3 Non-physical explanation

Consciousness is not physically reproducible in the reality. Namely it does not depend on physical reality and it is in some sense ‘outside’ it: thus, distinct conscious experiences and consciousnesses may even be attached to two physically identical brains, or only one of these brains may be conscious (and the other one may not; the duplicate brain, for example, might be a so-called zombie). If it is so, there is not much to do; as a matter of fact, consciousness would constitute a prime and alien property, to be added to the rest of the physical properties of the brain.

However, such hypothesis seems to be scientifically frustrating and, after all, not particularly reasonable. I guess that not many researchers would honestly feel up to deny any link between consciousness and the physical world (even if such link is not completely clear from a scientific point of view). Even if we assume that the existence of the brain is a necessary but not sufficient condition for the presence of consciousness, the latter must necessarily have a physical interaction with the brain, otherwise it would not even make sense to talk about brain as necessary condition, without speaking about the substantial amount of data achieved nowadays in the neurosciences on the neural correlates of consciousness. Moreover, matter can act on mind and consciousness, as any physical and chemical interferences on our state of consciousness can easily demonstrate. So, matter, mind and consciousness have to speak the same ‘language’. And therefore consciousness must be a physically characterizable, a physically tractable entity.

Such physical entity might be completely internal or external (partially or totally) to the brain. In the first case we just would have that the exact physical reproduction of one’s brain would give tout court the consciousness of that person (and hence, the brain duplication dilemma). If, on the contrary, this entity were external, it might be one and unique, or there might be many, one of them corresponding to each consciousness currently existing in the world. Yet, this further distinction is not important in our case: plausibly, two physically identical brains would always interact in the same way with the same external physical agent (like two equally tuned radio receivers ‘interact’ always with the same radio station, although there are a lot of different radio frequencies in the air), and again the brain duplication dilemma would be left untouched. However, I believe that now this dilemma is ready to be tackled with the possible physical explanation described in the next section.

## 4 A possible physical explanation

Let us suppose instead that consciousness depends on physical reality, namely on the peculiar organization of physical matter and energy, in the most general sense of these terms, and also let us do not exclude the possibility, suggested in the previous section, that consciousness originates not only in the physics of the brain but also through the interaction with an external physical process.

A possible explanation of the seemingly paradoxical picture originated in the brain duplication argument is that (the evolution of) the cerebral processes involved in the rise of consciousness might physically depend not only on the physics of the brain itself but actually also on all the things which physically surround such brain, in the sense that it depends on the organization of all the surrounding physical matter and energy.

In such framework, the physically exact copy of my brain, which is in the region of the space $`(x_c,y_c,z_c)`$, interacts with all the surrounding physical world, and thus also with my original brain which is in the region of space $`(x_o,y_o,z_o)`$. Similarly, my own original brain, which is in the region of space $`(x_o,y_o,z_o)`$, interacts with all the surrounding physical world, and thus also with the physically exact copy of my brain, which is in $`(x_c,y_c,z_c)`$. Therefore, it is clear (see fig. 1) that the boundary conditions for the two identical brains are not identical, and the evolution of cerebral states, as well as consciousness, might not be identical (although the memory of all the past experiences might be the same in every details for both brains). A similar conclusion on the possibility of an extensive physical interdependence of the human brain with the whole surrounding physical world in giving rise to mind has been drawn in another thought experiment (D’Abramo ) but through a different approach involving the notion of algorithmic complexity.

Note that any finite version of this approach, i.e. that the suggested physical dependence is only on a finite portion of the surrounding physical world, does not change much the point. As a matter of fact, the brain+‘finite portion of physical world’ must be a physically isolated system, i.e. a totally isolated system. For if it were not, two such systems within two different environments would eventually lead to the same dilemma on consciousness as before, since the same consciousness will eventually experience different environments at the same time. Hence, even in any finite approach, the finite portion of the surrounding physical world is in fact all the surrounding physical world. Now, whether there are spatial limits to be understood in the word “all ” used before, and whether they are posed by the concept of visible universe and finite speed of light, the discussion of this is beyond the scope of the present note.

Thus, the brain might be far from being a semi-closed box, opened only to the five senses, and the common perception that our brain does not continuously interact at a deeper physical level with the rest of the physical world might be simply wrong. As a loose analogy, if the physical reality were the water of a sea, our brain would not be like a submarine guided by sonar, rather it would be more like a soaked sponge.

So far, we have never made mention to how consciousness rises, or to which are the specific physical ‘mechanisms’ at the basis of consciousness, but this was not the topic of our paper, other than being a very complicated and long-standing issue (irresolvable at the moment, I think). Rather, I have proposed a possible physical characterization of consciousness and mind: the uniqueness of conscious experience and consciousness might depend also on the whole surrounding physical world, in the sense of the organization of its physical matter and energy. But, if this were the case, it would be very strange if the rise itself of consciousness and mind did not depend on the extensive interaction of the brain with whole surrounding physical world.

By the way, if we do not accept the picture of the extensive physical interaction between brain and the whole surrounding physical world, but want to maintain the dependence of mind and consciousness on physical reality, and thus that mind and consciousness completely arise from the physical organization of the brain alone, all this would inevitably result in another type of ‘extensive interdependence’, a sort of non-locality: as a matter of fact, after the brain duplication, the very same mind (the very same person) could be in two very distant places at the same time.

I believe that talking about the necessary physical mechanisms responsible for the suggested extensive interaction, and actually talking about the detailed physical dynamics and evolution of such interaction, is premature at the moment. It might be said (and it was actually said already) that the human brain may be assimilated to a complex dynamical system, extremely sensitive to all the surrounding physical conditions (in such case we should mention the so-called deterministic chaos; see for example Newman ), or that quantum mechanics may be deeply involved (many people have proposed various quantum mechanisms to explain consciousness, so to compile a complete list of references on this topic is hopeless; see for example Bohm  and references therein), or, most probably, that a new kind of physics is needed. However, although such aspect of the problem is obviously essential in the study of the origin of consciousness, it is secondary in the present context.

## Acknowledgments

I wish to thank Alessandro Silvestrini for having brought to my attention the ‘brain duplication’ argument, during a conversation around a coffee table. By the way, the name of the Cafè was ‘Miró’, thus for the sake of joke, I suggest to call the argument of this paper ‘the Miró hypothesis’. I am also grateful to William A. Adams for insightful comments on this paper. I wish to thank Herbert F. Muller for having posted a previous version of this paper on the Karl Jaspers Forum Website.


text:
# On the Speed of Gravity and Relativistic 𝑣/𝑐 Corrections to the Shapiro Time Delay

and

## Abstract

Recent papers by Samuel  declared that the linearized post-Newtonian $`v/c`$ effects are too small to have been measured in the recent experiment involving Jupiter and quasar J0842+1845  that was used to measure the ultimate speed of gravity defined as a fundamental constant entering in front of each time derivative of the metric tensor in the Einstein gravity field equations. We describe our Lorentz-invariant formulation of the Jovian deflection experiment and confirm that $`v/c`$ effects are do observed, as contrasted to the erroneous claim by Samuel, and that they vanish if and only if the speed of gravity is infinite.

Gravitational light-ray deflection experiments testing time-dependent components of the gravitational field are important for deeper understanding of the causal (null-cone) structure of the space-time manifold and the Einstein theory of general relativity. General relativity predicts that a light particle (photon) is deflected by the gravitational field of a moving body from its retarded, with respect to observer, position . This is a direct physical consequence of the relativistic nature of the gravitational field which is expressed mathematically in terms of the retarded Lienard-Wiechert solution of the linearized Einstein equations . We proposed relativistic VLBI experiment  to measure this effect of retardation of gravity by the field of moving Jupiter via observation of light bending from a quasar and successfully completed this experiment in September 2002. General relativistic experimental model of VLBI data processing algorithm is described in  and the results of the main experiment of September 2002 are published in the Astrophysical Journal .

Recently Samuel  made an attempt to review the results of the experiment by making use of a linearized Lorentz transformation of the static Shapiro time delay. This approach is physically meaningful but correct conceptual analysis and interpretation of the experiment requires matching of the first, $`v/c`$, and second $`v^2/c^2`$ order effects in the relativistic theory of the Shapiro time delay which has neither been conceived nor completed in . This led Samuel  to the confusion of the Lorentz transformation of the gravitational field and that of light used for observations. The same conceptual flaw is contained in papers by Asada  and Will  who were not able to separate a general-relativistic effect of the aberration of gravity from the special-relativistic effect of the aberration of light . Moreover, Samuel  has arrived to erroneous conclusion that the post-Newtonian effect of order of $`v/c`$ beyond the static part of the Shapiro effect can not be observed at all. In the present Letter we outline the basic equations of the exact Lorentz-invariant theory of the gravitational time delay of light and show mistakes in Samuel’s interpretation  of the gravitational light-ray deflection experiments.

For Jupiter, the perturbation of the metric tensor, $`h^{\mu \nu }`$, in harmonic coordinates $`x^\alpha =(ct,𝒙)`$ is obtained after solving the linearized Einstein equations in terms of the retarded Lienard-Wiechert gravitational potentials . It yields

$$h^{\mu \nu }=\frac{2GM_J}{c^4}\frac{2u^\mu u^\nu +\eta ^{\mu \nu }}{r_\alpha u^\alpha }.$$
(1)
Here 111Greek indices run from 0 to 3. Bold letters denote spatial vectors. Repeated indices mean the Einstein summation rule. Euclidean dot and cross products of two vectors are denoted as $`𝒂𝒃`$ and $`𝒂\times 𝒃`$ respectively. Partial derivative with respect to coordinate $`x^\alpha `$ is denoted as $`_\alpha /x^\alpha `$. $`\eta ^{\mu \nu }`$ is the Minkowski metric, $`M_J`$ is the mass of Jupiter, $`u^\alpha =\gamma (s)(c,𝒗_J(s))`$ is its four-velocity, with $`𝒗_J(s)=d𝒙(s)/ds`$ and $`\gamma (s)=(1v_J^2(s)/c^2)^{1/2}`$. The distance $`r^\alpha =x^\alpha x_J^\alpha (s)`$, and Jupiter’s worldline $`x_J^\alpha (s)=(cs,𝒙_J(s))`$ are functions of the retarded time $`s`$, determined as a solution of null-cone equation for gravity

$$\eta _{\mu \nu }r^\mu r^\nu =0.$$
(2)
Substituting $`r^\mu (s)`$ and solving (2) in the causality region, $`t>s`$, show that the retarded time $`s`$ obeys a (non-linear) equation

$$s=t\frac{1}{c}|𝒙𝒙_J(s)|,$$
(3)
which describes the null-cone characteristics of Jupiter’s gravity field with $`c`$ being understood as the ultimate speed of gravity. Retarded time $`s`$ appears as an argument in coordinates of Jupiter and reflects the obvious fact that gravity field is causal and propagates with finite speed $`c`$ as predicted by Einstein’s theory of general relativity . Despite that this derivation is well-known and given in the advanced course on general relativity  some researchers  are still confused about the nature of the retarded-time equation (3) in the gravitational Lienard-Wiechert potentials (1) and misinterpret it as a classic Römer time delay of light propagating from a massive body (Jupiter) to the field point.

VLBI measures the phase $`\varphi `$ of the radio wave front coming from a quasar. Diffferential equation for the phase is determined by the first integral of the light-ray geodesics that is known as the eikonal equation

$$g^{\mu \nu }_\mu \varphi _\nu \varphi =0,$$
(4)
where the metric tensor $`g^{\mu \nu }=\eta ^{\mu \nu }h^{\mu \nu }`$. Provided that one neglects orbital acceleration of Jupiter, solution of equation (4) is given by

$$\varphi =\varphi _0+\frac{\nu }{c}\left[k_\alpha x^\alpha +\frac{2GM_J}{c^2}\left(k_\alpha u^\alpha \right)\mathrm{ln}\left(k_\alpha r^\alpha \right)\right],$$
(5)
where $`\varphi _0`$ is a constant, $`\nu `$ is the radio frequency, $`k^\alpha =(1,𝒌)`$ is a null vector of the radio wave, and the unit vector $`𝒌`$ is directed from the quasar to VLBI station. Notice that the argument of the logarithm in equation (5) depends on the retarded time $`s`$ defined by the gravity null-cone equation (3). It makes evident that the phase of electromagnetic wave passing through time-dependent gravitational field of moving Jupiter is affected by this field not instantaneously but with a delay being required for gravity to propagate from the retarded position $`𝒙_J(s)`$ of Jupiter to the position $`𝒙`$ of a radio photon. Measuring this ‘retardation-of-gravity’ effect in the light-ray deflection experiments like that we had conducted in September 2002 , allows one to corroborate that gravity propagates on a null cone with finite speed .

Equation (5) is Lorentz-invariant and was obtained by making use of the Lienard-Wiechert gravitational potentials (1). It can be also derived by solving the eikonal equation (4) in a co-moving frame of Jupiter with subsequent Lorentz transformation of both gravitational and electromagnetic field from this frame to the barycentric frame of the solar system . This point is subtle and was not properly understood by Samuel  and some other researchers  who erroneously believe that the linearized (that is, $`v/c`$) approximation of the Lorentz transformation affects exclusively the direction of propagation of radio wave used for observation of the gravitational light-ray deflection but leaves the gravitational field of Jupiter in the barycentric frame of the solar system static (see also ). The message of this Letter is that in the barycentric frame of the solar system the gravitational field is characterized in the neigborhood of Jupiter not only by the components of the metric tensor $`g_{\alpha \beta }`$ but its space and time derivatives (the Christoffel symbols and the curvature tensor) as well. The eikonal equation (4) is the first integral of the light geodesic equations . This integral can be present in the form given by equation (4) if and only if the speed of gravity and the speed of light are equal or, in other words, the time derivatives of the metric tensor in the Einstein equations are normalized to the same fundamental speed as the time derivatives of electromagentic field in the Maxwell equations. Lorentz transformation from Jupiter’s frame to the barycentric frame of the solar system induces non-zero values of the time derivatives of the metric tensor if the speed of gravity is not infinite 222Notice that the speed of physical light does not play any role in the Lorentz transformation of the gravitational field.. In general relativity the effect of these first time derivatives of the metric tensor cancel with the aberration of light when one integrates equation of light geodesics to obtain the eikonal equation (4) which leaves the retardation of gravity effect present in Jupiter’s coordinates $`𝒙_J(s)`$ in equation (5). Therefore, the retardation of gravity effect appears in the linearized order of $`v/c`$ in the post-Newtonian expansion of the retarded coordinate $`𝒙(s)`$ of Jupiter around the time of observation $`t`$  as we demonstrate in the next paragraph. It can be easily confused with the aberration of light if the effect of the first time derivatives of the metric tensor is not pursued in derivation of equation (5) for the electromagnetic phase. The reason for the ‘aberration-of-light’ misinterpretation of the Jovian light-ray deflection experiment given in papers  is that the Nordtvedt-Will PPN formalism  has a faulty assumption leading to the replacement of the speed-of-gravity parameter $`c_g`$ in front of the first time derivatives of the metric tensor with the physical speed of light . This makes the PPN gravitational force acting on test particles, be dependent upon the speed of electromagnetic waves. Thus, the PPN supposes that the gravitational force is tightly connected with Maxwell’s theory which is obvously a wrong physical postulate leading to a conceptual failure of the Nordtvedt-Will PPN formalism in case of time-dependent gravitational fields.

Let us now consider the post-Newtonian expansion of the eikonal (5). Toward this end we approximate Eq. (5) by introducing two angles $`\mathrm{\Theta }`$ and $`\theta `$ between vector $`𝒌`$ and unit vectors $`𝒑=𝑹/R`$ and $`𝒍=𝒓/r`$, where vectors $`𝑹=𝒙𝒙_J(t)`$ and $`𝒓=𝒙𝒙_J(s)`$ connects the point of observation $`𝒙𝒙(t)`$ with the present, $`𝒙_J(t)`$, and retarded, $`𝒙_J(s)`$, positions of Jupiter, respectively. By definition $`\mathrm{cos}\mathrm{\Theta }=𝒌𝒑`$ and $`\mathrm{cos}\theta =𝒌𝒍`$. The product $`k_\alpha r^\alpha =r(1\mathrm{cos}\theta )r\theta ^2/2`$ for small angles which is the case in the gravitational light-ray deflection experiments. Hence, equation (5) tells us that the phase variation caused by space-time difference $`\delta x^\alpha =(c\delta t,\delta 𝒙)`$ between two antennas is

$$\delta \varphi =k_\alpha \delta x^\alpha +\frac{4GM_J}{c^2}\frac{\delta \theta }{\theta },$$
(6)
where $`\delta \theta =𝒏\delta 𝒙/r`$ with $`𝒏=𝒍\times (𝒌\times 𝒍)`$ as the impact unit vector of the light ray with respect to the retarded position of Jupiter $`𝒙_J(s)`$. Currently undetectable terms of order $`v_J/c`$ have been neglected. The quantity $`\delta t`$ is the measurable time delay and $`\delta 𝒙=𝑩`$ is a baseline between two VLBI stations. Since VLBI stations measure the same wave front, $`\delta \varphi =0`$. Thus, Eq. (6) gives

$$\delta t=c^1𝒌\mathbf{}𝑩+\mathrm{\Delta },$$
(7)
where

$$\mathrm{\Delta }=\frac{4GM_J}{c^3r}\frac{𝒏\mathbf{}𝑩}{\theta }.$$
(8)
The post-Newtonian expansion of the retarded position of Jupiter, $`𝒙_J(s)`$, in Eq. (8) with respect to the ratio $`v_J/c`$ (velocity of Jupiter/speed of gravity) around time of photon’s arrival $`t`$ yields

$$\mathrm{\Delta }=\mathrm{\Delta }_S+\mathrm{\Delta }_R.$$
(9)
Here

$$\mathrm{\Delta }_S=\frac{4GM_J}{c^3}\frac{(𝑵\mathbf{}𝑩)}{\mathrm{\Theta }},$$
(10)
is the static (instantaneous-like) part of the Shapiro time delay, and the gravity retardation $`v_J/c`$-correction

$$\mathrm{\Delta }_R=\frac{4GM_J}{c^4R\mathrm{\Theta }^2}\left[2(𝑵𝒗_J)𝑵𝒌\times (𝒗_J\times 𝒌)\right]\mathbf{}𝑩,$$
(11)
where $`𝑵=𝒑\times (𝒌\times 𝒑)`$ is the impact unit vector of the light ray with respect to the present position of Jupiter $`𝒙_J(t)`$ taken at the time of observation from JPL ephemeris . Eq. (11) is the same as Eq. (4) from . It describes the post-Newtonian $`v_J/c`$ correction to $`\mathrm{\Delta }_S`$ and can be detected because of the amplifying factor $`1/\mathrm{\Theta }^2`$.

Samuel  incorrectly assumed that the experiment directly compared the radio position of the quasar with that of Jupiter, and that the direction of Jupiter was determined by light reflected from its surface. This assumption would correspond to direct measurement of the angle $`\theta `$ and hence no $`v_J/c`$ terms would be observed since they are not evident in Eq. (8). The experiment, however, monitored the position of the quasar as a function of the atomic time by the arrival of the quasar’s photons at the radio telescope, while the Jupiter’s position was not directly observed 333Jupiter has too large angular size in the sky compared to VLBI directional diagram. but determined separately via a precise JPL ephemeris 444JPL ephemeris of Jupiter is constructed from radio-tracking of spacecrafts orbiting Jupiter. Hence, the experiment  effectively compared the retarded position of Jupiter $`𝒙_J(s)`$ derived from the gravitational time delay with that obtained from the spacecraft’s radio-tracking. We have used the same idea to measure the aberration of gravity in another VLBI experiment in October 2005 . , evaluated at the same atomic time as the arrival of a radio photon from the quasar (via standard transformations from barycentric time to atomic time ). Hence the actual angle used in VLBI correlator for measuring $`\mathrm{\Delta }`$ in equation (9) is $`\mathrm{\Theta }`$, not Samuel’s $`\theta _{obs}\theta `$. Thus, the $`v_J/c`$ correction $`\mathrm{\Delta }_R`$ was clearly separated in the VLBI correlator Calc/Solve program  from $`\mathrm{\Delta }_S`$ and measured with a precision of 20% . Additional critical analysis of the theoretical concepts involved to the interpretation of the Jovian light-ray deflection experiment can be found in . Their study can be useful for experts working on future projects like Gaia space mission  and Square Kilometer Array .

We thank S. Samuel, G. Schäfer, S. Carlip, C. M. Will for valuable comments and V. B. Braginsky and M. V. Sazhin for fruitful and constructive discussions.


text:
# Exchange-correlation energy of a multicomponent two-dimensional electron gas

## I Introduction

One of the most interesting discoveries in semiconductor heterostructures was the formation of a two-dimensional (2D) electron gasando1982 . For example, this has set the starting point for building nanostructures such as quantum dots or wires (for reviews see chakraborty1999 ; reimann2002 ; wiel2003 ). Many of the simplified models describing semiconductor nanostructures are based on the assumption of a homogeneous electron gas. The effects of the detailed band structure are included only by using an effective mass for the electrons and by screening the electron-electron interaction with a static dielectric constant. For the conduction electrons of a direct semiconductor, like GaAs for example, this approach describes fairly well the underlying electronic structurechakraborty1999 ; reimann2002 .

In the following, we call such a normal 2D electron gas as a two-component gas, the components being the spin-up and spin-down electrons. Similarly, a polarised electron gas is called a one-component Fermi gas. (In the literature, the latter is sometimes referred to as “spinless fermions”).

Layered semiconductors, however, provide several structures where the simple picture of an ideal (one- or two-component) two-dimensional electron gas is likely to fail. In elemental semiconductors, silicon and germanium, the conduction band minimum consists of four equivalent minima at nonzero values of the $`k`$-vector and, moreover, these minima have non-isotropic effective masses. The resulting conduction electron gas has an “internal” degeneracy of eight (two from spin and four from the four minima), and is a possible example of an eight-component Fermi gas. The correlations in such an electron gas are much more complicated and have been studied in connection with the 3D electron-hole plasma in semiconductorsvashista1974 .

The valence band maximum in most semiconductors is degenerate, consisting of two bands corresponding to heavy and light holes. In 3D structures the hole gas will then consist of two kinds of holes. In 2D structures there will be an energy separation between the band minima due to the fact that heavy holes can be more easily localised than the light hole. Nevertheless, if the mass difference is not too large, both kinds of holes can coexist and the correlations between them have to be considered. The hole gas of most semiconductors is one example of a four-component Fermi gas. Electron layers in semiconductor heterostructures provide yet another system of multicomponent electron (or hole) gas, intensively studied theoreticallyneilson1993 ; sarma1994 ; alatalo1994 ; conti1996 . Recently, electron addition spectra of vertical quantum dots based on layered structures have also been measuredausting1998 .

Density functional theory (DFT) in the local (spin) density approximation (LSDA) has shown to be a useful starting point for studies of the ground state properties of electrons in quantum dotsreimann2002 . The 2D conduction electron gas of the semiconductor is then treated as an interacting electron gas. In the LSDA, the exchange-correlation energy is approximated by that of the ideal 2D electron gas, the effective mass and the static dielectric constant appearing only as scaling factors (for energy and distance). Applying the density functional method to the above mentioned multicomponent structures calls for a parametrisation of the exchange-correlation energy ($`ϵ_{xc}`$) for the multicomponent system in question. Tanatar and Ceperleytanatar1989  computed the one and two-component (polarised and non-polarised 2D electron gas) using the quantum Monte Carlo technique. More recently Attaccalite et al.attaccalite2002  presented new Monte Carlo results for these systems. They also derived an interpolated expression for the exchange-correlation energy functional to be used in DFT calculations. The most important difference to previous parametrisations like for example the one by Tanatar and Ceperlytanatar1989  is that here the Monte Carlo-calculations are performed also at intermediate polarisations. The usual approximation to treat the polarisation dependence of the correlation energy on the same footing as the (exactly known) exchange partreimann2002  can thus be abandoned.

For a multicomponent system the situation becomes more complicated and the transferability of $`ϵ_{xc}`$ from one system to another is not so obvious. The purpose of this paper is to estimate the exchange-correlation energy of a general multicomponent 2D fermi gas. We will first discuss the general high and low density limits and the effect of mass difference between the components. We will then show that an analytic continuation of the parametrisation by Attaccalite et al.attaccalite2002  fits well the existing many-body calculations for multicomponent systems. As reference data we use the calculations of Conti and Senatoreconti1996  for a four-component 2D electron gas and the results of Apaja et al.apaja1997  for a 2D charged Bose gas which can be viewed as an upper bound for the energy of a Fermi gas with an infinite number of internal degrees of freedom.

Finally, we will study the phase diagram suggested for this energy functional and show that already a small mass difference will decrease the concentration of the lighter-mass components. When the density gets low enough, eventually all systems become one-component, as predicted by the Attaccalite et al.attaccalite2002  parametrisation for the normal 2D electron gas.

## II Exchange-Correlation energy: General discussion

We consider a homogenous gas consisting of $`\mathrm{\Lambda }`$ different kinds of fermions, or components, as we call them. For instance, the different spin states of the homogenous gas are treated as different components. For a normal electron gas $`\mathrm{\Lambda }=2`$ and for a fully polarised gas $`\mathrm{\Lambda }=1`$. We assume that all the components have the same charge, but they may have different masses. The total density $`n`$ of the gas is

$$n=\underset{i=1}{\overset{\mathrm{\Lambda }}{}}n_i=n\underset{i=1}{\overset{\mathrm{\Lambda }}{}}\nu _i,$$
(1)
where $`n_i`$ is the number density and $`\nu _i=n_i/n`$ is a dimensionless concentration of the component $`i`$. Note that the concentrations are normalised as $`_i^\mathrm{\Lambda }\nu _i=1`$. The density parameter $`r_s=1/\sqrt{\pi n}`$ is always related to the total density $`n`$ of the gas. The kinetic energy per particle of the noninteracting gas becomes

$$ϵ_k(r_s,\{\nu _i\},\{m_i\})=\frac{\mathrm{}^2}{r_s^2}\underset{i=1}{\overset{\mathrm{\Lambda }}{}}\frac{\nu _i^2}{m_i},$$
(2)
where $`m_i`$ is the mass of the component $`i`$. The total exchange energy is the sum of exchange energies of the different components. The total exchange energy per particle is then

$$ϵ_x(r_s,\{\nu _i\})=\frac{e^2}{4\pi \epsilon _0ϵ}\frac{8}{3\pi r_s}\underset{i=1}{\overset{\mathrm{\Lambda }}{}}\nu _i^{3/2}.$$
(3)
Note that while the kinetic energy depends on the masses of the components, the exchange energy is independent of the masses. The correlation energy per particle is defined as the difference between the exact energy and the Hartree-Fock energy, $`ϵ_c=ϵ_{tot}(ϵ_k+ϵ_x)`$. In true semiconductor systems the different components can represent, for example, different band minima. The exchange energy has then also components coming from the exchange between electrons belonging to these different minima. Generally such exchange interaction is found to be smallbeliaev1998 ; hada2003 . In the spirit of the effective mass approximation we will completely neglect such effects here.

In order to study the multicomponent correlation energy we will first consider the case where all the components have the same mass. In this case the exchange correlation energy can be written as a function of the density parameter $`r_s`$ and the concentrations $`\nu _i`$. The symmetry requires that the energy functional is symmetric with respect of the interchange of the different concentrations $`\nu _i`$, i.e., it does not depend on the indexing. We define numbers

$$Z_\gamma =\underset{i=1}{\overset{\mathrm{\Lambda }}{}}\nu _i^\gamma ,$$
(4)
which by construction have the desired symmetry. We notice that $`Z_1=1`$, $`ϵ_xZ_{3/2}`$ and the kinetic energy $`ϵ_kZ_2`$ (for equal masses). In general any function $`g(\{\nu _i\})`$ with the desired symmetry property can be written as a function of the $`Z_\gamma `$-numbers, for example, as a series expansion

$$g(\{\nu _i\})=f(\{Z_{\gamma _i}\})=a_1Z_{\gamma _1}+a_2Z_{\gamma _2}+a_{12}Z_{\gamma _1}Z_{\gamma _2}+\mathrm{},$$
(5)
where $`a_i`$ are constants and $`\gamma _i`$ can be noninteger. One way to approach the correlation energy is to use the suggestion of Wignerwigner1934  who interpolated the exchange-correlation energy of an electron gas between the two limiting cases, the pure exchange in the high density limit and the energy of the Wigner crystal in the low density limit:

$`\begin{array}{c}ϵ_{xc}(r_s,\{\nu _i\},\{m_i\})=ϵ_{WC}(r_s)\hfill \\ \\ +g(r_s,\{\nu _i\},\{m_i\})[ϵ_x(r_s,\{\nu _i\})ϵ_{WC}(r_s)]\hfill \end{array}`$ (9)
where $`ϵ_{WC}`$ is the classical Madelung energy of the Wigner crystal, which does not depend on the masses. The advantage of this approach is that the mass dependence appears only in the interpolation function $`g`$. Furthermore, if the masses are equal, the interpolation function $`g`$ can immediately be written as in Eq. (5).

In order to estimate the interpolating function $`g`$ one would need results of accurate many-body calculations for systems with different numbers of components, and for several values of $`r_s`$ and the concentrations $`\nu _i`$. We are not in the position to perform such computations at this point but use existing data for special cases. For the two-component system Attaccalite et al.attaccalite2002  have performed Monte Carlo calculations. They fitted an accurate interpolation formula to the results which is valid for any $`r_s`$ and polarisation.

For the four-component system, Conti and Senatoreconti1996  have presented Monte Carlo results for several values of $`r_s`$ in the case where all the concentrations $`\nu _i`$ are equal. For systems with more components there exist no accurate data. However, the limit of infinite components can be approximated by the results obtained for charged bosonsapaja1997 . If all the concentrations are the same ($`\nu _i=1/\mathrm{\Lambda }`$) in the limit of infinitely many components, all the particles are at different internal state and the Pauli exclusion principle can not prevent putting all the particles in the same orbital state. The symmetric boson state is then a legitime state for such a fermion system and provides an upper bound for the energy. However, the existing results beyond the one and two-component systems are still quite limited for an accurate interpolation of the function $`g`$.

### II.1 Mass dependence

As discussed above, the exchange energy is mass-independent and even in the multi-component electron gas only the kinetic and the correlation energies will depend on the masses. If all the components have the same mass $`M`$, the kinetic energy of Eq. (2) can be written as $`ϵ_k=\mathrm{}^2Z_2/(r_s^2M)`$. In the general case, the kinetic energy can still be written in the same form by defining an average mass $`M`$ as

$$\frac{1}{M}=\frac{1}{Z_2}\underset{i}{\overset{\mathrm{\Lambda }}{}}\frac{\nu _i^2}{m_i}.$$
(10)
We notice a scaling: By deviding the the kinetic and exchange energies by $`M`$, the resulting ratios will only depend on the product $`Mr_s`$ ($`ϵ_k/M=\mathrm{}^2Z_2/(Mr_s)^2`$ and $`ϵ_x/M=8e^2Z_{3/2}/(4\pi \epsilon _0ϵ3\pi Mr_s)`$). Moreover, considering the atomic units in the case that all masses are the same, we notice that (in this case) also the total energy, and thus the correlation energy has the same scaling. It is then a reasonable first approximation to take the mass as an average mass and use it as a scaling factor

$$ϵ_{xc}(r_s,\{\nu _i\},\{m_i\})=\frac{M}{m_e}ϵ_{xc}(Mr_s,\{\nu _i\},\{m_i=m_e\}),$$
(11)
where $`m_e`$ is the bare mass of the electron (or any suitably chosen effective mass). (Note that within this approximation, the interpolation function $`g`$ in Eq. (9) would be written as $`g(Mr_s,\{Z_\gamma \})`$).

## III Extension of the two-component function

For estimating an interpolation formula for the multicomponent correlation energy, we start by considering the parametrised form of the two-component 2D electron gas by Attaccalite *et al.* attaccalite2002 , in order to see if their function can be written in terms of the $`r_s`$ and the $`Z_\gamma `$-numbers, as discussed above. Attaccalite *et al.* use the common description where the energy is written as a function of $`r_s`$ and the polarisation $`\zeta =(n_1n_2)/n=\nu _1\nu _2`$ of the electron gas. They chose a correlation energy function $`e_c(r_s,\zeta )=ϵ_{xc}(r_s,\zeta )ϵ_x(r_s,\zeta )`$ which diminishes the contributions of $`ϵ_x`$ beyond fourth order in $`\zeta `$ as $`r_s`$ increases. The parametrisation satisfies known high- and low-density limits: For high densities i.e. small $`r_s`$ the exchange-correlation energy obeys the perturbation theory result $`ϵ_{xc}=ϵ_x+a_0(\zeta )+b_0(\zeta )r_s\mathrm{ln}r_s`$. For low densities the Wigner crystal limit is recovered, $`ϵ_{xc}c_1/r_s+c_2/r_s^{3/2}`$, where the constants $`c_1`$ and $`c_2`$ are independent of the number of components.

Using the concentrations $`\nu _i`$, the square of the polarisation $`\zeta `$ of the two-component gas becomes $`\zeta ^2=2(\nu _1^2+\nu _2^2)1=2Z_21`$. Now the exchange-correlation energy introduced by Attacalite et al.attaccalite2002  can be written as (using their notations except for the above replacement for $`\zeta `$)

$`\begin{array}{c}ϵ_{xc}(r_s,\{\nu _i\})=e^{\beta r_s}[ϵ_xϵ_x^{(6)}]\hfill \\ \\ +ϵ_x^{(6)}+\alpha _0(r_s)+\alpha _1(r_s)(2Z_21)+\alpha _2(r_s)(2Z_21)^2,\hfill \end{array}`$ (15)
where $`ϵ_x^{(6)}=(1+\frac{3}{8}(2Z_21)+\frac{3}{128}(2Z_21)^2)ϵ_x(r_s,\zeta =0)`$. The functions $`\alpha _i(r_s)`$ are parametrised by Attacalite et al.attaccalite2002  and are independent of the concentrations. This two-component exchange-correlation function of Attaccalite et al is thus already of the desired form: It depends only on $`r_s`$, $`Z_2`$ and $`Z_{3/2}`$ (through $`ϵ_x`$) and can be directly extended to any number of components.

For a given value of $`r_s`$ Attaccalite et alattaccalite2002  fitted their function to several values of polarisations, $`\zeta [0,1]`$, which correspond to $`Z_2[0.5,1]`$ and $`Z_{3/2}[\sqrt{2},1]`$. An analytic continuation to larger number of components means an extrapolation of the $`ϵ_{xc}`$ to the region $`Z_2[0,1]`$ and $`Z_{3/2}[0,1]`$ for a system with infinitely many components, but only to $`Z_2[0.25,1]`$ and $`Z_{3/2}[0.5,1]`$ for a four-component system. Moreover, since $`Z_{3/2}`$ only appears in the exchange energy, Eq. (3), which is exact for any number of components, errors in the extrapolation only result from the extension of the $`Z_2`$ space.

The multicomponent gas must obey the low density Wigner crystal limit and the high density exchange limit; both of these requirements are satisfied by Eq. (15). By taking a “paramagnetic” gas, i.e. by setting $`\nu _i=\mathrm{\Lambda }^1`$ for all the components and fixing $`r_s`$, it can be seen that the exchange-correlation energy above is also a monotonous function of the number of components.

Figure 1 shows the exchange-correlation energy derived from Eq. (15) for selected values of $`r_s`$ as a function of $`Z_2`$. For the interval $`Z_2[.5,1]`$ the lines correspond to the normal two-component electron gas with partial polarisation. For the interval $`Z_2[0,.5]`$ the lines are derived for the case where all the concentrations are equal, $`\nu _i=1/\mathrm{\Lambda }`$ and consequently $`Z_2=1/\mathrm{\Lambda }`$. The points are results of different many-body calculations. Note that the lines from $`Z_2=0.5`$ to $`Z_2=1`$ agree exactly to the results of Attaccalite et alattaccalite2002 . We can see that the extension to a four component system ($`Z_2=.25`$) is fairly accurate and the extension even to the boson case ($`Z_2=0`$) is still reasonable. In fact, the agreement of the Attaccalite et al.-function with the existing many-body results for four- and infinite-component systems is so good that there is no need to make a better interpolation formula for the multicomponent system, before more total energy calculations are available in the region $`Z_2<.5`$. It would be important to have results also for cases where the concentrations of all the components are not equal.

## IV Phase diagram of the multicomponent 2D electron gas

It was recently suggested by Attacalite et al.attaccalite2002  that the 2D electron gas shows a transition to a polarised gas at $`r_s25`$attaccalite2002 . We will first study the stability of the multicomponent gas assuming all the masses to be equal. The extended exchange-correlation energy functional suggests that all multicomponent cases have a transition to a one-component phase when $`r_s`$ increases. Table I gives the estimated transition points. For any number of components, the transition happens nearly at the same point and directly from the multicomponent to a one-component phase. It is caused by the slightly different $`r_s`$-dependence of the one-component total energy as compared to the others. For $`\mathrm{\Lambda }2`$ the total energy is a monotonously decreasing function of $`\mathrm{\Lambda }`$ (for any $`r_s`$).

Next we will study the effect of the mass difference, using the approximation of Eq. (11). We assume a four-component system and fix the masses pairwise equal so that two components are heavier than the two others, say $`m_1=m_2=m_l`$ and $`m_3=m_4=m_h`$. Note that due to the spin degeneracy of electrons (or holes) there will always be two components with the same mass, and the number of components will be even, except of the special case one. The phase diagram of the four-component electron gas obtained using Eq. (15) is shown in figure 2. We show the maximum of the concentrations ($`\nu _4`$) as a function of $`r_s`$ and the mass ratio $`m=m_h/m_l`$. When the masses are equal, also the concentrations remain equal ($`\nu _i=0.25`$) until $`r_s26.2`$ at which point the system becomes a one-component gas. When the mass ratio increases, the concentrations $`\nu _3`$ and $`\nu _4`$ increase due to the reduction of the kinetic energy, and eventually $`\nu _3=\nu _4=0.5`$ and the system has become two-component. A further increase of the mass ratio eventually changes the system to a one-component gas. We notice that as the mass is increased the polarisation occurs at smaller $`r_s`$. The reason is the scaling of the total energy and the dimensions with the mass.

In real semiconductor heterostructures the mass difference can arise from different constituents of the two layers in a double layer system or, in the case of holes, simply from the mass difference of the heavy and light hole. In both cases, however, there can be a constant potential energy difference from one component to another. This energy difference arises from the localisation of the particle in the direction perpendicular to the 2D-layer. For example, if we consider holes in a 1D harmonic potential, the ground state energy (the lowest perpendicular mode) has the energy $`ϵ_0(m)=\mathrm{}\omega _0/(2\sqrt{m/m_e})`$, where $`\omega _0`$ is the confinement strength for $`m_e`$. This energy is smaller for the heavy hole and, consequently, further favours the transition to a two-component system where only heavy holes exist. In order to have a 2D electron gas the Fermi energy has to be clearly below the first excited perpendicular state, $`ϵ_1(m)=3\mathrm{}\omega _0/(2\sqrt{m/m_e})`$ in a harmonic well. In order for the heavy (mass $`m_h`$) and light holes (mass $`m_l`$) to coexist, the ratio $`(ϵ_0(m_l)ϵ_0(m_h))/(ϵ_1(m_h)ϵ_0(m_h))`$ has to be clearly smaller than one. For example in the cases of Si, Ge and GaAs ($`m_h/m_l>5`$) this ratio is larger than 0.6. Consequently, a 2D layer with both heavy and light holes is not possible. However, there are semiconductors, like AlAs, where $`m_h/m_l2`$ and the above energy ratio is only about 0.2, making the four-component gas possible.

## V Conclusions

We have studied the possibility to formulate a general multicomponent 2D exchange-correlation energy to be used in density functional calculations. After general considerations of the concentration and mass dependence of the correlation energy, we have demonstrated that the recent parametrisation by Attacalite et al.attaccalite2002  of the two-component exchange-correlation energy has functional properties which allow a direct extension to multicomponent systems. Furthermore, although the function is fitted only to the two-component electron gas, it gives fairly accurate results for four-component and infinite-component systems.

We suggest that the mass dependence of the correlation energy can be estimated with a properly chosen average mass. The application for multicomponent 2D gas shows an interesting phase diagram. For equal masses all systems transform from a multicomponent directly to a one-component gas when $`r_s`$ increases to about $`25\mathrm{}27`$ depending on the number of components. Even a small mass difference decreases the concentration of the lighter component and it transforms first to a purely two-component and then to a one-component system when $`r_s`$ increases.

A more accurate interpolation formula requires extensive many-body calculations for multi-component gases. We hope that these considerations encourage such work.


text:
# Quantum phases of atomic boson-fermion mixtures in optical lattices

## Abstract

The zero-temperature phase diagram of a binary mixture of bosonic and fermionic atoms in an one-dimensional optical lattice is studied in the framework of the Bose-Fermi-Hubbard model. By exact numerical solution of the associated eigenvalue problems, ground state observables and the response to an external phase twist are evaluated. The stiffnesses under phase variations provide measures for the boson superfluid fraction and the fermionic Drude weight. Several distinct quantum phases are identified as function of the strength of the repulsive boson-boson and the boson-fermion interaction. Besides the bosonic Mott-insulator phase, two other insulating phases are found, where both the bosonic superfluid fraction and the fermionic Drude weight vanish simultaneously. One of these double-insulator phases exhibits a crystalline diagonal long-range order, while the other is characterized by spatial separation of the two species.

Following the seminal experiments on the superfluid to Mott-insulator transition in ultracold atomic Bose gases in optical lattice potentials GrMa02a ; GrMa02b , a new frontier in the field of degenerate, strongly correlated quantum gases has emerged. It has recently become possible to prepare degenerate mixtures of a bosonic and a fermionic atomic species in an optical lattice MoFe03 . These degenerate gases will facilitate the experimental study of quantum phase transitions in systems of mixed quantum statistics, which are extremely hard to access in the solid-state context. The superior experimental control available in ultracold atomic gases experiments promises detailed insights into the rich physics of quantum phase transitions. In this paper we explore the phase diagram of a binary boson-fermion mixture in an optical lattice and identify experimental signatures of the different phases.

A theoretical framework for degenerate boson-fermion mixtures in optical lattices can be constructed by a straight-forward generalization of the single-band Bose-Hubbard model, used successfully to describe the superfluid to Mott-insulator phase transition in purely bosonic systems JaBr98 . A complete single-particle basis within the lowest Bloch band is given by the Wannier functions for the vibrational ground state of each lattice well. We define creation operators $`\widehat{\mathrm{a}}_i^{}`$ ($`\widehat{\mathrm{c}}_i^{}`$) which create a boson (fermion) in the localized Wannier state at site $`i=1,\mathrm{},I`$. Using this basis we can translate the many-body Hamiltonian into a second quantized form. This leads to the Bose-Fermi-Hubbard (BFH) Hamiltonian for a one-dimensional translationally invariant lattice AlIl03 ; BuBl03 ; LeSa03

$$\begin{array}{cc}\hfill \widehat{\mathrm{H}}=& J_\text{B}\underset{i=1}{\overset{I}{}}(\widehat{\mathrm{a}}_{i+1}^{}\widehat{\mathrm{a}}_i+\text{h.a.})J_\text{F}\underset{i=1}{\overset{I}{}}(\widehat{\mathrm{c}}_{i+1}^{}\widehat{\mathrm{c}}_i+\text{h.a.})\hfill \\ & +\frac{V_{\text{BB}}}{2}\underset{i=1}{\overset{I}{}}\widehat{\mathrm{n}}_i^\text{B}(\widehat{\mathrm{n}}_i^\text{B}1)+V_{\text{BF}}\underset{i=1}{\overset{I}{}}\widehat{\mathrm{n}}_i^\text{B}\widehat{\mathrm{n}}_i^\text{F},\hfill \end{array}$$
(1)
where $`\widehat{\mathrm{n}}_i^\text{B}=\widehat{\mathrm{a}}_i^{}\widehat{\mathrm{a}}_i`$ and $`\widehat{\mathrm{n}}_i^\text{F}=\widehat{\mathrm{c}}_i^{}\widehat{\mathrm{c}}_i`$ are the occupation number operators for bosons and fermions, respectively. The first line describes the tunneling between adjacent lattice sites. The tunneling amplitudes $`J_\text{B}`$ and $`J_\text{F}`$ are connected to off-diagonal matrix elements of the kinetic energy and the lattice potential AlIl03 . The second line of Eq. (1) contains the on-site boson-boson and boson-fermion interactions with interaction strengths $`V_{\text{BB}}`$ and $`V_{\text{BF}}`$, respectively. Within the single-band approximation a fermion-fermion on-site interaction does not arise as the Pauli principle prevents two identical fermions from occupying the same site.

We exactly solve the eigenvalue problem of the BFH Hamiltonian within a complete basis of Fock states of the form $`|\{n_1^\text{B},\mathrm{},n_I^\text{B}\}_\alpha |\{n_1^\text{F},\mathrm{},n_I^\text{F}\}_\beta `$. By $`\{n_1^\text{B},\mathrm{},n_I^\text{B}\}_\alpha `$ we denote a set of boson occupation numbers for the individual sites. The index $`\alpha =1,\mathrm{},D_\text{B}`$ labels all possible compositions of occupation numbers which yield $`_{i=1}^In_i^\text{B}=N_\text{B}`$. For the set of fermion occupation numbers $`\{n_1^\text{F},\mathrm{},n_I^\text{F}\}_\beta `$ only those compositions with $`0`$ or $`1`$ fermions at a given site are allowed by the Pauli principle. The dimensions of the boson and fermion space are given by $`D_\text{B}=(N_\text{B}+I1)!/(N_\text{B}!(I1)!)`$ and $`D_\text{F}=I!/(N_\text{F}!(IN_\text{F})!)`$, respectively. For a lattice with $`I=8`$ sites, $`N_\text{B}=8`$ bosons, and $`N_\text{B}=4`$ fermions the basis dimensions are $`D_\text{B}=6435`$ and $`D_\text{F}=70`$. Clearly, the dimension of the fermion space is dramatically reduced due to the Pauli principle.

Within the basis of Fock-states we construct the matrix representation of the BFH Hamiltonian assuming periodic boundary conditions. The lowest eigenvalues and eigenvectors of the Hamilton matrix are computed using a refined Lanczos algorithm. The eigenvectors provide the coefficients for the representation of the eigenstates $`|\mathrm{\Psi }_\nu `$ in the Fock-state basis

$$|\mathrm{\Psi }_\nu =\underset{\alpha =1}{\overset{D_\text{B}}{}}\underset{\beta =1}{\overset{D_\text{F}}{}}C_{\alpha \beta }^{(\nu )}|\{n_1^\text{B},\mathrm{},n_I^\text{B}\}_\alpha |\{n_1^\text{F},\mathrm{},n_I^\text{F}\}_\beta .$$
(2)
From this representation of the ground state we can extract various observables such as mean occupation numbers and number fluctuations for both species. It also provides access to more elaborate quantities such as the one- or two-body density matrix and the static structure factor, which we will discuss below.

The dynamical response of the system to external perturbations is a crucial quantity in distinguishing between different quantum phases, for example, between the bosonic Mott-insulator and the superfluid phase. By imposing a linear variation of the phase onto the many-body wavefunction one can probe the mobility of the atoms in the lattice. On the single particle level, a spatial variation of the phase of the wavefunction is associated with a velocity field $`\stackrel{}{v}(\stackrel{}{x})=\frac{\mathrm{}}{m}\stackrel{}{}\theta (\stackrel{}{x})`$, and a linear phase variation thus corresponds to a constant velocity across the lattice. Roughly speaking, those particles free to move will respond to the phase twist and exhibit a homogeneous flow. This flow is then associated with an additional kinetic energy and hence an increase of the total energy. This provides a measure for the density of mobile particles in the lattice RoBu03c ; RoBu03a .

This phase variation can be realized by imposing twisted boundary conditions for the many-body wavefunction. In the case of the BFH model it is, however, more convenient to map the phase twist by means of a unitary transformation onto the Hamiltonian RoBu03c . This leads to a “twisted” Hamiltonian which contains Peierls phase factors in the hopping terms:

$$\begin{array}{cc}\hfill \widehat{\mathrm{H}}_{\mathrm{\Theta }_\text{B},\mathrm{\Theta }_\text{F}}=& J_\text{B}\underset{i=1}{\overset{I}{}}(\mathrm{e}^{\mathrm{i}\mathrm{\Theta }_\text{B}/I}\widehat{\mathrm{a}}_{i+1}^{}\widehat{\mathrm{a}}_i+\text{h.a.})\hfill \\ & J_\text{F}\underset{i=1}{\overset{I}{}}(\mathrm{e}^{\mathrm{i}\mathrm{\Theta }_\text{F}/I}\widehat{\mathrm{c}}_{i+1}^{}\widehat{\mathrm{c}}_i+\text{h.a.})\hfill \\ & +\frac{V_{\text{BB}}}{2}\underset{i=1}{\overset{I}{}}\widehat{\mathrm{n}}_i^\text{B}(\widehat{\mathrm{n}}_i^\text{B}1)+V_{\text{BF}}\underset{i=1}{\overset{I}{}}\widehat{\mathrm{n}}_i^\text{B}\widehat{\mathrm{n}}_i^\text{F}.\hfill \end{array}$$
(3)
In order to probe the response of the bosonic and the fermionic species separately, we have introduced different total twist angles $`\mathrm{\Theta }_\text{B}`$ and $`\mathrm{\Theta }_\text{F}`$, respectively.

For an isolated boson twist ($`\mathrm{\Theta }_\text{B}>0,\mathrm{\Theta }_\text{F}=0`$) the difference in the ground state energies of the boson-twisted Hamiltonian $`\widehat{\mathrm{H}}_{\mathrm{\Theta }_\text{B},0}`$ and the initial Hamiltonian $`\widehat{\mathrm{H}}_{0,0}`$ provides direct information on the superfluid density of the bosonic component. This is because the superfluid component responds to the imposed phase gradient producing flow and the associated energy change. The bosonic phase stiffness

$$f_\text{s}^\text{B}=\frac{I^2}{N_\text{B}J_\text{B}}\frac{E_{\mathrm{\Theta }_\text{B},0}E_{0,0}}{\mathrm{\Theta }_\text{B}^2},\mathrm{\Theta }_\text{B}\pi $$
(4)
is, therefore, a measure for the boson superfluid fraction RoBu03a ; RoBu03c ; FiBa73 ; Krau91 . We note that neither the stability of the superflow nor the reduction of the flow by the lattice itself are taken into account by this quantity. The superfluid weight $`f_\text{s}^\text{B}`$ vanishes for insulating phases such as the Mott-insulator phase and goes to $`1`$ for a perfect superfluid.

The energy change resulting from an isolated fermion twist ($`\mathrm{\Theta }_\text{B}=0,\mathrm{\Theta }_\text{F}>0`$) provides a measure for the mobility of the fermionic atoms. The fermionic phase stiffness

$$f_\text{d}^\text{F}=\frac{I^2}{N_\text{F}J_\text{F}}\frac{E_{0,\mathrm{\Theta }_\text{F}}E_{0,0}}{\mathrm{\Theta }_\text{F}^2},\mathrm{\Theta }_\text{F}\pi $$
(5)
is equivalent to the well-known Drude weight for charged fermions in a lattice FyMa91 ; ScWh93 ; ShSu90 . The vanishing of the Drude weight is an indicator of an insulating state, whereas a finite value indicates non-vanishing conductivity.

These two phase stiffnesses are an important tool to characterize the different quantum phases of the boson-fermion mixture in a lattice. In the following we map out the phase diagrams for different boson and fermion numbers, $`N_\text{B}`$ and $`N_\text{F}`$, as function of the strength of the repulsive boson-boson and the boson-fermion interaction. For simplicity we examine the case of equal tunneling rates for the two species, i.e. $`J=J_\text{B}=J_\text{F}`$, and use $`J`$ as energy unit. The only remaining parameters of BFH Hamiltonian are then the ratios between the two interaction strengths and the tunneling rate, $`V_{\text{BB}}/J`$ and $`V_{\text{BF}}/J`$. For each point in the $`V_{\text{BB}}`$-$`V_{\text{BF}}`$-plane we solve the eigenvalue problem of the BFH Hamiltonian (1) and the two twisted Hamiltonians $`\widehat{\mathrm{H}}_{\mathrm{\Theta }_\text{B},0}`$ and $`\widehat{\mathrm{H}}_{0,\mathrm{\Theta }_\text{F}}`$ (with $`\mathrm{\Theta }_\text{B},\mathrm{\Theta }_\text{F}=0.1`$) in order to compute the phase stiffnesses and several simple ground state observables.

First we consider a lattice with $`I=8`$ sites and equal filling factors $`N_\text{B}/I=1/2`$ and $`N_\text{F}/I=1/2`$ for the bosonic and the fermionic species, respectively. Although the system is rather small we could confirm that finite size effects are insignificant for the observables discussed in the following RoBu03c . Figure 1 shows contour plots of the boson stiffness $`f_\text{s}^\text{B}`$, i.e., the superfluid fraction within the bosonic species, and the fermion stiffness $`f_\text{d}^\text{F}`$, i.e., the Drude weight for the fermionic component. In addition the largest coefficient $`C_{\mathrm{max}}^2`$ of the expansion (2) of the ground state is given. Large values of $`C_{\mathrm{max}}^2`$ indicate that the ground state is dominated by one or a few particular Fock states.

Figures 1(a) and (b) reveal that in the large parameter region for $`V_{\text{BB}}/J5`$ and $`V_{\text{BF}}/J15`$ the boson and the fermion phase stiffnesses vanish. The bosonic and the fermionic components then behave like insulators, i.e., the superfluid fraction and the Drude weight are both zero. This double-insulator phase is generated by the combined action of the boson-boson and the boson-fermion repulsion. Any configuration which has two bosons or a boson and a fermion at the same site is energetically unfavorable, and the ground state is therefore dominated by Fock states with either one boson or one fermion per site. Only if one of the interaction strengths is small can a finite stiffnesses for bosons and fermions emerge.

Closer inspection of the largest coefficient of the expansion (2), displayed in Fig. 1(c), reveals that within the double-insulating region there are two distinct areas where the largest coefficient reaches $`C_{\mathrm{max}}^20.1`$. Hence there is a class of Fock states which provide the dominant contribution to the ground state. In region (A) for values $`2V_{\text{BB}}V_{\text{BF}}`$ the dominant Fock states exhibits an alternating occupation of the sites with one boson or one fermion. In region (B) the dominant Fock states contain a continuous block of bosons followed by a block of fermions. Note that although these states have the largest weight, other Fock states do also contribute to the ground state. Most important are those which can be generated from the dominant Fock state by exchange of two lattice sites or particle-hole excitations.

The reason for the existence of and the cross-over between these two different intrinsic structures is quite subtle. All the interaction matrix elements for the dominant Fock states are irrelevant, because all these states have exactly one particle (boson or fermion) at each site. The existence of the two different structures is triggered by states with double occupancies that enter into the ground state with much smaller weight but nevertheless couple to the dominant states via the tunneling term. This complex correlation dictates the cross-over between alternating and separated species.

A good observable to distinguish the different intrinsic configurations is the static structure factor RoBu03c . For the bosonic component it is given by

$$S^\text{B}(ka)=\frac{1}{N_\text{B}^2}\underset{i,j=1}{\overset{I}{}}\mathrm{e}^{\mathrm{i}ka(ij)}\mathrm{\Psi }_0|\widehat{\mathrm{n}}_i^\text{B}\widehat{\mathrm{n}}_j^\text{B}|\mathrm{\Psi }_0,$$
(6)
where $`a`$ is the lattice spacing. It measures the density-density correlations in the system and thus provides information on the presence or otherwise of diagonal long-range order.

Figure 2 shows an example for the structure factor within the region of alternating occupation (full line) and the region of phase separation (dashed line). Apart from the trivial peaks at integer multiples of $`2\pi `$ the structure factor differs substantially within the two phases. For the alternating occupation phase, $`S^\text{B}(ka)`$ exhibits an enhancement around $`ka=\pi `$ which signals an increased probability for a boson to be found at every other site. This is the signature of diagonal long-range order as it is found in a crystal. In a two-dimensional lattice the corresponding checkerboard configuration, which exhibits crystalline diagonal long-range order, was identified recently BuBl03 .

Within the separated phase the structure factor shows a rather different behavior (dashed line in Fig. 2): it is suppressed around $`ka=\pi `$ and enhanced at $`ka2\pi /I`$. Due to this difference in the structure factors it might be possible to distinguish the two types of double-insulator phases experimentally, e.g., through selective Bragg diffraction of light off the trapped atoms. We note that the interference pattern produced by the atoms after release from the lattice, along with ballistic expansion, will not distinguish the phases as their coherence properties, i.e., the off-diagonal long-range order, are too similar.

As a second example we consider a lattice with $`I=8`$ sites, a boson filling factor $`N_\text{B}/I=1`$, and a fermion filling factor $`N_\text{F}/I=1/2`$. Contour plots of the bosonic and fermionic stiffness and the largest coefficient as function of $`V_{\text{BB}}/J`$ and $`V_{\text{BF}}/J`$ are shown in Fig. 3.

The qualitative difference compared to the half-filling case is the appearance of a bosonic Mott-insulator phase (M) for $`V_{\text{BB}}>V_{\text{BF}}`$. In this region the boson stiffness, i.e. the superfluid fraction, vanishes whereas the fermionic stiffness remains at its noninteracting value as depicted in Figs. 3(a) and (b). Under the influence of the strong boson-boson repulsion the bosons form a Mott insulating layer. The fermions are not affected, they experience a homogeneous background of bosons on which they can move freely.

The situation changes if the boson-fermion interaction becomes comparable in strength to the boson-boson interaction. Through the boson-fermion repulsion the fermions are able to break up the bosonic Mott-insulator and partially restore bosonic superfluidity in the region $`V_{\text{BF}}V_{\text{BB}}`$. Further increase of the boson-fermion repulsion leads to the simultaneous reduction of both phase stiffnesses and to the formation of a double-insulator regime as observed in the case of half-filling. The contour plot for the largest coefficient $`C_{\mathrm{max}}^2`$ depicted in Fig. 3(c) reveals that both of the previously observed intrinsic configurations appear again: (A) an alternating occupation of the sites associated with crystalline diagonal long-range order and (B) continuous blocks of bosons and fermions which resemble a spatial separation of the two components. Due to the appearance of the bosonic Mott-insulator both phases have been shifted to the region $`V_{\text{BF}}>V_{\text{BB}}`$.

In summary, we have mapped out the phase diagrams of a binary mixture of bosonic and fermionic atoms in a one-dimensional lattice. By exact numerical solution of the eigenvalue problem for the Bose-Fermi-Hubbard Hamiltonian we are able to compute ground states as well as the dynamical response of the system. The stiffnesses of the bosonic and fermionic component under phase twists provides an important measure for the superfluid properties of the bosonic component (superfluid fraction) and the conductivity of the fermionic species (Drude weight). Using these quantities we have identified several distinct quantum phases, e.g., a bosonic Mott-insulator phase and two double-insulator phases which exhibit different degrees of diagonal long-range order—alternating occupation and component separation. This shows that atomic boson-fermion mixtures in optical lattices constitute a versatile tool for the study of fundamental quantum phase transitions.

This work was supported by the DFG, the UK EPSRC, and the EU via the “Cold Quantum Gases” network. K.B. thanks the Royal Society and Wolfson Foundation for support.