text:
# Oscillator strength distribution in C3H6 isomers studied with the time-dependent density functional method in the continuum

## 1 Introduction

The photoabsorption and photoionization cross sections of molecules are of significant interest in many fields of both fundamental and applied sciences. The oscillator strength distribution characterizing the optical response is the most important quantity in understanding the interaction of photon with electrons in atoms, molecules, and matters. The oscillator strength distribution in the whole spectral region has been extensively studied with the advanced synchrotron radiation and the high resolution electron energy loss spectroscopy . In order to see how the oscillator strength changes with varying molecular structures, it is useful to study isomer molecules. Since the isomers consist of the same kind and the same number of atoms, we expect similarity of the oscillator strengths at high photon energies. This is because the molecular structure has little influence on the excitation of inner core electrons. However, the valence photoabsorption may differ according to the difference of electronic states between the isomers. In fact, Koizumi et al. have observed a prominent distinction for the cross sections of simple hydrocarbon isomers, C3H6 (cyclopropane and propylene), at photon energies of $`1020`$ eV . The photoabsorption and photoionization data of cyclopropane in this energy region were later improved by a measurement with metallic thin film windows . These works clearly show that the oscillator strengths have different peaks and shoulders depending on the isomers in a continuous spectral region above the ionization potentials (IPs).

Theoretical investigation for the isomer effect of C3H6 has been demanded for a long time, however, none has been reported so far. This is due to difficulties in treatment of the electronic continuum in a non-spherical multicenter potential. There are several methods which are able to take into account correlations among valence electrons in the continuum . Nevertheless, some of them are not suitable for calculating detailed structures of the spectra, and some are difficult to be applied to large molecules. We have recently developed an alternative theoretical method for this purpose . The method is based on the time-dependent local-density approximation (TDLDA) in a grid representation of the three-dimensional Cartesian coordinate, and utilizes the Green’s function to take account of the continuum boundary condition. We use iterative methods to solve linear algebraic equations for construction of the dynamical screening field above the first IP. The theoretical background of our method is similar to the one of Ref.  in which the authors used a single-center expansion technique. However, the application was limited to small axially symmetric molecules, because of difficulties in the single-center expansion. Our method is based on direct calculation of the self-consistent screening potential in the three-dimensional grid representation, which does not rely on the expansion and requires no spatial symmetry.

In the present Letter, we report the valence photoabsorption of the C3H6 isomers studied with the continuum TDLDA method in Ref. , and would like to show the power of the method. Then, we give an interpretation of the continuous spectra and elucidate origins of the isomer effects.

## 2 Theory and computational method

Optical response of molecules is characterized by the oscillator strength, denoted as $`df/d\omega `$ in the followings, which is given by

$$\frac{df}{d\omega }=\frac{2m\omega }{3\pi }\mathrm{Im}\underset{\nu =x,y,z}{}d^3rr_\nu \delta n_\nu (𝐫,\omega ),$$
(1)
where the transition density $`\delta n_\nu `$ is related to the Fourier component of a time-dependent external dipole perturbation in $`\nu `$-direction, $`V_\nu (\omega )=r_\nu `$, through a complex susceptibility ($`\delta n_\nu (\omega )=\chi (\omega )V_\nu (\omega )`$).

The TDLDA describes a spin-independent $`N`$-electron system in terms of the time-dependent Kohn-Sham equations. Correlations among electrons are taken into account through deformations of the self-consistent Kohn-Sham potential. Linearizing the Kohn-Sham potential with respect to the transition density, we obtain

$$\begin{array}{cc}\hfill \delta n_\nu (𝐫,\omega )& =d^3r^{}\chi _0(𝐫,𝐫^{};\omega )\hfill \\ & \left\{r_\nu ^{}+d^3r^{\prime \prime }\frac{\delta V_{\mathrm{KS}}[n(𝐫^{})]}{\delta n(𝐫^{\prime \prime })}\delta n_\nu (𝐫^{\prime \prime })\right\}.\hfill \end{array}$$
(2)
The $`\chi _0(𝐫,𝐫^{};\omega )`$ is a complex susceptibility for a system without the correlations and is given by

$$\begin{array}{cc}\hfill \chi _0(𝐫,𝐫^{};\omega )=2\underset{i}{\overset{\mathrm{occ}}{}}& \varphi _i(𝐫)\{\left(G(𝐫,𝐫^{};ϵ_i\omega ^{})\right)^{}\hfill \\ & +G(𝐫,𝐫^{};ϵ_i+\omega )\}\varphi _i(𝐫^{}).\hfill \end{array}$$
(3)
Here, $`\varphi _i`$’s are the ground-state Kohn-Sham orbitals and $`G`$ is the Green’s function for an electron in the static Kohn-Sham potential. In order to properly treat the electronic continuum, the outgoing boundary condition must be imposed on $`G`$. Construction of the Green’s function is easily done for a rotationally invariant potential $`V_0(r)`$, using the partial wave expansion . Thus, we split the Kohn-Sham potential into two parts, a long-range spherical part, $`V_0(r)`$, and a short-range deformed part, $`\stackrel{~}{V}(𝐫)=V_{\mathrm{KS}}(𝐫)V_0(r)`$. First, we construct the Green’s function, $`G_0`$, for the spherical potential $`V_0`$, then, $`G`$ can be obtained from an identity

$$G=G_0+G_0\stackrel{~}{V}G.$$
(4)
We solve Eqs. (2), (3), and (4) simultaneously in the uniform grid representation of the three-dimensional real space. Equations (2) and (4) are linear algebraic equations with respect to $`\delta n`$ and $`G`$, respectively, for which an iterative method provides an efficient algorithm for the numerical procedure. We adopt the generalized conjugate residual method for these non-hermitian problems. We would like to refer the reader to Ref.  for detailed discussion of the methodology and the theoretical background.

The exchange-correlation potential is a sum of the local density part given by Ref. , $`\mu ^{(\mathrm{PZ})}[\rho ]`$, and the gradient correction of Ref. , $`\mu ^{(\mathrm{LB})}[\rho ,\rho ]`$, which will be abbreviated to LB potential. This gradient correction is constructed so as to reproduce the correct Coulomb asymptotic behavior of the potential ($`e^2/r`$) and to describe the Rydberg states. It was also pointed out that the LB potential is necessary to reproduce the excitation energies of high-lying bound states in the TDLDA . In our previous work , we have also found for simple molecules that the TDLDA with the LB potential reasonably accounts for resonances embedded in the continuum.

## 3 Results and discussion

### 3.1 Ground-state properties

We fix the geometry of nuclei optimized for the ground state. This is based on a semiempirical method known as PM3 . We only treat valence electrons in the TDLDA calculation. Thus, we use the norm-conserving pseudopotential  with a separable approximation  for the electron-ion potentials. The coordinate space is discretized in a square mesh of 0.3 Å and we adopt all the grid points inside a sphere of 6 Å radius. This results in a model space of 33,401 grid points.

First, we calculate the ground state of cyclopropane and propylene by solving the Kohn-Sham equations with the exchange-correlation potential of $`\mu ^{(\mathrm{PZ})}+\mu ^{(\mathrm{LB})}`$. The LB potential $`\mu ^{(\mathrm{LB})}`$ contains a parameter $`\beta `$ , and we adjust this value to make eigenvalues of the highest occupied molecular orbitals (HOMO) coincide with the empirical vertical IPs ($`10.54`$ eV for cyclopropane  and $`9.91`$ eV for propylene ). The occupied valence orbitals in the ground state calculated with $`\beta =0.015`$ are listed in Table 1. The HOMO eigenvalues are well reproduced for both isomers. Propylene has a geometry of the $`C_s`$ point group while cyclopropane has the $`D_{3h}`$ group. Although these isomers possess equal number of valence electrons (eighteen valence electrons), the electron configuration of cyclopropane is more degenerate in energy because of the higher symmetry.

### 3.2 Photoabsorption oscillator strength

Now we calculate the photoresponse of the isomers. We use complex frequencies, $`\omega +i\mathrm{\Gamma }/2`$ with $`\mathrm{\Gamma }=0.5`$ eV. The $`\mathrm{\Gamma }`$ plays a role of a smoothing parameter to make the energy resolution finite. This also helps a convergence of the numerical iteration procedure . In Figs. 1 and 2 respectively, we show the calculated photoabsorption oscillator strength, $`df/d\omega `$, for cyclopropane and propylene in a frequency (photon energy) range of $`850`$ eV. The calculations have been done with a frequency mesh of $`\mathrm{\Delta }\omega =0.25`$ eV for a region of $`8\omega 20`$ eV, and with $`\mathrm{\Delta }\omega =0.5`$ eV for the rest of frequencies. The oscillator strength distributions of the isomers are nearly identical at $`\omega 22`$ eV. This energy roughly corresponds to the ionization energy of the lowest-lying $`\sigma `$ orbital. The $`df/d\omega `$ monotonically decreases as the frequency increases but has a large tail at high frequency. This behavior of the high-frequency tail in $`df/d\omega `$ is one of the characteristics of the electronic excitations in the continuum, which was also found in our previous studies of simple molecules . The molecular structure has a little influence on the electronic continuum in this energy region ($`\omega 22`$ eV).

In contrast, in the frequency region below 22 eV, different structures are observed among the isomers. The $`df/d\omega `$ of propylene shows a single broad peak at $`\omega =1318`$ eV with small wiggles. On the other hand, distinctive three peaks at $`\omega =11.8`$, 13.5, and 15.7 eV, are found in cyclopropane. This difference exactly matches the experimental findings of the isomer effect (the thin solid line in Fig. 1 (a)). The energy positions of calculated peaks are lower than the experimental ones by about 1.5 eV. This is also true for propylene in Fig. 2 (a), in which the broad peak is shifted to lower energy by 1.5 eV compared to the experiment. We would like to mention that, if we treat the electrons as responding independently to the external dipole field, we cannot reproduce the main feature of the oscillator strength distributions. We call this approximation “independent-particle approximation (IPA)”, which corresponds to neglecting the induced screening potential, the second term in the bracket in Eq. (2). In Figs. 1 (a) and 2 (a), the IPA calculations are shown by dashed lines.

The Thomas-Kuhn-Reiche (TRK) sum rule tells us the integrated oscillator strength $`f(\mathrm{})=18`$ in our calculations, since we only treat valence electrons in C3H6. In Table 2, partial sum values of the oscillator strengths are listed and compared to the experiment. Again, the IPA calculation cannot account for the data, while the TDLDA in the continuum well agrees with the experiment. Calculated total sum values for $`8<\omega <60`$ eV are 16.2 for cyclopropane and 15.8 for propylene. These values correspond to about 90 % of the TRK sum rule for valence electrons.

### 3.3 Origin of the difference in photoabsorption between the isomers

We would like to discuss details of resonance peaks and an origin of the different behaviors between the isomers. First, let us compare the IPA results for the two isomers (See dashed lines in Figs. 1 (a) and 2 (a)). In the energy region of $`1020`$ eV, although the bulk structure is similar, cyclopropane shows a sharper main peak at 11.5 eV and an additional peak at 13 eV. This may be due to higher-fold degeneracies in electronic eigenstates in cyclopropane. These peak structures in cyclopropane remain after inclusion of the dynamical screening effects, while, for propylene, the strong peak at 12 eV seen in the IPA is diminished. Next, we shall examine this difference in the dynamical screening effects.

For this purpose, it is useful to calculate a partial oscillator strength  which corresponds to a contribution of each occupied orbital to the total oscillator strength. We display the partial $`df/d\omega `$ in the energy range of $`1025`$ eV in Figs. 1 (b) and 2 (b). One can see that the major contributions to $`df/d\omega `$ come from bound-to-continuum excitations of electrons near the Fermi level; the HOMO $`(3e^{})^4`$ and the second HOMO $`(1e^{\prime \prime })^4`$ in case of cyclopropane, and the second $`(10a^{})^2`$, the third $`(9a^{})^2`$ and the fourth HOMO $`(1a^{\prime \prime })^2`$ for propylene.

The sharp peaks at 11.8 eV and 13.5 eV in cyclopropane originate from bound-to-bound transitions of $`(3a_1^{})^2`$ and $`(1a_2^{\prime \prime })^2`$ electrons. These resonances are also seen in the IPA calculation. Then, the electron-electron correlation brings out coherent contributions of bound-to-continuum excitations of $`(3e^{})^4`$ and $`(1e^{\prime \prime })^4`$ electrons. The width of the resonances becomes slightly larger than that of the IPA, because of the autoionization process. The peak at 15.7 eV is also produced by coherent excitations of $`(2e^{})^4`$ (bound-to-bound) and $`(3e^{})^4`$ (bound-to-continuum) electrons. We only see a shoulder around 15.5 eV in the IPA calculation, however, the dynamical effect enhances the peak.

In the case of propylene, we find several small peaks of bound-to-bound excitations in the partial $`df/d\omega `$ in Fig. 2 (b). However, the bound-to-continuum excitations, which mostly contribute to the broad resonance in $`1318`$ eV, behave rather independently, and do not produce coherent enhancement of those peaks. As a result, the small peaks in the bound-to-bound transitions are mostly smeared out in the total oscillator strength.

We think that this difference in the continuum response could be attributed to the difference in strength of bound-to-bound transitions. The $`df/d\omega `$ of propylene shows a typical behavior of the dynamical screening effects. Namely, the oscillator strengths (and the peak at 12 eV) in the IPA calculation are significantly weakened by the induced screening field in Eq. (2). Conversely, those at energies above 16 eV are enhanced. This is because the real part of the dynamical polarizability changes its sign at $`\omega 16`$ eV, then the screening field changes into the “anti-screening field” at higher energies. In the case of cyclopropane, the situation is slightly more complicated. Generally speaking, the real part of the dynamical polarizability changes its sign from positive to negative at bound resonances. Since the degeneracies of electronic orbitals are higher in cyclopropane, the bound-to-bound transitions have large oscillator strengths. Then, there appears a reminiscence of bound resonance in the continuum region. The screening field suddenly drops down at energies corresponding to bound-to-bound transitions. This provides effective anti-screening effects to cause the peak structures in the bound-to-continuum transitions.

Finally, we would like to comment on dependence of our results upon the parameter $`\beta `$ in the LB potential $`\mu ^{(\mathrm{LB})}`$. A choice of this parameter is rather arbitrary, since the value of $`\beta `$ does not change the Coulomb asymptotic behavior. In fact, if we adopt $`\beta =0.05`$, the value proposed by the original paper , we obtain better agreement to the photoabsorption spectra for both the isomers. In the case of using $`\beta =0.05`$, the Kohn-Sham eigenvalues for occupied orbitals in cyclopropane are indicated as values in brackets in Table 1. A binding energy of each orbital becomes deeper by $`1.31.5`$ eV, though spacings between the orbitals almost stay constant. The calculated oscillator strength distribution is shown in Fig. 3 for cyclopropane. The disagreement on the peak positions are removed in the calculation. A bound peak at 10 eV is also well reproduced in the calculation. This peak consists of the excitations of $`(3e^{})^4`$ and $`(1e^{\prime \prime })^4`$ electrons. These excitations have an almost identical energy when we use $`\beta =0.05`$, while the excitation of $`(1e^{\prime \prime })^4`$ electrons is shifted to lower energy by 1 eV when using $`\beta =0.015`$. Except for the bound peak at 10 eV, the characteristic features of the oscillator strength distribution are not changed, and we maintain the interpretation given above.

## 4 Conclusion

The oscillator strength distributions of C3H6 isomer molecules are studied with the TDLDA in the continuum utilizing the three-dimensional Cartesian coordinate representation. The calculation shows good agreement with experiments. The oscillator strength in the energy region above 22 eV is almost identical among the isomers, however, different peaks appear below that. This isomer effect is analysed by calculating the partial oscillator strength of each occupied orbital. In addition to the difference in properties of bound electronic orbitals, it turns out that bound-to-continuum excitations of electrons near the Fermi level behave differently between the isomers. The bound-to-bound transitions in cyclopropane possess large strengths, and the bound-to-continuum transitions exhibit coherent peak structures because of the anti-screening effects. On the other hand, in propylene, the bound-to-bound transitions are too weak to produce the anti-screening peaks for the continuum excitations. Although the molecular structure directly has minor influence on the continuum, the difference in bound-to-bound transitions leads to variation in the dynamical screening effects to change the continuum excitations.

Acknowledgements

This work is supported in part by Grants-in-Aid for Scientific Research (No.1470146 and 14540369) from the Japan Society for the Promotion of Science. Calculations were performed on a NEC SX-5 Super Computer at Osaka University and a HITACHI SR8000 at Institute of Solid State Physics, University of Tokyo.


format: text

text:
# Tunneling broadening of vibrational sidebands in molecular transistors

## I Introduction

In recent years, the field of molecular electronics has renewed the interest in transport through mesoscopic systems with strong electron-phonon coupling. A number of experiments have been reported to demonstrate transport through single molecules,reed97 ; park00 ; park02 ; lian02 ; smit02 ; zhit02  some of which show signs of vibrational sidebands. In the pioneering experiments by Park et al.park00  it was shown that electron transport through a single $`C_{60}`$ molecule was strongly influenced by a single vibrational mode. The single phonon mode was associated with the motion of the molecule in the confining potential created by the van der Waals interactions with the electrodes. Later similar devices with more complicate molecules were investigatedpark02 ; zhit02  and they also showed excitation spectra which possibly could be associated with emission of vibrational quanta, so-called Frank Condon peaks in the differential conductance. Some of these devicespark02 ; lian02  furthermore exhibited a peak stemming from the Kondo effect proving that the tunnel coupling to the leads was rather strong.

Theoretically there has been a large amount of work on the problem of tunneling through a single level with coupling to phonon modes. The approaches fall into two categories. The first category is the kinetic equation approach, which is relevant in the weak tunneling limit.rate ; flen03  In the kinetic equation approach, it is essential that the excited vibrational levels are allowed to relax through coupling to a bath. For a large phonon relaxation rate one can assume an equilibrium phonon distribution, otherwise a non-equilibrium distribution function must be determined from the kinetic equations.rate  The coupling to the dissipative environment of the molecule leads to additional broadening of phonon sidebands.flen03

The second category, in which also this paper falls, deals with the strong tunneling limit. The first approaches were an exact solution of a simplified situation, where only one electron is present,glaz88shek ; wing89  i.e., the presence of the Fermi sea is ignored. The result of this exact solution is that the phonon sidebands in the tunneling density of states all are Lorentzians with the same width. The exact solution amounts to decoupling the electron and oscillator displacement operators. This decoupling approximation has been used in a number of papers.kuochang02 ; lundin02 ; alex02  Below we refer to this approximation as the “single particle approximation” (SPA) and we demonstrate explicitly that this approximation is only valid for high energies, i.e., for electrons or holes sufficiently far from the Fermi surface, or in the weak tunneling limit where the rate equation is valid. For the many-body system many decay channels are in fact blocked by the Pauli principle. Other many-body works that include Fermi sea effects have been donekral97 ; koni96 ; keil02  and the blocking of the line width has also been alluded to.kral97 ; koni96

In this paper, we focus on the broadening of the Frank-Condon peaks in the differential conductance due to tunneling broadening, but in the non-Kondo regime. While this case cannot be solved exactly, an approximation, which is believed to be valid above the Kondo temperature and which includes the Pauli principle, but ignores correlations between lead electrons, is developed. The main result is that the combination of the Pauli principle and the splitting of the spectral weight into phonon sidebands severely limits the phase space for tunnel broadening. Hence the peaks are much narrower than the bare tunneling rate would suggest. We find a simple analytic result that describes this, see Eq. (38) and compare with the SPA formula in Eq. (25). The approximation is exact in both the single particle limit and in the weak tunneling limit. Furthermore, we also support our conclusion by a perturbation theory in the electron-phonon coupling derived in Appendix B.

## II Model Hamiltonian and
current formula

We consider a model of a single electron level coupled to two leads. The single level is coupled to a vibrational mode of the molecule through the onsite energy. For simplicity we ignore the spin degree of freedom, which is not relevant unless Kondo-type effects are important. The model Hamiltonian is

$$H=H_k+H_D+H_{DB}+H_B+H_T,$$
(1)
where

$`H_k`$ $`={\displaystyle \underset{k,\alpha =L,R}{}}\xi _{k,\alpha }c_{k\alpha }^{}c_{k\alpha },H_D=\xi _0d^{}d,`$ (2)
$`H_\mathrm{B}`$ $`={\displaystyle \frac{p_0^2}{2m_0}}+{\displaystyle \frac{1}{2}}m_0\omega _0^2x_0^2,H_{DB}=\lambda x_0d^{}d,`$ (3)
$`H_T`$ $`={\displaystyle \underset{k,\alpha =L,R}{}}t_{k\alpha }\left(c_{k\alpha }^{}d+d^{}c_{k\alpha }\right).`$ (4)
Here $`x_0`$ is the oscillator degree of freedom and $`\xi _0`$ is the bare onsite energy. The lead electron energies are given by $`\xi _{k\alpha }=\epsilon _{k\alpha }\mu _\alpha `$, where $`\mu _\alpha `$ is the chemical potential of the lead $`\alpha `$. The tunneling amplitude between lead $`\alpha `$ and the molecule is $`t_{k\alpha }`$. The amplitudes could in principle also depend on the oscillator position. However, in the experimental realizations so far this is a negligible dependence and it is omitted here. Because of this, we can apply the current formula derived by Wingreen and Meirmeir92  who expressed the current in terms of the non-equilibrium Green’s function of the molecule. The current through the system is

$$I=\frac{e}{h}\frac{d\xi }{2\pi }\frac{\mathrm{\Gamma }_L\mathrm{\Gamma }_R}{\mathrm{\Gamma }_L+\mathrm{\Gamma }_R}\left[n_L(\xi )n_R(\xi )\right]A(\xi ),$$
(5)
where the left and right distribution functions are

$$n_\alpha (\epsilon )=n_F(\epsilon eV_\alpha ),$$
(6)
and where $`n_F(ϵ)=(e^{\beta ϵ}+1)^1`$ is the usual Fermi-Dirac function. The spectral function, $`A`$, is given by

$$A(\xi )=2\mathrm{Im}G^R(\xi ),$$
(7)
where the retarded Green’s, $`G^R`$, functions are

$$G^R(t,t^{})=i\theta (tt^{})\{d(t),d^{}(t^{})\}.$$
(8)
The widths of the level due to coupling to the leads are

$$\mathrm{\Gamma }_\alpha (\xi )=2\pi \underset{k}{}t_{ka}^2\delta (\xi \xi _{k\alpha }).$$
(9)
The Green’s function $`G^R`$ must be calculated in non-equilibrium and in the presence of the leads.

It is often useful to eliminate the linear coupling terms of the Hamiltonian Eq. (1) by a unitary transformation similar to one used in the independent boson model.mahan  The transformation is

$$\stackrel{~}{H}=SHS^{},S=e^{ip_0\mathrm{}d^{}d},\mathrm{}=\frac{\lambda }{m_0\omega _0^2},$$
(10)
and with this choice the Hamiltonian transforms into

$$\stackrel{~}{H}=H_k+H_B+\stackrel{~}{H}_D+\stackrel{~}{H}_T,$$
(11)
where

$$\stackrel{~}{H}_D=\epsilon _0d^{}d,\epsilon _0=\xi _0\frac{1}{2}\lambda \mathrm{},$$
(12)
and

$$\stackrel{~}{H}_T=\underset{k\alpha }{}t_{k\alpha }\left(c_{k\alpha }^{}e^{ip_0\mathrm{}}d+d^{}e^{ip_0\mathrm{}}c_{k\alpha }\right).$$
(13)
The retarded Green’s function becomes after this transformation

$$G^R(tt^{})=i\theta (tt^{})\{e^{ip_0(t)\mathrm{}}d(t),d^{}(t^{})e^{ip_0(t^{})\mathrm{}}\}_S,$$
(14)
where the average of course should be taken with respect to $`\stackrel{~}{H}`$, which is indicated by $`\mathrm{}_S.`$

## III Calculation of the
spectral function

In this section, we calculate the spectral function, $`A(\xi )`$ that enters into the current formula, Eq. (5). As mentioned in introduction, the tunneling broadening of phonon assisted side bands has been considered before in the case when the presence of the Fermi sea is ignored(SPA). In this case the model can be solved exactly, at least in the so-called wide band limit (WBL). glaz88shek ; wing89  For latter reference we start by quoting the SPA result

$$G^{R,\mathrm{SPA}}(t)=\mathrm{exp}\left(t\mathrm{\Gamma }/2\right)G_0^R(t),$$
(15)
where $`G_0^R`$ is the Green’s function in the absence of tunneling. We see from this expression that the SPA implies that all conductance steps are smeared by the same amount.

### III.1 Dyson equation for $`G^R`$

In the following, we develop a method to calculate the broadened Green’s function using a truncated equation of motion technique. Our starting point is the retarded Green’s function in Eq. (14). We expand it in terms of eigenstates of the boson systems

$`G^R(t)`$ $`={\displaystyle G_{nn^{},mm^{}}^R(t)f_{nn^{}}f_{mm^{}}^{}},`$ (16)
$`G_{nn^{},mm^{}}^R(t)`$ $`=i\theta (t)[(|nn^{}|d)(t),d^{}|m^{}m|]_S,`$ (17)
where we have defined the overlap function between the oscillator states

$$f_{nn^{}}=n|e^{ip_0\mathrm{}}|n^{}.$$
(18)
An expression for $`f_{nn^{}}`$ is given in Eq. (51).

We now generate a series of equations of motion, starting with one for $`G_{nn^{},mm^{}}^R`$

$`(i_t+E_{nn^{}}\epsilon _0)G_{nn^{},mm^{}}^R(t)=L_{nn^{}}^R(t)`$
$`+\delta (t)(dd^{}|nm|_S\delta _{n^{}m^{}}+d^{}d|m^{}n^{}|_S\delta _{nm}),`$ (19)
where $`E_{nn^{}}=(nn^{})\omega _0`$, and $`L_{nn^{}}^R`$ is the contribution stemming from the tunneling Hamiltonian. The function $`L^R`$ is (see Appendix A for details of the derivations)

$`L_{nn^{}}^R(t)`$ $`={\displaystyle \underset{k\alpha ,l}{}}t_{k\alpha }f_{ln^{}}^{}(G_{k\alpha ,nl,mm^{}}^R(t)G_{n_dk\alpha ,nl,mm^{}}^R(t))`$
$`+{\displaystyle \underset{k\alpha ,l}{}}t_{k\alpha }f_{nl}^{}G_{n_dk\alpha ,ln^{},mm^{}}^R(t),`$ (20)
where we encounter two new Green’s functions

$`G_{k\alpha ,nl,mm^{}}^R(t)`$ $`=i\theta (t)[(|nn^{}|c_{k\alpha })(t),d^{}|m^{}m|]_S,`$ (21)
$`G_{n_dk\alpha ,nl,mm^{}}^R(t)`$ $`=i\theta (t)[(|nn^{}|d^{}dc_{k\alpha })(t),d^{}|m^{}m|]_S.`$ (22)
The equation of motion for the first one is

$`(i_t+E_{nl}\xi _{k\alpha })G_{k\alpha ,nl,mm^{}}^R(t)=L_{k\alpha ,nl}^R(t)`$
$`+\delta (t)(c_{k\alpha }d^{}|nm|_S\delta _{lm^{}}+d^{}c_{k\alpha }|m^{}l|_S\delta _{nm}),`$ (23)
where again $`L_{k\alpha ,nl}^R`$ is the contribution from the tunneling term. Again the tunneling term generates new higher order Green’s functions, but at this point we truncate it using the following physical principles: we neglect correlations involving lead electrons, which means that we decouple terms lead electron operators using the Hartree-Fock approximation. Furthermore, we set $`c_{k\alpha }d^{}0`$ in the equation of motion for (21) and (22). With these approximations (see Appendix A for more details)

$`L_{k\alpha ,nl}^R(t)`$ $`t_k{\displaystyle \underset{j}{}}n_{k\alpha }f_{jn}G_{jl,mm}^R`$
$`+t_{k\alpha }{\displaystyle \underset{j}{}}(1n_{k\alpha })f_{lj}G_{nj,mm}^R`$ (24)
It is now worth to note that in the case with only one electron, which is the SPA, the term $`n_k`$ is exactly zero and Eq. (24) becomes exact for this case. In the same limit, $`G_{n_dk\alpha ,nl,mm^{}}^R=0,`$ and the equations are easily solved. After setting the result back into Eq. (16) we get

$$G^{R,\mathrm{SPA}}(\omega )=\underset{nn^{}}{}|f_{nn^{}}|^2\frac{(1\overline{n})N_n+\overline{n}N_n^{}}{\omega +E_{nn^{}}\epsilon _0+i\mathrm{\Gamma }/2},$$
(25)
which is nothing but the SPA result in Eq. (15) written in the oscillator eigenstate basis. In doing this we have furthermore evaluated the last term in Eq. (19) as $`dd^{}|nm|_S\delta _{n^{}m^{}}+d^{}d|m^{}n^{}|_S\delta _{nm}=[(1\overline{n})N_n+\overline{n}N_n^{}]\delta _{n^{}m^{}}\delta _{nm},`$ where $`N_n=|nn|`$ is the occupation of the $`n`$’th oscillator level and $`\overline{n}=d^{}d`$ is the average level occupation.

Here we want to go beyond the SPA and we continue by looking at the equation of motion for the last Green’s function in Eq. (22). We have

$`(i_t+E_{nl}\xi _{k\alpha })G_{n_dk\alpha ,nl,mm^{}}^R(t)=`$
$`L_{n_dk\alpha ,nl}^R(t)+\delta (t)c_{k\alpha }d^{}|nm|\delta _{lm^{}}.`$ (26)
The function $`L_{n_dk\alpha ,nl}^R(t)`$ becomes after doing the same line of approximations as was done for $`L_{k\alpha ,nl}^R(t)`$ (see Appendix A)

$$L_{n_dk\alpha ,nl}^R(t)t_{k\alpha }n_{k\alpha }\underset{j}{}f_{jn}G_{jl,mm^{}}^R.$$
(27)
The set of equations now close. After integrating out the Green’s function involving the lead electrons, we arrive at linear equations for the molecule Green’s function

$`(\omega +E_{nn^{}}\epsilon _0)G_{nn^{},mm^{}}^R(\omega )=`$
$`dd^{}|nm|_S\delta _{n^{}m^{}}+d^{}d|m^{}n^{}|_S\delta _{nm}+`$
$`{\displaystyle \underset{lj}{}}\left(\mathrm{\Lambda }_{ln^{}}^ef_{nl}^{}f_{jl}G_{jn^{},mm^{}}^R(\omega )+\mathrm{\Lambda }_{nl}^hf_{ln^{}}^{}f_{lj}G_{nj,mm^{}}^R(\omega )\right),`$ (28)
where

$$\mathrm{\Lambda }_{nl}^{e,h}(\omega )=\underset{\alpha }{}𝑑\xi \frac{\mathrm{\Gamma }_\alpha }{2\pi }\left(\begin{array}{c}n_\alpha (\xi )\\ 1n_\alpha (\xi )\end{array}\right)\frac{1}{\omega +E_{nl}\xi +i\eta }.$$
(29)
The functions $`\mathrm{\Lambda }^{e,h}`$ are in the case of energy independent $`\mathrm{\Gamma }_\alpha `$’s (which we assume henceforth) given by

$`\mathrm{\Lambda }_{nl}^e(\omega )`$ $`={\displaystyle \underset{\alpha }{}}\left[{\displaystyle \frac{i\mathrm{\Gamma }_\alpha }{2}}n_\alpha (\omega +E_{nl})+{\displaystyle \frac{\mathrm{\Gamma }_\alpha }{2\pi }}\psi (\omega +E_{nl})\right],`$
$`\mathrm{\Lambda }_{nl}^h(\omega )`$ $`={\displaystyle \underset{\alpha }{}}\left[{\displaystyle \frac{i\mathrm{\Gamma }_\alpha }{2}}n_\alpha (\omega E_{nl}){\displaystyle \frac{\mathrm{\Gamma }_\alpha }{2\pi }}\psi (\omega +E_{nl})\right],`$ (30)
where

$$\psi (\epsilon )=𝑑\xi n_F(\xi )\left[\frac{1}{\epsilon \xi }+\frac{1}{\xi }\right],$$
(31)
where the last term has been added as a regularization of the integral at large energies. This is allowed because it cancels when adding the two $`\psi `$-terms in (28) (which is easily seen by noting that $`_lf_{nl}^{}f_{jl}=\delta _{n,l}`$). Note that $`\psi (\epsilon )`$ is a function of $`\epsilon /kT`$ only and that $`\psi (\epsilon )=\psi (\epsilon )`$. At large values of $`\epsilon /kT`$ the $`\psi `$-function has the asymptotic form $`\psi (\epsilon )\mathrm{ln}(\epsilon /kT)`$, while for small $`\epsilon `$ it goes as $`\psi (\epsilon )(\epsilon /kT)^2\times 0.213139\mathrm{}`$.

The only remaining question is how to evaluate the last term in Eq. (19). In absence of tunneling or in the single particle approximation, we can decouple the electron and phonon degrees of freedom, such that $`dd^{}|nm|dd^{}N_n\delta _{nm},`$ where $`N_n=|nn|`$ is the occupation of the $`n`$’th oscillator level and $`\overline{n}=d^{}d`$ is the average level occupation. To leading order in $`\mathrm{\Gamma }`$ this is a reasonable approximation. However, in the more general case we should in principle solve the terms $`dd^{}|nm|`$ and $`d^{}d|m^{}n^{}|`$ on first right hand of (28) self-consistently. In the general non-equilibrium case, it is, however, not possible to calculate these expectation values from the retarded Green’s function alone. One must therefore invoke new approximations for this purpose. Instead of pursuing that, we concentrate on a case where this is not necessary, namely the situation of a strongly asymmetric device. This is done in the following section.

### III.2 Asymmetric transistor, $`\mathrm{\Gamma }_L\mathrm{\Gamma }_R`$

In the strongly asymmetric thermal case, things simplify considerably, because the dot states are in thermal equilibrium with one of the leads. For the case $`\mathrm{\Gamma }_L\mathrm{\Gamma }_R`$, the current is

$$I=\frac{e}{h}\frac{d\xi }{2\pi }\mathrm{\Gamma }_R\left[n_L(\xi )n_R(\xi )\right]A_L(\xi ),$$
(32)
where $`A_L`$ is the equilibrium spectral function in equilibrium with the left lead. Therefore, we can use the standard equilibrium expression for the self-consistent equations

$`\delta _{nm}d^{}d|m^{}n^{}|`$ $`={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}A_{nn^{},mm^{}}(\omega )n_F(\omega ),`$ (33)
$`\delta _{n^{}m^{}}dd^{}|nm|`$ $`={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}A_{nn^{},mm^{}}(\omega )n_F(\omega ),`$ (34)
where $`A_{nn^{},mm^{}}(\omega )=i(G_{nn^{},mm^{}}^R(\omega )G_{d,nn^{},mm^{}}^A(\omega ))`$.

We have numerically solved these equations and in Fig. 1 we plot the result for the spectral function, $`A_L(\xi )`$, for $`\epsilon _0=0`$, and different values of the tunnel broadening. It is clearly seen how the Fermi surface effects sharpen the lines as compared to the SPA. This effect is more pronounced for the first peaks, which is also evident from Fig. 2, where we show the spectral function for the off-resonance condition.

### III.3 Approximate solution of the Dyson equations

In order to gain more physical understanding of the narrowing of the lines seen above and to obtain analytical results, we now proceed to solve the equations (28) approximately. As is shown in the numerical plots, the approximate analytic solution, that is derive below, is in fact close to the full solution.

We look for a solution of Eq. (28) near one of the resonances, $`\omega \epsilon _0+n\omega _0`$. Because the Green’s function $`G_{d,nn^{},mm^{}}^R(\omega )`$ peaks at $`\omega =\epsilon _0E_{nn^{}}`$, the two terms in the sum in Eq. (28) are dominated by, $`E_{jn^{}}E_{nn^{}}`$ and $`E_{nj}E_{nn^{}}`$, respectively. Hence we can set $`j=n`$ in the first one and $`j=n^{}`$ in second one, which then gives the following diagonal equations

$`(\omega +E_{nn^{}}\epsilon _0\mathrm{\Sigma }_{nn}(\omega ))G_{nn^{},mm^{}}^R(\omega )=`$
$`\delta _{nm}\delta _{n^{}m^{}}\{(1\overline{n})N_n+\overline{n}N_n^{}\},`$ (35)
where the self-energy is

$$\mathrm{\Sigma }_{nn^{}}(\omega )=\underset{l}{}\left(\mathrm{\Lambda }_{ln^{}}^e(\omega )|f_{nl}|^2+\mathrm{\Lambda }_{nl}^h(\omega )|f_{ln^{}}|^2\right).$$
(36)
We have, furthermore, approximated

$`dd^{}|nm|_S\delta _{n^{}m^{}}+d^{}d|m^{}n^{}|_S\delta _{nm}`$
$`\delta _{nm}\delta _{n^{}m^{}}\{(1\overline{n})N_n+\overline{n}N_n^{}\},`$ (37)
where $`N_n=\mathrm{exp}(n\beta \omega _0)`$ and $`\overline{n}=d^{}d`$. This approximation is valid for not too large $`\mathrm{\Gamma }`$ and we have verified numerically, in the self-consistent calculation performed above, that it is indeed reasonable accurate for $`\mathrm{\Gamma }/\omega _01`$. With these approximations, we obtain a closed expression for the Green’s function $`G^R`$

$$G^R(\omega )=\underset{nn^{}}{}\frac{(1\overline{n})N_n+\overline{n}N_n^{}}{\omega +E_{nn^{}}\epsilon _0\mathrm{\Sigma }_{nn}(\omega )}|f_{nn^{}}|^2.$$
(38)
The result in Eq. (38) thus generalizes the SPA approximation, Eq. (25), by including the leading order effect of the Fermi sea. The new self-energy in Eq. (36) has a simple physical interpretation: the broadening, which is caused by tunneling out and tunneling in processes, can only occur if the state in the lead is either empty or occupied, respectively. For the non-interacting case, where we in Eq. (36) set $`E_{nn^{}}=0`$ and use that $`_l|f_{nl}|^2=1`$, the spectral functions reduces correctly to an Lorentzian with width $`\mathrm{\Gamma }`$.

We can now in detail study the narrowing of the on-resonance line seen in Fig. 1. At resonance we have $`\overline{n}=\frac{1}{2}`$ and we then straightforwardly find

$$A(\omega )\frac{|f_{00}|^2\mathrm{\Gamma }_{00}^2}{\omega ^2+(\mathrm{\Gamma }_{00}/2)^2},\mathrm{\Gamma }_{00}=|f_{00}|^2\mathrm{\Gamma }.$$
(39)
This result shows that the width of the resonance at $`\epsilon _0=0`$ is considerably narrower than the bare tunnel broadening would suggest. For large values of $`g`$ $`|f_{00}|^2`$ becomes smaller and the lines narrows in. However, also the weight of the line goes down with $`|f_{00}|^2`$ as is evident from Eq. (39).

### III.4 Comparison with perturbation theory in $`g`$

It is interesting to compare the results where the tunneling is treated approximate but the electron-phonon coupling exactly, to the results of a perturbation theory, where the tunneling is treated exactly but the electron-phonon coupling to lowest order. This situation is solved in Appendix B and Fig. 3 shows the result of the calculation. Note that both methods give qualitative the same result, namely that the zero energy line is narrowed and that spectral function is suppressed between the two peaks.

## IV Conclusions

We have presented results for transport through quantum dots with strong electron-phonon coupling and with strong tunneling broadening of the phonon sidebands. This is not a solvable problem even when the spin degree of freedom is ignored, and we resorted to an approximation which incorporates the Fermi sea, but ignores correlation effects. The approximation is exact in both the single particle and weak tunneling cases. However, physical conclusions can be drawn from the approximate approach, namely that the tunnel broadening is much weaker than expected from a model where the Pauli principle is not incorporated. We have, as mentioned, neglected the spin degree of freedom which is not expected to be important at temperatures above the Kondo temperature. It is, however, an interesting question to ask how the Kondo effect is influenced by the coupling between the electron occupancy and the vibrational mode of the molecule.

###### Acknowledgements.

The author is thankful for stimulating discussions with P. Brouwer, A.-P. Jauho, P. McEuen, and J. Sethna and for financial support from the Danish National Research Council.

## Appendix A Derivation of the equations of motion

Here we present the commutators and the methods used to derive the set of equations of motion in Section III.1. We write the Hamiltonian in Eq. (11) as $`H=H_0+H_T`$. In the oscillator eigenstate representation $`H_T`$ is

$$H_T=\underset{k\alpha ,ll^{}}{}t_{k\alpha }\left(d^{}c_{k\alpha }|l^{}l|f_{ll^{}}^{}+|ll^{}|f_{ll^{}}c_{k\alpha }^{}d\right).$$
(40)
Then Eqs. (19),(20) follow from

$`[H_0,d|nn^{}|]`$ $`=(E_{nn^{}}\epsilon _0)d|nn^{}|,`$ (41)
$`[H_T,d|nn^{}|]`$ $`={\displaystyle \underset{k\alpha ,ll^{}}{}}t_kf_{ll^{}}^{}\left(d^{}dc_{k\alpha }|l^{}n^{}|\delta _{nl}dd^{}c_{k\alpha }|nl|\delta _{n^{}l^{}}\right).`$ (42)
The equation of motion in Eq. (23) follows from the commutators

$`[H_0,c_{k\alpha }|nl|]`$ $`=(E_{nl}\xi _{k\alpha })c_{k\alpha }|nl|,`$ (43)
$`[H_T,c_{k\alpha }|nl|]`$ $`={\displaystyle \underset{k^{}\alpha ^{}jj^{}}{}}t_{k^{}\alpha ^{}}[f_{jj^{}}(|nj^{}|\delta _{k\alpha ,k^{}\alpha ^{}}\delta _{jl}+c_{k^{}\alpha ^{}}^{}c_{k\alpha }(|jl|\delta _{j^{}n}|nj^{}|\delta _{jl}))d`$
$`+f_{jj^{}}^{}(|j^{}l|\delta _{jn}|nj|\delta _{j^{}l})c_{k^{}\alpha ^{}}c_{k\alpha }d^{}].`$ (44)
At this point we truncate the equations of motion approximating the terms with 3 electron operators by the Hartree-Fock decomposition for the lead electrons. We have

$`c_{k^{}\alpha ^{}}^{}c_{k\alpha }d`$ $`c_{k^{}\alpha ^{}}^{}c_{k\alpha }dc_{k^{}\alpha ^{}}^{}dc_{k\alpha },`$ (45)
$`c_{k^{}\alpha ^{}}c_{k\alpha }d^{}`$ $`c_{k^{}\alpha ^{}}c_{k\alpha }d^{}c_{k\alpha }c_{k^{}\alpha ^{}}d^{}.`$ (46)
As mentioned in the main text, we will at this level neglect the terms $`c_{k\alpha }d^{}`$ and therefore only $`c_{k^{}\alpha ^{}}^{}c_{k\alpha }d`$ contributes and gives $`c_{k^{}\alpha ^{}}^{}c_{k\alpha }d\delta _{k\alpha ,k^{}\alpha ^{}}c_{k\alpha }^{}c_{k\alpha }d`$. The last approximate sign here means up to order $`\mathrm{\Gamma }`$, which is consistent with neglecting $`c_{k\alpha }d^{}.`$ When this is inserted back into (44), we get

$$[H_T,c_{k\alpha }|nl|]\underset{jj^{}}{}t_{k\alpha }\left[f_{jj^{}}\left(|nj^{}|\delta _{jl}+c_{k\alpha }^{}c_{k\alpha }(|jl|\delta _{j^{}n}|nj^{}|\delta _{jl})\right)d\right],$$
(47)
which then leads to Eq. (24). Finally, the last equation of motion in Eq. (23) follows from

$`[H_0,d^{}dc_{k\alpha }|nl|]`$ $`=(E_{nl}\xi _{k\alpha })d^{}dc_{k\alpha }|nl|,`$ (48)
$`[H_T,d^{}dc_k|nl|]`$ $`={\displaystyle \underset{k^{}\alpha ^{},jj^{}}{}}t_{k^{}\alpha ^{}}\left[f_{jj^{}}|jl|c_{k^{}\alpha ^{}}^{}c_{k\alpha }d\delta _{j^{}n}f_{jj^{}}^{}|nj|c_{k\alpha }c_{k^{}\alpha ^{}}d^{}\delta _{j^{}l}\right],`$ (49)
and after using the same arguments as before, we obtain

$$[H_T,d^{}dc_{k\alpha }|nl|]\underset{j}{}t_kf_{jn}|jl|c_{k\alpha }^{}c_{k\alpha }d,$$
(50)
which leads to Eq. (27). Furthermore, when the oscillator states are chosen real the function $`f_{nn^{}}`$ is

$$f_{nn^{}}=\sqrt{\frac{1}{2^n2^n^{}n!n^{}!}}e^{g/2}\left[\mathrm{sign}(n^{}n)\sqrt{\frac{g}{2}}\right]^{|nn^{}|}2^{\mathrm{max}(n,n^{})}L_{|nn^{}|}^{\mathrm{min}(n,n^{})}(g).$$
(51)
## Appendix B Perturbation theory in $`g`$

In this appendix we derive a perturbative result for the spectral function to first order in the coupling constant but to all orders in the tunneling matrix element. We do, however, also assume a constant bare tunneling density of states and also take the asymmetric case, $`\mathrm{\Gamma }_R\mathrm{\Gamma }_L=\mathrm{\Gamma }`$. We want to compute $`G^R`$ in Eq. (8) and for the perturbative calculation it is more convenient to use the Hamiltonian (1) instead of the transformed Hamiltonian. We will start from the equations of motion for the Green’s function $`G_d^R`$. We use here the notation $`G_y^R=i\theta (t)[y(t),d^{}].`$ We obtain

$$\left(\omega \xi _0+i\frac{\mathrm{\Gamma }}{2}\right)G_d^R(\omega )=1+\lambda G_{xd}^R(\omega ).$$
(52)
The equation of motion for the function $`G_{xd}^R`$ is

$$\left(\begin{array}{cc}\omega \xi _0& i/m\\ im\omega _0^2& \omega \xi _0\end{array}\right)\left(\begin{array}{c}G_{xd}^R\\ G_{pd}^R\end{array}\right)=\left(\begin{array}{c}x\\ 0\end{array}\right)+\underset{k}{}t_k\left(\begin{array}{c}G_{xk}^R\\ G_{pk}^R\end{array}\right)+\frac{\lambda }{2}\left(\begin{array}{c}\mathrm{}_0^2\\ i\end{array}\right)G_{d,\lambda =0}^R(\omega ),$$
(53)
and

$$\left(\begin{array}{cc}\omega \xi _k+i\eta & i/m\\ im\omega _0^2& \omega \xi _k+i\eta \end{array}\right)\left(\begin{array}{c}G_{xk}^R\\ G_{pk}^R\end{array}\right)=t_k\left(\begin{array}{c}G_{xd}^R\\ G_{pd}^R\end{array}\right)i\lambda \left(\begin{array}{c}0\\ 1\end{array}\right)G_{nk,\lambda =0}^R(\omega ),$$
(54)
where we used that since $`G_{xd}^R`$ in (52) is multiplied by $`\lambda `$, the last term of (54) should be calculate without electron phonon coupling. Therefore, for temperatures $`kT\omega _0`$ the following replacements were made in that term $`\lambda x^2d=\lambda x^2_0d=\frac{1}{2}\mathrm{}_0^2d`$, and $`\lambda pxd=\lambda px_0d=\frac{i}{2}d`$. The average value $`x`$ appearing in Eq. (53) follows from $`\ddot{x}=\omega _0^2xn_d\lambda /m=0`$. Expressing the functions $`G_{xk}^R`$ and $`G_{pk}^R`$ in Eq. (54) in terms of $`G_{xd}^R`$ and $`G_{pd}^R`$ and inserting this back into Eqs. (53) and (52), we obtain

$$G_d^R(\omega )=\frac{1g}{N_0}\frac{g\omega _0}{N_0^2}+\frac{gn}{N_{+1}}+\frac{g\left(1n\right)}{N_1}+\frac{\lambda ^2}{m}\frac{1}{N_{+1}N_1N_0}\underset{k}{}\frac{t_k\left(N_0+\omega \xi _k\right)}{\left(\omega \xi _k+i\eta \right)^2\omega _0^2}G_{nk,\lambda =0}^R(\omega ),$$
(55)
where $`N_p=\omega \xi _0+p\omega _0+i\frac{\mathrm{\Gamma }}{2}`$. The first four terms in (55) are nothing but the SPA result in Eq. (25) expanded to lowest order in $`g`$. *Therefore, it is evident that the last term of Eq. (55) constitutes the correction to the SPA*. The Green’s function $`G_{nk,\lambda =0}^R`$ can be evaluated using Wick’s theorem

$$G_{nk,\lambda =0}^R(\omega )=n_d_0G_{k,\lambda =0}^R(\omega )d^{}c_k_0G_{d,\lambda =0}^R(\omega ),$$
(56)
where

$$G_{k,\lambda =0}^R(\omega )=t_k\frac{1}{\omega \xi _k+i\eta }G_d^R(\omega ),d^{}c_k_0=i\frac{d\omega ^{}}{\pi }n_F(\omega ^{})\left(G_{k,\lambda =0}^R(\omega ^{})G_{k,\lambda =0}^A(\omega ^{})\right).$$
(57)
The integrant in Eq. (55) is analytic in the upper half plane of the complex $`\xi _k`$ plane, and when performing the integral by a contour in the lower half-plane only one term in Eq. (57) contributes. The last term in Eq. (55), which was the correction to the SPA result, then reads

$$\delta G_d^R(\omega )=\frac{g\omega _0^2\mathrm{\Gamma }}{N_{+1}N_1N_0^2}\gamma (\omega )$$
(58)
where

$$\gamma (\omega )=\omega _0_{\mathrm{}}^{\mathrm{}}\frac{d\omega ^{}}{\pi }n_F(\omega ^{})\frac{1}{\omega ^{}\xi _0i\mathrm{\Gamma }/2}\frac{2\omega \xi _0\omega ^{}}{\left(\omega ^{}\omega i\eta \right)^2\omega _0^2}.$$
(59)
The SPA is restored when the energy is either far above or far below the Fermi surface, because the model is then similar to a single electron or a single hole. This is clearly seen in the limit $`|\xi _0|,\omega _0\mathrm{max}(\mathrm{\Gamma },kT)`$, where we can replace $`n_F`$ by either one or zero. For the first cases it is clear that $`\delta G`$ is 0. When $`n_F=1`$, we have also get zero because the integrant is analytic in one half plane.

Finally, at $`T=0,`$ we have

$`\gamma (\omega )`$ $`={\displaystyle \frac{\omega _0}{2\pi }}{\displaystyle \frac{2\omega 2\xi _0i\mathrm{\Gamma }/2}{(\omega \xi _0i\mathrm{\Gamma }/2)^2\omega _0^2}}\left(2i\pi n_0+\mathrm{ln}\left[(\xi _0/\omega _0)^2+(\mathrm{\Gamma }/2\omega _0)^2\right]\right)`$
$`+{\displaystyle \frac{\omega \xi _0\omega _0}{2\pi (\omega \xi _0+\omega _0i\mathrm{\Gamma }/2)}}\left(i\pi \theta (\omega _0\omega )+\mathrm{ln}\left|(\omega _0+\omega )/\omega _0\right|\right)`$
$`{\displaystyle \frac{\omega \xi _0+\omega _0}{2\pi (\omega \xi _0\omega _0i\mathrm{\Gamma }/2)}}\left(i\pi \theta (\omega _0\omega )+\mathrm{ln}\left|(\omega _0\omega )/\omega _0\right|\right),`$
where

$$n_0=\frac{1}{\pi }\mathrm{tan}^1\left(\frac{2\xi _0}{\mathrm{\Gamma }}\right)+\frac{1}{2}.$$
(60)
The zero temperature expression for the SPA correction is seen to diverge at $`|\omega |=\omega _0`$. This divergence is, however, cut-off at finite temperatures. If Fig. 3 we show the result of the perturbation theory compared to the approximate spectral function derived in the main text. The plot is at the symmetric point $`\epsilon _0=\xi _0g\omega _0=0`$. One should note that when expanding $`\epsilon _0`$ to linear order it cancels off the term $`g\omega _0/N_0^2`$ in Eq. (55).


text:
# Constraints on the equation of state of dark energy and the Hubble constant from stellar ages and the CMBPartly based on observations collected at ESO (La Silla and Paranal) under programmes 62.O-0592, 63.O-0468, and 64.O-0281, 65.O-0446 and 66.A-0362, and at the W. M. Keck Observatory, which is operated jointly by the California Institute of Technology and the University of California.

## 1 Introduction

There is now solid observational evidence that the universe is accelerating. The absolute ages of high-redshift galaxies (Dunlop et al., 1996; Spinrad et al., 1997) and the value of the Hubble constant $`H_0`$ (Freedman et al., 2001), the Hubble diagram of Type Ia supernovae (Riess et al., 1998; Perlmutter et al., 1999), the angular power spectrum of the cosmic microwave background (CMB; Benoit et al., 2002; Ruhl et al., 2002; Goldstein et al., 2002) in combination with galaxy surveys (Efstathiou et al., 2002; Verde et al., 2002) and lensing (Wang et al., 2002), all indicate that the universe is flat and dominated at present by some form of dark energy with negative pressure. Strong confirmation for this comes from the newly released results from the Wilkinson Microwave Anisotropy Probe (WMAP ) first year data (Spergel et al., 2003).

The equation of state of the dark energy, $`p=w\rho `$, expresses the ratio between the pressure, $`p`$, and the mass density, $`\rho `$, of the dark energy in terms of the parameter $`w`$ (in units of $`c=1`$). The value of $`w`$ could either be constant, as in the case of a cosmological constant ($`w=1`$), or be time-dependent, as in the case of a rolling scalar field or “quintessence” (Peebles & Ratra, 1988; Caldwell et al., 1998). Any such behavior would have far-reaching implications for particle physics.

In this paper we use the ages of the oldest stellar populations between $`z=0`$ and $`z=1.5`$ to derive new and independent constraints on two key cosmological parameters, $`w`$ and $`H_0`$. Both parameters are relatively poorly constrained by CMB data alone. For example, in order to get accurate constraints on $`w`$, Spergel et al. (2003) combined the WMAP data with six external data sets: two additional CMB data sets (CBI and ACBAR; Pearson et al., 2002; Kuo et al., 2002), large scale structure measurements from the 2dF galaxy redshift survey (2dFGRS; Percival et al., 2001; Verde et al., 2002), measurements of the Lyman $`\alpha `$ forest power spectrum (Croft et al., 2002; Gnedin & Hamilton, 2002), the luminosity distance–redshift relation from Type Ia supernovae measurements (Garnavich et al., 1998; Riess et al., 2001), and $`H_0`$ from the Hubble Space Telescope (HST ) Key Project (Freedman et al., 2001). Therefore it is particularly important to obtain accurate and independent measurements of cosmological parameters relying on different physics.

In Section 2 we focus on $`w`$, summarizing the basis of our method and deriving tight constraints on the equation of state using the ages of globular clusters in combination with the location of the first acoustic peak in the CMB power spectrum. Then, in Section 3, we determine stellar ages for a sample of early-type galaxies between $`z=0`$ and $`z=1.5`$ and use the age–redshift relationship defined by the old envelope to determine $`H_0`$. The prospects of using this method to derive a measurement of $`w(z)`$ are also discussed. Section 4 summarizes and discusses the results.

Throughout the paper we adopt the following notation. The Hubble constant is $`H_0=100h`$ km s-1 Mpc-1 and $`c`$ is the speed of light. The density of baryons, matter, radiation, photons, neutrinos, and dark energy, at present, in critical units, are indicated as $`\mathrm{\Omega }_b`$, $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_{\mathrm{rad}}`$, $`\mathrm{\Omega }_\gamma `$, $`\mathrm{\Omega }_\nu `$, and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ respectively ($`\mathrm{\Omega }_{\mathrm{rad}}\mathrm{\Omega }_\gamma +\mathrm{\Omega }_\nu `$). The physical quantities, independent of $`h`$, are indicated with $`\omega `$; e.g., $`\omega _m=\mathrm{\Omega }_mh^2`$. When needed, a standard, flat, reference $`\mathrm{\Lambda }`$CDM model with $`\mathrm{\Omega }_m=0.27`$, $`h=0.71`$, $`\mathrm{\Omega }_bh^2=0.023`$, and $`w=1`$ is assumed (Spergel et al., 2003).

## 2 The value of $`w`$ from stellar ages

In this section we use the ages of the oldest globular clusters in combination with the location of the first acoustic peak of the CMB power spectrum (Page et al., 2003) to constrain the value of $`w`$. The method is briefly described in Section 2.1 and the results are given in Section 2.2.

### 2.1 Method

If the universe is assumed to be flat, the position of the first acoustic peak ($`\mathrm{}_1`$ in the standard spherical harmonics notation) depends primarily on the age of the universe and on the effective value111I.e., the average value over redshift. We use $`w(z)`$ to indicate when $`w`$ is allowed to vary over time. of $`w`$ (Caldwell et al., 1998; Hu et al., 2001; Knox et al., 2001). As noted by these authors, for a fixed $`w`$ value, a change in the physical density parameter $`\omega _m=\mathrm{\Omega }_mh^2`$ that keeps the characteristic angular scale of the first acoustic peak fixed will also leave the age approximately unchanged. Thus an independent estimate of the absolute age of the universe at $`z=0`$ combined with a measurement of $`\mathrm{}_1`$ yields an estimate of $`w`$ largely independent of other cosmological parameters (e.g., $`\mathrm{\Omega }_m`$ and $`h`$).

This can be better understood by considering that for a constant $`w(z)`$, the age of a flat universe is given by

$$t_0=H_0^1_0^{\mathrm{}}(1+z)^1\left[\mathrm{\Omega }_{\mathrm{rad}}(1+z)^4+\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }(1+z)^{3(1+w)}\right]^{1/2}𝑑z.$$
(1)
The position of the first acoustic peak is fixed by the quantity

$$\theta _A=r_s(a_{dec})/d_A(a_{dec}),$$
(2)
where $`a_{dec}`$ is the scale factor at decoupling, $`r_s(a_{dec})`$ is the sound horizon at decoupling, and $`d_A(a_{dec})`$ is the angular diameter distance at decoupling. For a flat universe,

$$r_s(a_{dec})=\frac{c}{H_0\sqrt{3}}_0^{a_{dec}}\left[\left(1+\frac{3\mathrm{\Omega }_b}{4\mathrm{\Omega }_\gamma }\right)(\mathrm{\Omega }_\mathrm{\Lambda }x^{13w}+\mathrm{\Omega }_mx+\mathrm{\Omega }_{\mathrm{rad}})\right]^{1/2}𝑑x$$
(3)
and

$$d_A(a_{dec})=\frac{c}{H_0}_{a_{dec}}^1\left[\mathrm{\Omega }_\mathrm{\Lambda }x^{13w}+\mathrm{\Omega }_mx+\mathrm{\Omega }_{\mathrm{rad}}\right]^{1/2}𝑑x.$$
(4)
See Verde et al. (2003) for more details. We use equations 1 to 4 and fix $`\theta _A`$ to be consistent with the WMAP determination ($`\mathrm{}_1220`$; Page et al., 2003). Figure 1 shows the allowed region in the age-$`w`$ plane obtained for $`\mathrm{\Omega }_bh^2=0.02`$, $`0.05<\mathrm{\Omega }_m<0.4`$, and $`0.5<h<0.9`$. It is clear from the plot that an independent and accurate age determination of the universe will provide a measurement of $`w`$. Thus by using only the WMAP observation of $`\mathrm{}_1`$ and an independent estimate of the age of the universe, we can place constraints on $`w`$, largely independent of other cosmological parameters.

### 2.2 Results from the ages of globular clusters

The ages of the oldest globular clusters (GCs) provide a lower limit to the total age of the universe. Numerous star-forming galaxies have now been observed up to redshift $`z=6.6`$ (e.g., Kodaira et al., 2003). Since the oldest GCs contain the oldest stellar populations in galaxies, it is reasonable to assume that GCs too have formed at redshift $`z>6`$. Furthermore, since small halos are easily destroyed after reionization (e.g., Bullock et al., 2001), GCs likely formed prior to the reionization epoch. Under this assumption, the WMAP recent estimate of the reionization redshift $`z_{\mathrm{reion}}17`$ (Kogut et al., 2003) implies that the ages of the oldest globular cluster do not significantly underestimate the age of the universe.

Globular cluster ages can reliably be estimated in several different ways: from the luminosity of the main-sequence turn-off (if the distance is known), from detached binaries, from the termination of the white-dwarf sequence, and from the luminosity function method which provides a joint constraint on the age and distance of a GC. Renzini et al. (1996) determine the distance to NGC 6752 by means of accurate photometry of GC white dwarfs and obtained an age of $`14.5\pm 1.5`$ Gyr. Using the same method they find 47 Tucanae to have an age of $`13\pm 2.5`$ Gyr (Zoccali et al., 2001). Gratton et al. (1997), using the new Hipparcos distance scale to determine main-sequence turn-offs, estimate an age of $`11.8_{2.5}^{+2.1}`$ Gyr (95% confidence level) for the oldest Galactic GCs. Jimenez & Padoan (1998) apply the luminosity function method to M55 and obtain an age of $`12.5\pm 1.0`$ Gyr (95% confidence level). Hansen et al. (2002), using very deep HST images, find the white dwarf cooling sequence for M4 and estimate its age to be $`12.7\pm 0.7`$ Gyr at the $`95`$% confidence level. Using the theoretical age–mass relation and applying it to the detached eclipsing binary OGLEGC-17, Chaboyer & Krauss (2002) determine Omega Centauri to have an age of $`11.2\pm 1.1`$ Gyr (95% confidence level). These age estimates are all remarkably consistent. Krauss & Chaboyer (2003) perform the most careful analysis to date of the effects of systematics in GC age determinations. They estimate the age of the oldest Galactic GCs using the main-sequence turn-off luminosity and evaluate the errors with Monte-Carlo techniques, paying careful attention to uncertainties in the distance, to systematics, and to model uncertainties. They find a best-fit age for the oldest Galactic GCs of $`12.5_{2.2}^{+3.4}`$ Gyr (95% confidence limits). We adopt their probability distribution for the oldest GC age, which we find can be accurately described by

$$P(t)=\frac{A}{\sigma (tT)}\mathrm{exp}\left[\frac{\mathrm{ln}[(tT)/m]^2}{2\sigma ^2}\right]$$
(5)
where $`t`$ denotes the age of the oldest GCs, $`A=1.466`$, $`\sigma =0.25`$, $`T=6.5`$ and $`m=5.9`$. Since reionization took place at $`z17`$ (Kogut et al., 2003) and we assume that oldest GCs must have formed before then, for all reasonable cosmologies we only need add $`0.20.3`$ Gyr to the GC ages to obtain an estimate of the age of the universe. We conservatively add 0.3 Gyr222The uncertainty in this correction is negligible compared to the uncertainty in the determination of GC ages..

In order to constrain $`w`$, we assume a flat universe and Monte Carlo simulate the distribution of $`\mathrm{}_1`$ subject to only weak constraints on the other cosmological parameters. For different values of $`w`$, we generate $`10^5`$ models, randomly sampling the cosmological parameters $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_bh^2`$ and $`h`$ with uniform priors, $`0.5<h<0.9`$, $`0.05<\mathrm{\Omega }_m<0.45`$, and $`0.01<\mathrm{\Omega }_bh^2<0.03`$. We then impose an age of the universe constraint by randomly sampling these models so that the age of the universe has a probability distribution whose shape is given by equation 5, but offset by $`0.3`$ Gyr. We use the publicly available code CMBFAST (Seljak & Zaldarriaga, 1996) to compute the resulting distribution of $`\mathrm{}_1`$. This is shown in Figure 2 where the two solid lines are the 68% and 90% joint confidence levels.

As expected, the age alone does not constrain $`w`$; there is a degeneracy between $`w`$ and $`\mathrm{}_1`$. If we now use the observed position of the first acoustic peak as recently measured, in a model independent way, from WMAP (Page et al., 2003), we can constrain $`w`$ with high accuracy. The filled contours in Figure 2 are marginalized values for $`w`$ at the 68% and 90% confidence levels. Thus we find $`w<0.8`$ ($`w<0.67`$) at the 68% (90%) confidence level. If we drop the assumption that the oldest GC must have formed before reionization, but still impose that they must have formed by $`z=6`$, this constraint is slightly weakened. In this case the age distribution of equation 5 needs to be offset by $`1`$ Gyr and we obtain $`w<0.5`$ at the 90% confidence level. We have also repeated the calculation using Gaussian priors for the cosmological parameters, $`h=0.7\pm 0.2`$, $`\mathrm{\Omega }_m=0.25\pm 0.2`$, $`\mathrm{\Omega }_bh^2=0.02\pm 0.01`$, and find similar, but slightly more stringent constraints on $`w`$.

We stress here that this determination depends solely on the GC determination of the age of the universe and on the observed position of the first acoustic peak in the CMB power spectrum. Our constraint is slightly less stringent than that obtained by Spergel et al. (2003) from a joint likelihood analysis of WMAP with six external data sets (WMAP \+ CBI + ACBAR + 2dFGRS + Lyman$`\alpha `$ forest power spectrum + Type Ia supernovae + $`H_0`$ constraint from the HST key project), but is tighter than the CMB-only (WMAP \+ CBI + ACBAR) determination and comparable to the WMAP \+ ACBAR + CBI + HST constraint.

## 3 $`H_0`$ and $`w(z)`$ from the ages of unresolved stellar populations at $`z=01.5`$

### 3.1 Method

We now turn our attention to the age of the oldest objects in the universe at $`z0`$ as a function of redshift. As recently suggested by Jimenez & Loeb (2002), these can be used to determine $`H(z)=\frac{1}{1+z}\frac{dz}{dt}`$ and therefore to measure $`w(z)`$ and $`H_0`$. The Hubble parameter is given by $`H^2(z)=H_0^2[\rho _{\mathrm{tot}}(z)/\rho _{\mathrm{tot}}(0)]`$, where $`\rho _{\mathrm{tot}}`$ is the total density of the universe. For a flat universe composed of matter and dark energy, $`H(z)`$ is related to the equation of state by

$$H_0^1\frac{dz}{dt}=(1+z)\frac{H(z)}{H_0}=(1+z)^{5/2}\left[\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }\mathrm{exp}\left\{3_0^z\frac{w(z^{})}{(1+z^{})}𝑑z^{}\right\}\right]^{1/2},$$
(6)
where we have used the conservation equation for dark energy density, $`\dot{\rho _\mathrm{\Lambda }}=3H(z)[1+w(z)]\rho _\mathrm{\Lambda }`$. Thus a measurement of $`dz/dt`$ is also a measurement of $`w(z)`$; for further details, see Jimenez & Loeb (2002).

To obtain $`dz/dt`$, Jimenez & Loeb (2002) proposed to use the old envelope of the age–redshift relation of E/S0 galaxies. We assume that E/S0 at different redshifts are drawn from the same parent population with the bulk of their stellar populations formed at relatively high redshift (e.g., Bower et al., 1992; Stanford et al., 1998). Then, at relatively low redshift, they are evolving passively and may be used as “cosmic chronometers”.

In the rest of this section we will apply this method to a sample of spectra (described in §3.2) to identify the age–redshift relation defined by the old envelope (§3.3) and derive an independent measurement of $`H_0`$ (§3.4). We show that at the level of accuracy of our measurement, the use of the age–redshift relation allows us to apply this method even if a fraction of the stars in E/S0 galaxies is assembled at $`z<1`$ (e.g., Treu et al., 2002). We also briefly discuss the prospects of using this method to derive $`w(z)`$.

### 3.2 Sample and measurement

In order to apply this method we need high-quality spectra of a sample of old stellar populations covering the largest possible range in redshift. To this aim we combine the following data sets: (i) the luminous red galaxy (LRG) sample from the Sloan Digital Sky Survey (SDSS) early data release (Eisenstein et al., 2001); (ii) the sample of field early-type galaxies from Treu et al. (1999, 2001, 2002, hereafter the Treu et al. sample); (iii) a sample of red galaxies in the galaxy cluster MS1054$``$0321 at $`z=0.83`$; and (iv) the two radio galaxies 53W091 and 53W069 (Dunlop et al., 1996; Spinrad et al., 1997; Nolan et al., 2003, Dey et al., in preparation).

The LRG sample is the most heterogenous, since it includes red galaxies in general, regardless of morphology. The Treu et al. sample contains high quality spectra with average $`S/N=16`$ per Å. The spectra in the Treu et al. sample have higher $`S/N`$ than the typical LRG spectrum, but the wavelength coverage is better for the LRG sample (see example in Figure 3). To select potentially passively evolving galaxies, we include only galaxies with $`S/N>9`$ per Å and for which passively evolving synthetic stellar population models (Jimenez et al., 1998) provide a good fit to the continuum. This selection criterion discards dusty or star-forming objects. The spectra of seven red galaxies in the cluster MS1054$``$0321 were obtained by one of us (DS) on UT 11 March 2002 using the dual-beam Low Resolution Imaging Spectrometer (LRIS; Oke et al., 1995) at Keck Observatory. On the blue side, we used the 300 lines mm-1 grism ($`\lambda _{\mathrm{blaze}}=5000`$ Å; spectral resolution $`\mathrm{\Delta }\lambda _{\mathrm{FWHM}}13.5`$ Å) and on the red side, we used the the 400 lines mm-1 grating ($`\lambda _{\mathrm{blaze}}=8500`$ Å; spectral resolution $`\mathrm{\Delta }\lambda _{\mathrm{FWHM}}8`$ Å). These observations, which were obtained at the parallactic angle, totaled one hour of integration. In order to ensure no gaps in the spectra, half the data were obtained with a D560 dichroic and the remainder were obtained with a D680 dichroic. These spectra have a $`S/N`$ of 8. We also add the red radio galaxies 53W091 and 53W069 for which accurate ages have already been determined (Dunlop et al., 1996; Spinrad et al., 1997; Nolan et al., 2003, Dey et al., in preparation). Examples of the spectra used are given in Figure 3.

The age of the dominant stellar population in the galaxies is obtained by fitting single stellar population models (Jimenez et al., 1998, 2003a) to the observed spectrum. These models have been extensively tested (Kotoneva et al., 2002; Jimenez et al., 2003b) and have two free parameters: age and metallicity. By construction, star formation is assumed to occur in a single burst of duration much shorter ($`<10`$%) than the current age of the galaxy (Jimenez et al., 1999). Therefore, the derived ages are single-burst equivalent ages. The best fit model is found by standard $`\chi ^2`$ minimization. Then, stellar ages and their errors are obtained by marginalizing the resulting likelihood surface over metallicity. Some examples for the parameter confidence regions, at the 1 and 2$`\sigma `$ level, are shown in Figure 3. Note that for galaxies with high $`S/N`$, like the Treu et al. sample, the constraints in both age and metallicity are very tight and thus the so-called age-metallicity degeneracy is dramatically reduced and does not affect the age determination at a significant level.

Before we can use the estimated ages to constrain $`H_0`$ we have to estimate the effects of the potential systematic uncertainties related to adopting a single-burst stellar population model and a single-metallicity model. As far as the first uncertainty is concerned, we know that even massive field E/S0 galaxies can experience some level of star formation activity at $`z>0.2`$ (Menanteau et al., 2001; Treu et al., 2002; Willis et al., 2002). The presence of a young and luminous stellar population superimposed on a dominant old stellar population can bias the single-burst equivalent ages of the spectrum to be much younger than the mass-weighted average age. Although this is a serious concern for the interpretation of individual cases, our method relies on the old envelope of the age distribution at any given redshift. Therefore, the effects of this bias on the determination of $`H_0`$ and $`w`$ are strongly reduced, as we discuss further in § 3.2. Concerning single metallicity, the uncertainty connected to metallicity gradients is negligible with respect to other sources of uncertainty. In fact, our spectra typically cover $`6`$ kpc (for the LRG) to $`10`$ kpc (for the higher redshift samples). Given that typical metallicity gradients in E/S0s follow $`d\mathrm{log}[\mathrm{Fe}/\mathrm{H}]/d\mathrm{log}r=0.2`$ (Davies et al., 1993), metallicity only changes by 0.2 dex within this region. This is fairly small and the single-metallicity approximation affects our results well within the 1$`\sigma `$ value of the error in the recovered metallicity.

### 3.3 An age–redshift relation and $`w(z)`$

Figure 4 shows the derived single-burst equivalent ages of the galaxies as a function of redshift. The circles correspond to galaxies in the SDSS LRG sample, triangles to the Treu et al. sample, diamonds are galaxies in MS1054$``$0321, and crosses are 53W091 and 53W069. For clarity, the SDSS LRG points have not been plotted for ages $`<7`$ Gyr. Typical errors on the ages of LRG galaxies are 10% and are not plotted. An age–redshift relation (i.e., an “edge” or “envelope” of the galaxy distribution in the age–redshift plane) is apparent from $`z=0`$ to $`z=1.5`$.

Although an envelope can be identified, it is important to notice that there is a large spread of ages at any given redshift. In particular, the measured age for many galaxies makes them significantly younger than the age–redshift envelope. This is in part due to a real spread in the star formation histories of early-type galaxies, especially in the field, but mostly a result of using single-burst equivalent ages. Single-burst equivalent ages effectively measure a luminosity-weighted age. Therefore, even a relatively minor episode of star formation is sufficient to significantly underestimate the age of the bulk of (old) stars. Indeed, many of the galaxies in the Treu et al. sample show \[O II\] $`\lambda `$3727 emission, consistent with minor ongoing star formation activity (Treu et al., 2002). Is this contamination enough to bias our results, or is the use of the relative ages of the old-envelope sufficient to reduce this bias to useful levels?

To answer this question we consider a scenario where the bulk of stars in E/S0 galaxies formed at high redshift, while a relatively small fraction (up to 10% in mass) of stars is formed at later times ($`z<1.5`$). This scenario is approximately consistent with the observed evolution of the number density (Im et al., 2002), colors (Jimenez et al., 1999; Menanteau et al., 2001), and spectral properties of early-type galaxies (Trager et al., 2000; Treu et al., 2002). Within this scenario, we use Monte Carlo simulations to estimate the effects of this bias on the age–redshift relation and ultimately on $`H_0`$ and $`w(z)`$ determinations. We generate galaxy spectra assuming that after the initial burst of star formation at $`z>3`$ an extra 10% of the stellar mass is formed at a time $`t_1`$ between redshift $`0<z<1.5`$ with uniform probability for $`t_1`$. For simplicity, we use discrete values for the ages with spacing 0.1 Gyr and assume a flat $`\mathrm{\Omega }_m=0.27`$, $`w=1`$ and $`h=0.71`$ cosmology to convert time to redshift. The metallicity for a given galaxy is chosen to be uniformly distributed in the range $`0.1<Z/Z_{}<4`$. Finally, we add random noise to the spectra so that $`S/N=15`$ and determine the age with the same procedure used for the data.

The results of our Monte Carlo simulations are shown in Figure 5. The top panel shows the simulated age–redshift distribution if only six early-types are found per redshift interval333We have drawn galaxies ages randomly from a grid that spans the age range 0.1 to 14 Gyr, with a spacing 0.1 Gyr. For our standard, baseline cosmology, this age spacing translates to the redshift interval seen in Figure 4.. The edge is blurred and smeared due to recent star formation contamination and because of the error in the age determination; nevertheless, it is clearly recognizable. The galaxies observed within a few 100 Myr of bursts have very small single-burst equivalent ages, as observed for some galaxies in the LRG and Treu et al. sample, but they do not interfere with the location of the old envelope.

The lower panel shows the distribution for the case where 60 early-types are observed per redshift bin (this is less than the number of galaxies in the SDSS LRG sample). The solid line shows the theoretical position of the edge, which is remarkably close to observed edge. Note also that increasing the number of galaxies has improved the agreement between the observed and theoretical edge. Thus we conclude that small, recent bursts of star formation do not introduce a significant bias on the edge position; our primary uncertainty in determining $`w(z)`$ to high redshift is the sparse size of the sample of distant, early-types considered here. We can now attempt to recover $`dz/dt`$ from this simulation by looking at the shift in the edge position as a function of $`z`$ (see § 3.4 and Jimenez & Loeb, 2002) and thus deduce $`w(z)`$. We recover $`w(z)=1\pm 0.2`$ for the whole redshift range. This shows that, with a sample similar to the LRG sample extended to higher redshift, we should be able to measure $`w(z)`$ within a few 10%. This is similar to what was already concluded in Jimenez & Loeb (2002).

Having established that the age–redshift relation for the old envelope is not significantly affected by the bias caused by later episodes of star formation, we now consider the observational limits that we can impose on $`w(z)`$. As a consistency check we note that the age at $`z=0`$ obtained here is in good agreement with GC ages (see §2) and that the Treu et al. galaxy ages agree with those from of SDSS LRG sample where they overlap. The solid line in Figure 4 corresponds to the age–redshift relation for our flat, reference $`\mathrm{\Lambda }`$CDM model: $`\mathrm{\Omega }_m=0.27`$, $`h=71`$, and $`w=1`$. This is consistent with the observed age–redshift envelope, and it seems to indicate that we live in a universe with a classical vacuum energy density and that, on average, stars in the oldest galaxies formed about 0.7 Gyr after the Big Bang (e.g., at redshifts $`z710`$). To make the observed age–redshift relation consistent with a model with widely different $`w(z)`$ behavior, we would have to infer that early-type galaxies at different redshifts have very different formation epochs for their stellar populations. For example, the dotted line in Figure 4 indicates the age-redshift relation for an evolving population in a universe where $`w(z)=2`$ for $`z>1`$ and grows linearly from $`w=2`$ at $`z=1`$ to $`w=0`$ at $`z=0`$. For this model to work, we would have to conclude that the oldest galaxies at $`z1.5`$ formed about 0.7 Gyr after the Big Bang, but that the oldest galaxies at $`z0`$ formed $`3`$ Gyr after the Big Bang. Since the oldest galaxies in the LRG sample are as old as the oldest GCs, we consider this to be an unlikely explanation. Furthermore, not only would this scenario require that *a*ll local galaxies either formed their stars more recently than galaxies at high redshift or had a recent episodes of recent star formation, but would also require a remarkable fine-tuning for them all to obtain exactly the same age. We therefore conclude that this extreme model for $`w(z)`$ is unlikely given the data.

Unfortunately, the small number of galaxies in our samples at $`z>0.2`$ does not allow us to compute $`dz/dt`$ with enough accuracy to constrain $`w(z)`$, and therefore this measurement will have to await better data. However, the $`z<0.2`$ region is well populated by LRG galaxies (Figure 4) and we can determine $`dt/dz`$ at $`z0`$. This is illustrated in the bottom panel of Figure 4 by the clear shift in the upper envelope of the age histogram between the two redshift ranges $`0<z<0.04`$ and $`0.08<z<0.12`$. In the next section we will investigate this phenomenon in more detail and use it to determine $`H_0`$.

### 3.4 The value of $`H_0`$

By concentrating on the age–redshift relation for SDSS LRG galaxies at $`z0.2`$, we now determine the Hubble constant $`H_0`$. This procedure relies upon determining the “edge” of the galaxy age distribution in different redshift intervals.

If the ages of the galaxies were known with infinite accuracy, for each galaxyi at redshift $`z_i`$, one could associate an age $`a_i`$. The probability for the age of the oldest stellar population in galaxyi, $`P(t|a_i)`$, would be given by a step function which jumps from 0 to 1 at $`a_i`$. An age–redshift relation edge could then be obtained by dividing the galaxy sample into suitably-large redshift bins and multiplying the $`P(t|a_i)`$ for all the galaxies in each bin.

In practice, the age of each galaxy is measured with some error $`\delta a_i`$. We thus assume $`P(t|a_i)=1`$ if $`t>a_i+\delta _{a_i}`$ and $`\mathrm{ln}P(a_i)=x^2`$ where $`x=(a_i+\delta _{a_i}t)/(\sqrt{2}\delta _{a_i})`$ otherwise. We have divided the $`z<0.2`$ portion of the LRG sample into 51 redshift bins. For each bin, we obtain $`P(z,t)`$ by multiplying the $`P(t|a_i)`$ of the galaxies in that bin. We define the “edge” of the distribution $`t(z)`$ to be where $`\mathrm{ln}P(z,t)`$ drops by $`0.5`$ from its maximum (i.e. $`\mathrm{\Delta }\mathrm{ln}P=0.5`$. We associate an error to this determination $`\delta _{t(z)}`$ given by $`t_2(z)t(z)`$ where $`t_2(z)`$ corresponds to where $`\mathrm{ln}P(z,t)`$ drops by $`2`$ from its maximum (i.e. $`\mathrm{\Delta }\mathrm{ln}P=2`$, this approximately corresponds to the 68% confidence level). This procedure makes the determination of the “edge” less sensitive to the outliers (see also Raychaudhury et al. (1997)). Figure 6 presents the resultant $`t(z)`$ relation.

For $`z<z_{\mathrm{max}}`$, we fit $`t(z)`$ with a straight line whose slope $`dt/dz`$ is related to the Hubble constant at an effective redshift by $`H(z_{\mathrm{eff}})=\frac{1}{1+z_{\mathrm{eff}}}(\frac{dt}{dz})^1`$. This fit is performed by standard $`\chi ^2`$ minimization. We also compute $`P_{\chi ^2}`$, the probability of obtaining equal or greater value of the reduced $`\chi ^2`$ if the $`t(z)`$ points were truly lying on a straight line. Values of $`P_{\chi ^2}<0.1`$ means that a straight line is not a good fit to the points.

For $`w=1`$, we obtain $`H_0=H(z_{\mathrm{eff}})\left[\mathrm{\Omega }_m(1+z_{\mathrm{eff}})^3+\mathrm{\Omega }_\mathrm{\Lambda }\right]^{\frac{1}{2}}`$ where $`\mathrm{\Omega }_m=0.27`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. Figure 7 shows how the $`H_0`$ measurement depends on $`z_{\mathrm{max}}`$. For $`z_{\mathrm{max}}>0.17`$ not only the value for $`H_0`$ drifts but also a straight line is not a good fit to the points ($`P_{\chi ^2}`$ suddenly drops below $`0.1`$). For our $`H_0`$ determination, $`z_{\mathrm{max}}=0.17`$, $`P_{\chi ^2}=0.32`$, $`z_{\mathrm{eff}}=0.09`$, and the correction from $`H(z_{\mathrm{eff}})`$ to $`H_0`$ is a $`4\%`$ effect. We obtain $`H_0=69\pm 12`$ km s-1 Mpc-1.

Our $`H_0`$ determination is in good agreement with the Hubble Key Project measurement ($`h=0.72\pm 0.03\pm 0.07`$; Freedman et al., 2001), with the valued derived from the joint likelihood analysis of WMAP \+ 2dFGRS + Lyman-$`\alpha `$ forest power spectrum ($`h=0.71_{0.08}^{+0.04}`$; Spergel et al., 2003), with gravitational lens time delay determinations ($`h=0.59_{0.07}^{+0.12}\pm 0.03`$; Treu & Koopmans, 2002), and Sunyaev-Zeldovich measurements ($`h=0.60_{0.040.18}^{+0.04+0.13}`$; Reese et al., 2002). We emphasize that our method is fully independent from any of the above measurements; it is reassuring that so many different methods yield such a consistent value for $`H_0`$.

## 4 Conclusions

We have obtained tight constraints on the effective value of $`w`$ using solely the determination of the ages of the oldest Galactic globular clusters and the observed position of the first acoustic peak of the CMB power spectrum. Our constraint on $`w`$ is in good agreement with that of Spergel et al. (2003), but is subject to different possible systematics and is very weakly model dependent. Our constraint is slightly less stringent than that obtained by Spergel et al. (2003) from a joint likelihood analysis of WMAP with six external data sets (WMAP \+ CBI + ACBAR + 2dFGRS + Lyman$`\alpha `$ forest power spectrum + Type Ia supernovae + HST Key Project $`H_0`$ determination), but is tighter than the CMB-only (WMAP+CBI+ACBAR) determination and comparable to the WMAP \+ ACBAR + CBI + HST constraint (Spergel et al., 2003). We speculate that the joint constraint from WMAP and all the external data sets considered in Spergel et al. (2003) could be further improved by adding stellar populations age determination as an additional external data set. We have restricted our analysis to values of $`w>1`$ since theoretical models for $`w<1`$ seem difficult to construct (Carroll et al., 2003). The use of external data sets (see e.g. Kujat et al. (2002) and references therein for other approaches to measure $`w`$) could help to remove this restriction and constrain $`w`$ for values $`<1`$.

We have dated a large sample of galaxies and shown that a clear age–redshift envelope exists. By computing the relative ages of SDSS red galaxies at $`z<0.2`$ we obtain an estimate of $`H_0`$ independent of and in good agreement with determinations from the Hubble Key Project, gravitational lensing time delays, Sunyaev-Zeldovich, and WMAP. We also show that when surveys with a substantial number of high $`S/N`$ spectra of old galaxies at $`z>>0`$ become available, a determination of the possibly time-dependent equation of state parameter $`w(z)`$ will be obtainable.

We thank David Spergel for discussions that stimulated this investigation. We gratefully thank Fiona Harrison and Megan Eckart for providing the time during which the spectra of MS1054$``$0321 were obtained. We also thank Peter Eisenhardt and Charles Lawrence for useful discussions and Hiranya Peiris for insightful discussion on the position of the peak and $`w`$. LV thanks David Schelegel for help with the SDSS early data release catalogue. The work of RJ is supported in part by NSF grant AST-0206031. LV is supported by NASA through Chandra Fellowship PF2-30022 issued by the Chandra X-Ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-39073. The work of DS was carried out at Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. The work of TT is partially supported by STScI grant AR-09222. Finally, the authors wish 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.