text:
# Temperature Dependence of Hall Response in Doped Antiferromagnets

(March 14, 2024)

## Abstract

Using finite-temperature Lanczos method the frequency-dependent Hall response is calculated numerically for the $`t`$-$`J`$ model on the square lattice and on ladders. At low doping, both the high-frequency $`R_\text{H}^{}`$ and the d.c. Hall coefficient $`R_\text{H}^0`$ follow qualitatively similar behavior at higher temperatures: being hole-like for $`T>T_s!!1.5J`$ and weakly electron-like for $`T<T_s`$. Consistent with experiments on cuprates, $`R_\text{H}`$ changes, in contrast to $`R_\text{H}^{}`$, again to the hole-like sign below the pseudogap temperature $`T^{}`$, revealing a strong temperature variation for $`T0`$.

The anomalous behavior of the Hall constant $`R_\text{H}`$ in the normal state of cuprates ong  remains the challenge for theoreticians for over a decade. Two aspects, possibly interrelated, are evident and should be understood: a) the d.c. $`R_\text{H}^0`$ at low temperatures $`T0`$ is clearly doping dependent. In the prototype material La2-xSrxCuO4 (LSCO) it changes from positive $`R_\text{H}^01/x`$ at low doping $`x<x^{}0.3`$, consistent with the picture of hole-doped (Mott-Hubbard) insulator, to the electron-like $`R_\text{H}^0<0`$ at $`x>x^{}`$ in agreement with the usual band picture. b) $`R_\text{H}^0`$ is also strongly temperature dependent, both at low doping and optimum doping. At optimum doping, the attention has been devoted to the anomalous variation of the Hall angle $`\theta _\text{H}T^2`$ in YBa2Cu3O7 ande . On the other hand, at low hole concentration $`c_h<0.15`$, $`R_\text{H}(T)`$ in LSCO has been shown to follow an universally behaved hwan  decrease with $`T`$ in which $`R_\text{H}^0(T0)`$ and the characteristic temperature $`T^{}`$ of vanishing $`R_\text{H}^(T^{})0`$ both scale with $`c_h`$. In underdoped cuprates, the same $`T^{}(c_h)`$ has been in fact associated with the (large) pseudogap crossover scale in uniform susceptibility $`\chi _0(T)`$, in-plane resistivity $`\rho (T)`$, specific heat $`c_v(T)`$, and some other quantities imad .

A number of theoretical investigations have addressed the first question, i.e. the doping dependence of $`R_\text{H}`$ in models of strongly correlated electrons, in particular within the $`t`$-$`J`$ model and the Hubbard model on a planar lattice. The advantage is that one can study the dynamical Hall response and the d.c. Hall constant as a ground state ($`T=0`$) property, in particular in systems with finite transverse dimension cast ; rojo  and in the ladder geometry tsun ; prel . It has been also shown that within the $`t`$-$`J`$ model the change from a hole-like to an electron-like Hall response can be qualitatively reproduced by studying the high-frequency $`R_\text{H}^{}=R_\text{H}(\omega \mathrm{})`$ shas , analytically tractable at $`T\mathrm{}`$. Recently, a connection of the reactive $`R_\text{H}^0(T=0)`$ to the charge stiffness has also been found zoto .

The anomalous temperature dependence of $`R_\text{H}(T)`$, being the main subject of this work, has been much less clarified in the literature, The Hall mobility $`\mu _\text{H}(T)`$ of a single charge carrier in the Mott-Hubbard insulator has been first evaluated within the generalized retraceable path approximation brin . The high-frequency $`R_\text{H}^{}(T)`$ has been calculated using the high-$`T`$ expansion shas . At low doping, $`c_h<0.15`$, it has been observed that on decreasing temperature $`R_\text{H}^{}`$ is also decreasing instead of approaching presumed (larger) semiclassical and experimentally observed d.c. result $`R_\text{H}^c=1/c_he_04R_\text{H}^{}(T=\mathrm{})`$. Related are the conclusions of the quantum Monte-Carlo study of the planar Hubbard model assa , where close to the half-filling electron-like $`R_\text{H}^{}<0`$ has been found at low $`T`$. The same has been claimed generally for $`R_\text{H}(\omega )`$ even for low $`\omega `$ assa . Quite controversial are also results for $`R_\text{H}^0(T)`$ on ladders tsun . In regard to that, we should also mention the questionable relation of the off-diagonal $`\sigma _{xy}`$ to the orbital susceptibility $`\chi _d`$ wido ; rojo , potentially useful as an alternative route to the understanding of $`R_\text{H}^0(T)`$ vebe .

In the following we present numerical results for the dynamical $`R_\text{H}(\omega )`$, as obtained within the low doping regime of the $`t`$-$`J`$ model using the finite-temperature Lanczos method (FTLM) jplan ; jprev . The aim of this letter is to approach the low-$`\omega `$ and low-$`T`$ limit as much as possible and to investigate the relation between $`R_\text{H}^{}(T)`$ and $`R_\text{H}^0(T)`$. We find these two quantities essentially different for $`T<T^{}`$, establishing the pseudogap scale $`T^{}<J`$ both in the ladder and planar systems.

We study the $`t`$-$`J`$ model in an external homogeneous magnetic field $`𝐁=\mathrm{curl}𝐀`$,

$`H(𝐀)=`$ $`t{\displaystyle \underset{ijs}{}}(\mathrm{e}^{\mathrm{i}\theta _{ij}}\stackrel{~}{c}_{is}^{}\stackrel{~}{c}_{js}^{}+\text{H.c.})+`$
$`+J{\displaystyle \underset{ij}{}}(𝐒_i𝐒_j\frac{1}{4}n_in_j),`$ (1)
where the (inhomogeneous) vector potential enters the phases $`\theta _{ij}=e𝐀(𝐫_i)𝐫_{ij}`$. The hopping is only between the nearest neighbors $`ij`$. Projected fermionic operators $`\stackrel{~}{c}_{is}^{}`$, $`\stackrel{~}{c}_{is}^{}`$ do not allow for the double occupancy of sites.

In order to calculate the dynamical Hall coefficient

$$R_\text{H}(\omega )=\frac{\rho _{xy}(\omega )}{B}|_{B0}=\frac{\sigma _{xy}(\omega )}{B\sigma _{xx}(\omega )\sigma _{yy}(\omega )}|_{B0},$$
(2)
the conductivity tensor is evaluated within the linear response theory,

$$\sigma _{\alpha \beta }(\omega )=\frac{\mathrm{i}e^2}{N\omega ^+}\left[\tau _{\alpha \beta }\mathrm{i}_0^{\mathrm{}}\text{d}t\mathrm{e}^{\mathrm{i}\omega t}[j_\alpha (t),j_\beta ]\right],$$
(3)
where in the presence of $`B0`$ the particle current $`𝐣`$ and the stress tensor $`\underset{¯}{\tau }`$ operators are given by

$`𝐣`$ $`=t{\displaystyle \underset{ijs}{}}𝐫_{ij}(\mathrm{ie}^{\mathrm{i}\theta _{ij}}\stackrel{~}{c}_{is}^{}\stackrel{~}{c}_{js}^{}+\text{H.c.}),`$
$`\underset{¯}{\tau }`$ $`=t{\displaystyle \underset{ijs}{}}𝐫_{ij}𝐫_{ij}(\mathrm{e}^{\mathrm{i}\theta _{ij}}\stackrel{~}{c}_{is}^{}\stackrel{~}{c}_{js}^{}+\text{H.c.}).`$ (4)
On a square lattice with $`N`$ sites and periodic boundary conditions (b.c.) one cannot apply arbitrary magnetic field $`B`$ since only quantized $`B=B_m=m\mathrm{\Phi }_0a^2/N`$ can be made compatible with the periodic b.c. assa . Therefore the smallest but finite $`B=B_1`$ is used in calculations. The square lattices used are in general Euclidean (tilted) $`N=l^2+n^2`$, in particular we investigate systems $`N=10`$, 16, 18. On the other hand, the ladder geometry of $`N=L\times M`$ sites with the periodic b.c. in the $`L`$ direction and open b.c. in the perpendicular $`M`$ direction allows for any finite $`B0`$, the fact already used in several $`T=0`$ calculations cast ; rojo ; prel . The advantage of ladder systems is also the existence of the reference ground-state results $`R_\text{H}^0(T=0)`$ which seem to be better understood prel ; zoto . Furthermore, at low doping they reproduce the simple semiclassical behavior $`R_\text{H}^0(T=0)R_\text{H}^\text{c}=1/c_he_0`$.

Dynamical components $`\sigma _{\alpha \beta }(\omega )`$ are evaluated using the FTLM jplan , employed so far for various dynamic and static quantities within the $`t`$-$`J`$ model jprev , among them also the $`B=0`$ optical conductivity $`\sigma (\omega )=\sigma _{\alpha \alpha }(\omega )`$ on a square lattice. Comparing to the diagonal $`\sigma _{\alpha \alpha }`$, the evaluation of the off-diagonal $`\sigma _{xy}(\omega )`$ is more demanding for several reasons: a) the introduction of $`B>0`$ in the model (1) breaks the translational invariance and prevents the reduction of the basis states in the Lanczos procedure, hence available finite-size systems are somewhat smaller, b) we expect $`\sigma _{xy}(\omega )B`$ while $`\sigma _{xy}(B=0,\omega )`$ does not vanish identically within the FTLM; consequently larger sampling over initial wavefunctions jplan ; jprev  are needed to reduce the statistical error, c) on a finite square lattice the reference result $`R_\text{H}^0(T=0)`$ is not meaningful for $`B_m>0`$, while in ladder systems it is quite sensitive to the introduction of an additional flux prel . Nevertheless, in general, restrictions for the validity of the FTLM results are similar to other quantities. Through the thermodynamic partition function $`Z(T_{\text{fs}})=Z^{}`$, we can define the marginal finite-size $`T_{\text{fs}}`$ below which too few levels contribute to the average and results loose the thermodynamic validity jprev . In the following, we analyze results for $`J=0.4t`$ at low hole doping $`c_h=N_h/N`$ ($`N_h=1`$, 2). In this regime we can estimate $`T_{\text{fs}}/t0.150.20.5J/t`$.

Let us first present results for the dynamical $`R_\text{H}(\omega )`$. In Fig. 1 we show the normalized real part $`r_\text{H}=e_0c_h\mathrm{Re}R_\text{H}`$ for systems with a single hole $`N_h=1`$. In the evaluation of $`R_\text{H}(\omega )`$ from Eq. (2) we insert complex $`\sigma _{\alpha \alpha }`$ at $`B=0`$ and the most sensitive quantity remains $`\sigma _{xy}(\omega )`$ calculated at $`B=B_1`$ on a square lattice and $`B0.3B_1`$ on ladders. In the presentation of results an additional frequency smoothening $`\delta =0.2t`$ is used. The normalization of $`R_\text{H}`$ is chosen such that at low doping $`r_\text{H}=1`$ would show up in the case of the semiclassical result.

In Fig. 1 several common features of $`R_\text{H}`$ in the ladder geometry and in the 2D systems are recognized:

a) $`r_\text{H}(\omega )`$ is quite smoothly varying function of $`\omega `$, at least in contrast to strongly $`\omega `$-dependent $`\mathrm{Re}\sigma (\omega )`$ on a 2D system, which is found jprev  to decay with an anomalous relaxation rate $`1/\tau (\omega )\omega +\xi T`$.

b) At high temperatures $`T>t`$ we get a hole-like $`r_\text{H}>0`$ for all systems. In this regime $`r_\text{H}(\omega )`$ is very smooth, in particular for the $`M=4`$ ladder and the 2D lattice.

c) For low temperatures $`T<t`$, $`r_\text{H}(\omega )`$ is less smooth and the dependence is more pronounced for the 2-leg ladder. On the other hand, $`M=4`$ ladder clearly approaches the behavior of the 2D system, whereby both of the latter show quite a modest variation of $`r_\text{H}(\omega )`$. In all systems the resonances (and the variation) visible in $`r_\text{H}(\omega )`$ at high $`\omega >t`$ reflect the predominantly local physics of the hole motion and are thus not related to a current relaxation rate deduced from $`\sigma (\omega )`$.

Results for $`r_\text{H}(\omega )`$ are the basis for the calculation of high-frequency $`r_\text{H}^{}=r_\text{H}(\omega =\mathrm{})`$ as well as the d.c. limit $`r_\text{H}^0=r_\text{H}(\omega 0)`$. The latter is more sensitive since in a finite system (even at $`T>0`$) $`\sigma _{\alpha \beta }(\omega 0)`$ can be singular due to the coherent charge transport in a system with periodic b.c.. The coherent transport shows up in a finite (but small) charge stiffness jprev , which should be omitted in the evaluation of Eq. (2). In any case, one should take into account proper $`\omega 0`$ behavior of dissipative systems at $`T>0`$ which is different in ladders and in 2D lattices, respectively: a) On a ladder we get in the leading order of $`\omega 0`$ a normal conductance along the $`x`$-direction, i.e. $`\sigma _{xx}(\omega 0)\sigma _0`$, but a finite polarizability along the $`y`$-direction, $`\sigma _{yy}(\omega 0)\omega \chi _{yy}^0`$. Hence, we expect $`\sigma _{xy}\omega `$ and finite $`r_\text{H}^0`$. b) For a macroscopic isotropic 2D system we get $`\sigma _{\alpha \alpha }(\omega 0)\sigma _0`$ and we expect as well $`\sigma _{xy}\sigma _{xy}^0`$, leading to finite $`r_\text{H}^0`$.

In Fig. 2 we present results for $`r_\text{H}^{}(T)`$ and $`r_\text{H}^0(T)`$ for the ladder systems with $`N_h=1`$. Results are shown for 2-leg ladders with various lengths $`L=9`$, 10, 11 and for $`M=4`$ ladder with $`L=4`$. Since $`r_\text{H}^0`$ and $`r_\text{H}^{}`$ are properly scaled, for given $`M`$ curves are expected to approach a well defined macroscopic limit at $`L\mathrm{}`$. In fact, $`r_\text{H}^{}`$ are nearly independent of $`L`$ (as well as of $`M`$) down to $`TT_{\text{fs}}`$. A crossover at $`T_s0.6t`$ from a hole-like $`r_\text{H}^{}>0`$ into a electron-like $`r_\text{H}^{}<0`$ can be explicitly observed. $`r_\text{H}^0`$ results are more size ($`L`$) dependent, nevertheless they reveal a crossover nearly at the same $`TT_s`$. In contrast to $`r_\text{H}^{}`$ which remains negative for the whole regime $`T<T_s`$, $`r_\text{H}^0`$ changes sign again at $`T=T^{}0.2t`$. Although our data for $`T^{}`$ are more scattered the crossover into the hole-like $`r_\text{H}^0(T<T^{})>0`$ is expected. Namely, from the ground state calculations in same systems prel  we know that $`r_\text{H}^0(T=0)1.5`$ and $`r_\text{H}^0(T=0)1.2`$ for $`M=2`$ and $`M=4`$ ladders, respectively. Therefore, it is not surprising that the observed dependence $`r_\text{H}^0(T<T^{})`$ is very steep.

Corresponding results for the planar lattice in Fig. 3 are both qualitatively and quantitatively similar. Note, that at low doping the limiting value $`r_\text{H}^{}(T\mathrm{})=1/4`$ agrees with the analytical result shas , while obtained $`r_\text{H}^0(T\mathrm{})0.3`$ is also quite close. Again, the crossover into an electron-like regime appears at $`T_s0.6t`$. For larger sizes $`N16`$ the lower crossover $`T^{}0.2t`$ is visible as well. In finite 2D systems a reference numerical result at $`T=0`$ does not exist, however, the analytical theory prel1  indicates that in a macroscopic limit with a single hole ($`N_h=1`$) in an ordered antiferromagnet one should get $`r_\text{H}^0=1`$.

In numerically available systems, $`N_h=2`$ represents already a substantial doping. Therefore, results for $`r_\text{H}^0`$ and $`r_\text{H}^{}`$ shown in Fig. 3 should be interpreted in relation with the corresponding finite doping $`c_h`$. Main message of Fig. 3 is that upper crossover $`T_s`$, still nearly the same in both $`r_\text{H}^0(T)`$ and $`r_\text{H}^{}(T)`$, shifts down quite systematically with increasing $`c_h`$, i.e. with decreasing size $`N`$ at given $`N_h`$. At least in ladder systems at $`c_h<0.3`$, we still find $`r_\text{H}^0(T=0)>0`$ in the ground state prel , therefore also the lower crossover $`T^{}<T_s`$ is expected. However, we cannot detect such a crossover in $`r_\text{H}^0(T)`$ down to $`T_{\text{fs}}0.15t`$, not surprisingly since also the experimental value, e.g. in LSCO at $`c_h>0.1`$, is $`T^{}<600\text{K}0.15t`$ (assuming $`t0.4e`$V).

Let us finally comment on the relation of the d.c. $`\sigma _{xy}^0`$ to the orbital susceptibility $`\chi _d`$ in a macroscopic 2D system. Namely, $`\stackrel{~}{\sigma }_{xy}^0=eB\chi _d/\mu =eB(\chi _d/c_h)(c_h/\mu )`$, (where $`\mu `$ denotes chemical potential) was derived using seemingly quite general thermodynamic relations wido ; rojo , but at the same time put under question rojo . Since the d.c. $`\sigma _{\alpha \alpha }^0(T)>0`$ is quite a smooth function the above relation seems to yield also a qualitative connection between $`\chi _d(T)`$ and $`R_\text{H}^0(T)`$. The situation should be particularly simple at low doping (but not too low $`T`$), where $`c_h/\mu c_h/T`$ and $`\chi _dc_h`$ is expected, and consequently $`\stackrel{~}{\sigma }_{xy}B\chi _d/T`$. Indeed, results for $`N_h=1`$ indicate vebe  that both crossovers $`T_s`$ and $`T^{}`$ appear also as a change of sign in $`\chi _d(T)`$ nearly at the same values. Here, the intermediate regime $`T^{}<T<T_s`$ corresponds to an anomalous paramagnetic response $`\chi _d>0`$. On the other hand, it is quite evident from our results that the relation is not valid at high $`Tt`$. Namely, in this regime $`\sigma _{\alpha \alpha }^01/T`$ and $`\sigma _{xy}^0B/T^2`$ brin  is obtained, leading to $`R_\text{H}^0(T\mathrm{})\text{const}`$. On the other hand, from the high-$`T`$ expansion $`\chi _d1/T^3`$ is acquired vebe , so that the assumed relation would demand $`\stackrel{~}{\sigma }_{xy}^0B/T^4`$, in conflict with previous $`\sigma _{xy}^0B/T^2`$.

In conclusion, we have presented results for both dynamical and d.c. Hall constant within the $`t`$-$`J`$ model on ladders and on square lattices. The main novel point is the observation of two crossover temperatures $`T_s`$ and $`T^{}`$ which are at low doping generally present in all systems. Both $`R_\text{H}^{}`$ and $`R_\text{H}^0`$ are positive at $`T>T_s`$ and change sign at $`T_s`$. While $`R_\text{H}^{}(T<T_s)`$ stays negative, $`R_\text{H}^0`$ reveals a sign change into a hole-like behavior at $`T=T^{}<T_s`$ as well as steep variation of $`R_\text{H}^0(T<T^{})`$. This reconciles some seemingly controversial theoretical results shas ; assa . Our results are in agreement with high-$`T`$ expansion results for $`R_\text{H}^{}(T)`$ which at low $`c_h`$ also show decreasing positive values with decreasing $`T`$. Quantum Monte Carlo results within the Hubbard model for $`R_\text{H}(\mathrm{i}\omega )`$ correspond effectively to high (imaginary) frequencies and low $`T`$, and being negative they are in agreement with our findings for $`R_\text{H}^{}`$.

How should we understand the above numerical results? At high $`Tt`$ and low doping $`c_h1`$, $`R_\text{H}^{}`$ as well as $`R_\text{H}^0`$ are governed by a loop motion (that is where the dependence on $`B0`$ comes from) of a hole within a single plaquette brin ; shas . One expects $`R_\text{H}^{}>0`$, but $`r_\text{H}^{}=1/4`$ is a non-universal value which e.g. depends on the lattice coordination shas . The electron-like $`r_\text{H}^{}(T=0)<0`$ represents an instantaneous Hall response within the ground state near half filling is harder to explain, but is clearly the signature of strong correlations. On the other hand at low $`T`$, $`R_\text{H}^0`$ tests the (low energy) quasiparticle properties. Evidently, at low doping and $`T<T^{}`$ at least a single hole in an antiferromagnetic spin background behaves as a well defined hole-like quasiparticle leading to $`r_\text{H}^0(T0)1`$ both in 2D prel1  and in ladders prel . Our results for vanishing $`R_\text{H}^0(T^{})0`$ indicate that the quasiparticle character is essentially lost at quite low $`TT^{}<J`$, with the pseudogap scale $`T^{}(c_h)`$ decreasing with doping. Such phenomenon is possibly consistent with the scenario of electrons being effectively composite particles (spinons and holons) in strongly correlated systems ande ; laug , at least at $`T>T^{}(c_h)`$, whereby $`T^{}(c_h)`$ vanishes at optimum doping.

Finally, let us note that our results for $`R_\text{H}^0`$ are in several aspects consistent with experiments on cuprates, and with LSCO in particular. At low doping $`c_h<0.1`$ we find $`T^{}J/2`$, close to the observed $`T^{}800`$K. At the same time, we find a very steep dependence in $`R_\text{H}^0(T<T^{})`$. With increasing $`c_h`$, $`T^{}(c_h)`$ seems to have desired decreasing tendency, although to establish this beyond a reasonable doubt more work is needed.

Authors acknowledge the support of the Ministry of Education, Science and Sport of Slovenia.


format: text

text:
# Evolution of a coherent array of Bose-Einstein Condensates in a magnetic trap

## I Introduction

The development of the technologies of the laser trapping and evaporative cooling has yielded intriguing Bose-Einstein condensates (BECs) , a state of matter in which many atoms are in the same quantum mechanical state. The realization of BECs has made remarkable theoretical and experimental advances on this exotic quantum system . Recently, optical lattices are used to investigate further the unique character of the ultra-cold atoms . The application of the optical lattices to the ultra-cold atoms is very promising, such as the quantum computing scheme proposed in . Experimentally, the macroscopic quantum interference effect  and thermodynamic properties  of the BECs in the optical lattices have been investigated thoroughly. In addition, the superfluid and dissipative dynamics  of a BEC in the optical lattices are investigated in the experiments too. In particular, through the application of the optical lattices, the quantum phase transition from a superfluid to a Mott insulator in Bose-Einstein condensed gases is observed in a recent experiment .

Recently, the expansion of a coherent array of BECs is carried out to illustrate the interference effect in the experiment by Pedri et al. , when both the magnetic trap and optical lattices are switched off. In another experiment by Morsch et al. , the expansion of the condensates is investigated when only the magnetic trap is switched off. In this paper, we give a theoretical investigation on the expansion of the condensates when only the optical lattices are switched off. Due to the confinement of the harmonic potential, we find that there are several interested phenomena not encountered in the experiment by Pedri et al.  and Morsch et al. . In the presence of the harmonic potential, the interference effect would be much stronger than the case when the magnetic trap is switched off too. In this situation, researches show that the interaction between atoms would give important correction, in contrast to the case investigated by Pedri et al. .

## II Wave function in harmonic potential after optical lattices are switched off

In the experiment conducted by Pedri et al. , the external potential of the Bose gas is given by :

$$V=\frac{1}{2}m\left(\omega _x^2x^2+\omega _{}^2\left(y^2+z^2\right)\right)+sE_R\mathrm{cos}^2\left(\frac{2\pi x}{\lambda }+\frac{\pi }{2}\right),$$
(1)
The last term represents the external potential due to the presence of the optical lattices. In the above expression, $`\omega _x`$ and $`\omega _{}`$ are the axial and radial frequencies of the harmonic potential, respectively. In addition, $`\lambda `$ is the wavelength of the retroreflected laser beam, and $`sE_R`$ denotes the depth of the optical lattices. For the optical lattices created by the retroreflected laser beam, the last term has a period $`d=\lambda /2`$. In other words, $`d`$ can be regarded as the distance between two neighboring wells induced by the optical lattices. In this paper, the experimental parameters in  are used to calculate various physical quantities, thus it would be useful to give them here. The experimental parameters in  are $`\omega _x=2\pi \times 9`$ $`\mathrm{Hz}`$, $`\omega _{}=2\pi \times 92`$ $`\mathrm{Hz}`$, $`\lambda =795`$ $`\mathrm{nm}`$, $`E_R=2\pi \mathrm{}\times 3.6`$ $`\mathrm{kHz}`$, and $`s=5`$.

Due to the presence of the optical lattices, there are an array of condensates formed in the combined potential, when the temperature is lower than the critical temperature. In this work, we shall investigate the case of the strong tunneling between neighboring BECs, which holds in the experiment by Pedri et al. . In this situation, all the condensates are fully coherent, and can be described by a single order parameter. To emphasis the role of the optical lattices, our researches are carried out mainly on the character of the coherent array of condensates in the $`x`$direction.

In the presence of a magnetic trap, the number of atoms in each well should be different. Based on the analysis of the 3D model of the condensates in the combined potential , the ratio between the number of condensed atoms in $`k`$th and central wells is given by $`N_k/N_0=\left(1k^2/k_M^2\right)^2`$, where $`2k_M+1`$  represents the total number of condensates induced by the optical lattices. In this situation, using the Gaussian approximation in the $`x`$direction for each well, the normalized wave function in coordinate space takes the form

$$\phi _0\left(x\right)=A_n\underset{k=k_M}{\overset{k_M}{}}\left(1\frac{k^2}{k_M^2}\right)\mathrm{exp}[\left(xkd\right)^2/2\sigma ^2],$$
(2)
where $`\sigma `$ denotes the width of the condensate in each well. It can be calculated by numerical minimization of the energy of the condensates . In the following calculations in this paper, $`\sigma =0.25d`$  is used for $`s=5`$. In the above expression, the normalized constant $`A_n`$ takes the form

$$A_n=\frac{1}{\sqrt{\left(16k_M^41\right)/15k_M^3}\pi ^{1/4}\sigma ^{1/2}}.$$
(3)
It is well known that once the wave function at an initial time is known, the wave function at a later time can be obtained through the following integral equation :

$$\phi _0(x,t)=_{\mathrm{}}^{\mathrm{}}K(x,t;y,t=0)\phi _0\left(y,t=0\right)𝑑y,$$
(4)
where $`\phi _0\left(y,t=0\right)`$ is the wave function at the initial time $`t=0`$ which is given by Eq. (2), and $`K(x,t;y,t=0)`$ is the well known propagator. For the atoms in the harmonic potential, the propagator can be expressed as :

$$K(x,t;y,t=0)=\left[\frac{m\omega _x}{2\pi i\mathrm{}\mathrm{sin}\omega _xt}\right]^{1/2}\mathrm{exp}\left\{\frac{im\omega _x}{2\mathrm{}\mathrm{sin}\omega _xt}\left[\left(x^2+y^2\right)\mathrm{cos}\omega _xt2xy\right]\right\}.$$
(5)
From the formulas (2), (4) and (5), after a straightforward calculation, one obtains the analytical result of the wave function confined in the magnetic trap

$`\phi _0(x,t)=A_n\sqrt{{\displaystyle \frac{1}{\mathrm{sin}\omega _xt\left(\mathrm{ctg}\omega _xt+i\gamma \right)}}}{\displaystyle \underset{k=k_M}{\overset{k_M}{}}}(1{\displaystyle \frac{k^2}{k_M^2}})\mathrm{exp}[{\displaystyle \frac{\left(kd\mathrm{cos}\omega _xtx\right)^2}{2\sigma ^2\mathrm{sin}^2\omega _xt\left(\mathrm{ctg}^2\omega _xt+\gamma ^2\right)}}]\times `$

$$\mathrm{exp}\left[\frac{i\left(kd\mathrm{cos}\omega _xtx\right)^2\mathrm{ctg}\omega _xt}{2\gamma \sigma ^2\mathrm{sin}^2\omega _xt\left(\mathrm{ctg}^2\omega _xt+\gamma ^2\right)}\right]\mathrm{exp}\left[\frac{i\left(x^2\mathrm{cos}\omega _xt+k^2d^2\mathrm{cos}\omega _xt2xkd\right)}{2\gamma \sigma ^2\mathrm{sin}\omega _xt}\right],$$
(6)
where we have introduced a dimensionless parameter $`\gamma =\mathrm{}/m\omega _x\sigma ^2`$. Assuming $`N`$ denotes the total number of particles in the condensates, the density distribution in $`x`$direction is $`n(x,t)=N\left|\phi _0(x,t)\right|^2`$.

## III Periodicity of the Density Distribution and the Motion of $`n=\pm 1`$ Peak

Due to the confinement of the harmonic potential, the density distribution in $`x`$direction should exhibit a periodic character. From Eq. (6), it is easy to find that the period of the density distribution $`n(x,t)`$ is determined by $`\omega _xT=\pi `$. For the experiment in , $`\omega _x=2\pi \times 9`$ $`\mathrm{Hz}`$, this means that the period of the density distribution is given by $`500/9`$ $`\mathrm{ms}`$.

In Fig. 1, displayed is $`n\left(x=0,t\right)`$ when only the optical lattices are switched off. The periodicity of the density is clearly shown in the figure and in agreement with the analytical result given by $`T=\pi /\omega _x`$. We see that the density at $`x=0`$ reaches a maximum value at time $`t_m=\left(2m1\right)\pi /2\omega _x`$, with $`m`$ positive integer. At time $`t_m`$, the wave function takes the form

$$\phi _0(x,t_m)=A_n\sqrt{\frac{1}{i\gamma }}\underset{k=k_M}{\overset{k_M}{}}\left(1\frac{k^2}{k_M^2}\right)\mathrm{exp}\left[\frac{x^2}{2\sigma ^2\gamma ^2}i\frac{xkd}{\gamma \sigma ^2}\right].$$
(7)
The maximum density at $`x=0`$ is then

$$n\left(x=0,t_m\right)=\frac{NA_n^2}{\gamma }\left[\underset{k=k_M}{\overset{k_M}{}}\left(1\frac{k^2}{k_M^2}\right)\right]^2.$$
(8)
In the case of $`k_M>>1`$, the above expression can be approximated as:

$$n\left(x=0,t_m\right)N\alpha _{xideal}^2,$$
(9)
where $`\alpha _{xideal}^2`$ is given by

$$\alpha _{xideal}^2=\frac{5k_Mm\omega _x\sigma }{3\pi ^{1/2}\mathrm{}}.$$
(10)
In the experiments, the depth of the optical lattices can be changed through the variation of the parameter $`s`$. From , we know $`\sigma 1/s^{1/4}`$ and $`k_M1/\sigma ^{1/5}`$. Therefore, $`n\left(x=0,t_m\right)1/s^{1/5}`$. This shows that in the non-interacting model, the maximum density at $`x=0`$ decreases with the increase of the depth of the optical lattices.

In Figures 2 (a)-(d), we show the evolution of the density distribution of the condensates confined in the magnetic trap, after the optical lattices are switched off. The density distributions are shown at $`t=0`$, $`0.1\pi /\omega _x`$, $`0.3\pi /\omega _x`$, and $`0.5\pi /\omega _x`$ in these figures. The motion of the $`n=\pm 1`$ peaks is clearly shown in these figures. The oscillation motion of the $`n=\pm 1`$ peaks is due to the confinement of the harmonic potential. In fact, the motion of the $`n=\pm 1`$ peaks can be described very well using the classical harmonic motion. Using the classical harmonic motion, the motion of the $`n=\pm 1`$ peaks is determined by the following expression:

$$x_{n=\pm 1}\left(t\right)=\pm \frac{2\pi \mathrm{}}{m\omega _xd}\mathrm{cos}\left(\omega _xt\frac{\pi }{2}\right).$$
(11)
When obtaining the above formula, we have used the fact that the momentum distribution is characterized by sharp peaks at the values $`p_x=n2\pi \mathrm{}/d`$ . The solid line in Fig. 3 shows the harmonic motion of the $`n=1`$ peak using the above formula, while the circles show the result given by Eq. (6). We see that the classical harmonic motion agrees quite well with the result given by Eq. (6). In a sense, Eq. (11) describes the motion of the center of mass of the $`n=\pm 1`$ peak. Thus, we anticipate that the interaction between atoms will not change the motion of the $`n=\pm 1`$, although it will affect the density and width of the $`n=\pm 1`$ peak.

From Fig. 2(d) we see that at time $`t_m`$ the density distribution in $`x`$direction exhibits a very sharp peak at the center of the magnetic trap. The maximum density of the central peak is given by Eq. (9). As a comparison, assume there are $`N`$ atoms confined in an identical magnetic trap, but there are no optical lattices to induce the interference effect. In this situation, in $`x`$direction, the density distribution at $`x=0`$ is given by

$$n_{mag}\left(x=0\right)=\frac{\pi \mu _{mag}^2}{gm\omega _{}^2},$$
(12)
where $`\mu _{mag}`$ is the chemical potential of the Bose gas confined in the magnetic trap. The ratio between $`n\left(x=0,t_m\right)`$ and $`n_{mag}\left(x=0\right)`$ is then

$$\frac{n\left(x=0,t_m\right)}{n_{mag}\left(x=0\right)}=\frac{Ngm\omega _{}^2\alpha _{xideal}^2}{\pi \mu _{mag}^2}.$$
(13)
For the experimental parameters in , $`n\left(x=0,t_m\right)/n_{mag}\left(x=0\right)=29.6`$. This shows clearly that there is a very strong interference effect for the case considered here.

If both the magnetic trap and optical lattices are switched off, there is also a sort of interference effect. The maximum value of the density distribution in this case is given by

$$n_{bs}\left(x=0\right)=\frac{N_0}{\sqrt{\pi }\sigma }.$$
(14)
From Eq. (14), $`n_{bs}\left(x=0\right)s^{1/4}`$. In contrast to the case where only the optical lattices are switched off, $`n_{bs}\left(x=0\right)`$ decreases with the decrease of the depth of the optical lattices.

The ratio between $`n_{bs}\left(x=0\right)`$ and $`n_{mag}\left(x=0\right)`$ is then

$$\frac{n_{bs}\left(x=0\right)}{n_{mag}\left(x=0\right)}=\frac{N_0gm\omega _{}^2}{\pi ^{3/2}\sigma \mu _{mag}^2}.$$
(15)
For the experimental parameters in , we have $`n_{bs}\left(x=0\right)/n_{mag}\left(x=0\right)=2.2`$. From Eqs. (9) and (14), $`n\left(x=0,t_m\right)/n_{bs}\left(x=0\right)=13.7`$. Therefore, when only the optical lattices are switched off, the interference effect would be much stronger than the case when the magnetic trap is switched off too. When only the optical lattices are switched off, at $`t_m`$ the density of the central peak is very high, and we anticipate that in this situation the interaction between atoms would give important correction.

## IV Decay and revival of the density oscillation

From Fig. 1, we see that there is a phenomenon of decay and revival of the density oscillation at $`x=0`$. We now turn to discuss this unique character. To proceed, it is useful to introduce an important time scale which determines when the interference between two neighboring condensates begins to occur. Before the magnetic trap and optical lattices are switched off, from the Gaussian approximation of the condensates in each well, the width of the condensates in each well is given by

$$\mathrm{\Delta }x_0^2=\frac{x^2\mathrm{exp}\left[x^2/\sigma ^2\right]𝑑x}{\mathrm{exp}\left[x^2/\sigma ^2\right]𝑑x}.$$
(16)
It is easy to get $`\mathrm{\Delta }x_0=\sigma /\sqrt{2}`$ from the above formula. When the optical lattices are switched off, for time much smaller than $`\pi /\omega _x`$, the condensates can be approximated as a free expansion and the width of the condensates would increase in this situation. Based on the analysis of the spreading of the wave packet, the width of each condensate is given by

$$\mathrm{\Delta }x\left(t\right)=\mathrm{\Delta }x_0\sqrt{1+\frac{\mathrm{}^2t^2}{m^2\mathrm{\Delta }x_0^4}}.$$
(17)
When $`\mathrm{\Delta }x\left(t\right)=d`$, the condensates in neighboring wells begin to interfere with each other. By setting $`\mathrm{\Delta }x\left(t\right)=d`$ in the above formula, we obtains a time scale $`t_w`$ which determines when the interference between neighboring condensates begins to occur. From Eq. (17), it is easy to find that the following analytical result of $`t_w`$ can give a rather well approximation

$$t_w=\frac{\sigma dm}{\sqrt{2}\mathrm{}}.$$
(18)
From Eqs. (17) and (18), for $`t>t_w`$, one gets the following useful result:

$$\mathrm{\Delta }x\left(t\right)=\frac{t}{t_w}d.$$
(19)
As illustrated in Fig. 1, the oscillation of the density at $`x=0`$ will cease ultimately when $`t>k_Mt_w`$. However, when the time approaches $`\pi /\omega _x`$, the density oscillation will reappear. Note that the time for the revival of the density oscillation is determined solely by the axial frequency of the harmonic potential. This shows that the confinement of the harmonic potential plays a crucial role for the revival phenomenon of the density oscillation. To verify further the decay and revival of the density oscillation, Fig. 4 shows the density at $`x=k_Md/2`$. Analogous decay and revival of the density oscillation are illustrated clearly in the figure. However, the oscillation of the density disappears at a longer time $`t=1.47k_Mt_w`$, in comparison with the case at $`x=0`$.

In Fig. 5, we display the time of the disappearance of the density oscillation for different locations in the region $`0<t<\pi /\omega _x`$. We can give a rather simple interpretation for this result. When the optical lattices are switched off, the width of the expanding condensates in each well will increase. For the location at $`x=0`$, there are more and more expanding BECs interfere at this point with the development of the time. This is the reason why the density at the point $`x=0`$ will oscillate intensely. When $`t>k_Mt_w`$, however, all expanding BECs have participated in the interference at the point $`x=0`$. Therefore, the oscillation of the density at $`x=0`$ will cease at time longer than $`k_Mt_w`$. Generalizing this result, the time for the disappearance of the density oscillation for different locations is given by the following simple expression:

$$t=(k_M+x/d)t_w.$$
(20)
The solid line in Fig. 5 displays the above analytical result. We see from Fig. 5 that this simple expression agrees well with the result given by Eq. (6). Maybe the slight difference from the result given by Eq. (6) lies in the fact that we do not account for the effect of non-uniform atom distribution in each well when obtaining Eq. (20).

It is worth pointing out that the disappearance of the density oscillation do not mean the disappearance of the interference effect between the expanding condensates. When $`0<t<<\left(2k_M+1\right)t_w`$, there are only several expanding condensates interfering with each other, and there are interference fringes (or density oscillation) in this situation. For $`\left(2k_M+1\right)t_w<t<\pi /\omega _x\left(2k_M+1\right)t_w`$, however, all expanding condensates will interfere with each other, and this means the emergence of the diffraction fringes. In fact, $`n=0`$ and $`n=\pm 1`$ peaks in Fig. 2 (a)-(d) should be regarded as the diffraction fringes, rather than the interference fringes. Note that the phenomena of diffraction and interference are basically equivalent. Different from the interference phenomenon, however, the diffraction phenomenon should be regarded as a consequence of interference from many coherent wave sources. In a sense, Eq. (20) gives the time for the emergence of the diffraction fringes, which means the disappearance of the density oscillation.

## V Interaction correction to the central peak at $`t_m`$

In the case of non-interacting model, we have shown that, at time $`t_m`$, there would be a sharp central peak in the magnetic trap. In this case, the interaction between atoms can not be simply omitted, in contrast to the case when the magnetic trap is switched off too. At time $`t_m`$, the central peak of the density in $`x`$direction can be approximated as a Gaussian distribution. After the optical lattices are switched off, using Thomas-Fermi approximation in the radial direction, the square of the modulus of the 3D wave function at time $`t_m`$ takes the form

$$\left|\phi _0(x,r_{},t_m)\right|^2=\alpha _x^2\alpha _{}^2\mathrm{exp}\left[\frac{2x^2}{R_x^2}\right]\left(1\frac{r_{}^2}{R_{}^2}\right),$$
(21)
where $`R_{}^2=\sqrt{2\mu _0/m\omega _{}^2}`$ with $`\mu _0=m\omega _x^2k_M^2d^2/2`$ . In the above expression, $`\alpha _x=(2/\pi )^{1/4}/\sqrt{R_x}`$ and $`\alpha _{}=\sqrt{2/\pi R_{}^2}`$ are the normalized constants in the axial and radial directions, respectively. Obviously, $`N\alpha _x^2`$ represents the density $`n\left(x=0,t_m\right)`$ in the $`x`$direction. Due to the repulsive interaction between atoms, we anticipate that $`\alpha _x^2<\alpha _{xideal}^2`$.

Assume $`E_{int}`$ and $`E_{kin}`$ are the interaction energy and kinetic energy of the central peak at $`t_m`$, respectively. The interaction energy of the central peak is given by

$$E_{int}=\frac{gN^2}{2}\left|\phi _0(x,r_{},t_m)\right|^4𝑑V=\frac{\sqrt{2}gN^2\alpha _x^2}{3\pi R_{}^2}\text{.}$$
(22)
Assuming the total energy of the condensates is $`E_{all}`$, we have

$$E_{kin}+E_{int}+E_{ho}=E_{all},$$
(23)
where $`E_{ho}`$ is the potential energy of the condensates. For the central peak, $`E_{ho}`$ can be omitted safely. In the case of non-interacting model, $`E_{all}`$ will transform fully to the kinetic energy and potential energy of the condensates once the optical lattices are switched off. Due to the presence of the repulsive interaction between atoms, the maximum density of the central peak would be smaller than the result of the non-interacting model.

Note that the kinetic energy of the central peak can not be calculated through $`N\frac{\mathrm{}^2}{2m}\left(\sqrt{\left|\phi _0(x,r_{},t_m)\right|^2}\right)^2𝑑V`$, because the phase factor is different for different well (See Eq. (7)). From the uncertainty relation, assume $`E_{kin}1/R_x^2\alpha _x^4`$. In the presence of repulsive interaction, we have

$$E_{kin}=\frac{\alpha _x^4}{\alpha _{xideal}^4}E_{all}.$$
(24)
From Eqs. (22), (23) and (24), one obtains the following equation to determine $`\alpha _x^2`$:

$$\beta ^2+\theta \beta 1=0,$$
(25)
where $`\beta =(\alpha _x/\alpha _{xideal})^2`$. The value of $`\beta `$ reflects how the repulsive interaction between atoms reduces the density of the central peak. In the above expression, the dimensionless parameter $`\theta =E_{int}\left(\alpha _x=\alpha _{xideal}\right)/E_{all}`$. From Eq. (25), we have

$$\beta =\frac{\theta +\sqrt{\theta ^2+4}}{2}.$$
(26)
Now let us turn to discuss the total energy of the condensates which is necessary to calculate the value of $`\beta `$. The total energy of the condensates can be obtained through the sum of the energy of the condensates in each well before the optical lattices are switched off. Before the optical lattices are switched off, the normalized wave function in $`k`$th well takes the form

$$\phi _{0k}(x,r_{})=\phi _{0k}\left(x\right)\phi _{0k}\left(r_{}\right)$$
(27)
where

$$\phi _{0k}\left(x\right)=\left(\frac{1}{\sqrt{\pi }\sigma }\right)^{1/2}\mathrm{exp}\left[\frac{\left(xkd\right)^2}{2\sigma ^2}\right],$$
(28)
and $`\phi _{0k}\left(r_{}\right)`$ is the normalized wave function in the radial direction. From Eq. (27), we have

$`E_{all}={\displaystyle \underset{k=k_M}{\overset{k_M}{}}}\{N_k{\displaystyle }\phi _{0k}(x,r_{})[{\displaystyle \frac{\mathrm{}^2}{2m}}^2+{\displaystyle \frac{1}{2}}m\omega _{xe}^2x^2+{\displaystyle \frac{1}{2}}m\omega _{}^2r_{}^2]\phi _{0k}(x,r_{})dV`$

$$+\frac{gN_k^2}{2}\left|\phi _{0k}(x,r_{})|^4dV\right\},$$
(29)
where $`\omega _{xe}`$ is the effective harmonic frequency in the $`x`$direction of the well induced by the optical lattices. The last term in the above expression represents the interaction energy of the condensates in each well, and it is easy to verify that it can be omitted safely. In addition, the kinetic energy and potential energy in the radial direction can be also omitted because $`\omega _{xe}>>\omega _{}`$, and $`\mathrm{}\omega _{xe}>>\mu _0`$. In this case, one gets

$$E_{all}N\left(\frac{\mathrm{}^2}{4m\sigma ^2}+\frac{1}{4}m\omega _{xe}^2\sigma ^2\right).$$
(30)
For the experimental parameters in , from the formulas (22), (26), and (30), the calculation shows that $`\beta =0.80`$. This shows clearly that the repulsive interaction between atoms would reduce the density of the central peak at $`t_m`$.

With the variation of the parameter $`s`$,, the maximum density of the central peak can be also calculated based on the method given here. In the interacting model, the solid line in Fig. 6 shows the ratio between $`n\left(x=0,t_m,s\right)`$ and $`n\left(x=0,t_m,s=5\right)`$. The result of the non-interacting model is also shown in the figure.

## VI Discussion and Conclusion

In brief, the evolution process of the condensates is investigated after the optical lattices are switched off. We find that the density oscillation exhibits a phenomenon of decay and revival, based on the numerical result of the evolution of the density distribution. The decay of the density oscillation is interpreted as the emergence of the diffraction phenomenon, which is regarded as a consequence of interference from a lot of coherent expanding condensates. Due to the confinement of the harmonic potential, there is a periodic character of the density distribution, and it is this periodic character which leads to the revival of the density oscillation. In contrast to the condensates in the magnetic trap, there is no revival of the density oscillation, when both the magnetic and optical lattices are switched off. In addition, in the case of non-interacting model, it is shown that the maximum value of the density distribution at $`x=0`$ would be approximately $`30`$ times larger than the case when there are no optical lattices to induce the interference effect.

It is shown here that the repulsive interaction between atoms has an effect of reducing the maximum density of the central peak. In a real experiment in future, maybe the experimental result of the maximum density of the central peak would be smaller than the theoretical prediction given here, because it is possible that there is a loss of the total energy of the condensates during the process of the removal of the optical lattices. For the attractive interaction such as $`Li`$, the role of interaction would become very important when the atoms are confined by the combined potential. For example, when only the optical lattices are removed, based on the non-interacting model, the density at $`x=0`$ would increase largely due to the interference and confinement of the magnetic trap. In addition, due to the attractive interaction between atoms, the density of the central peak would increase rapidly. In this situation, it is possible that the condensates would collapse and even explode in a subsequent time, in analogous with the dynamic process of the collapsing and exploding atoms  by switching the interaction from repulsive to attractive.

## Acknowledgments

This work was supported by the Science Foundation of Zhijiang College, Zhejiang University of Technology, and Natural Science Foundation of Zhejiang Province. One of us (G. H.) is indebted to National Natural Science Foundation of China under Grant No. 19975019, and the French Ministry of Research for a visiting grant at Universitè Paris 7.

## Figure Caption

Fig. 1 Displayed is the density of the condensate at $`x=0`$ vs time $`t`$, after the optical lattice is switched off. Here the density $`n(x=0,t)`$ is in units of $`NA_n^2`$. We can see clearly in the figure that there is a decay and revival of the density oscillation. In addition, there is a periodicity of the density due to the confinement of the magnetic trap.

Fig. 2 (a)-(d) show the evolution of the density distribution with time $`t`$, after the optical lattices are switched off. The density distributions are shown at $`t=0`$, $`0.1\pi /\omega _x`$, $`0.3\pi /\omega _x`$, and $`0.5\pi /\omega _x`$. The emergence and motion of the $`n=\pm 1`$ peaks are clearly shown in these figures. Here the density distribution $`n(x,t)`$ is in units of $`NA_n^2`$, while the location $`x`$ is in units of $`d`$, i.e. the distance between two neighboring condensates.

Fig. 3 Displayed is the motion of the $`n=1`$ peak, after the optical lattices are switched off. Here the location $`x\left(t\right)_{n=1}`$ of $`n=1`$ peak is in units of $`d`$. The solid line is the result calculated from the classical harmonic motion given by Eq. (11). The squares show the motion of $`n=1`$ peak obtained from the numerical result given by Eq. (6). We see that the classical harmonic motion agrees quite well with the numerical result.

Fig. 4 Displayed is the density of the condensate at $`x=k_Md/2`$ vs time $`t`$, after the optical lattices are switched off. Here the density $`n(x,t)`$ is in units of $`NA_n^2`$. There is a phenomenon of the decay and revival of the density oscillation.

Fig. 5 Displayed is the time of the disappearance of the density oscillation for different locations. Here the location $`x`$ is in units of $`d`$, while the time $`t`$ is in units of the time scale $`t_w`$. The solid line is obtained from the formula (20), while the squares show the result obtained directly from the numerical result given by Eq. (6). We see that the analytical formula (20) can give a well description for the disappearance of the density oscillation.

Fig. 6 Displayed is the ratio between $`n\left(x=0,t_m,s\right)`$ and $`n\left(x=0,t_m,s=5\right)`$ for the interacting and non-interacting models.


text:
# PG1115+166 – a long period DA+DB binary.

## 1 Introduction

The chemical composition of the atmosphere of a white dwarf star reflects a balance between gravitational settling, which tends to produce a pure hydrogen atmosphere, and processes which “pollute” the atmosphere such as radiative levitation, convection and accretion (Fontaine & Wesemael 1987). This balance will change during the lifetime of a white dwarf as it cools and this is seen as a change in the relative numbers of hydrogen-rich (DA) and helium rich (non-DA) white dwarfs as a function of effective temperature. In white dwarfs with effective temperatures $`\mathrm{T}_{\mathrm{eff}}`$12 – 60 000K, gravitational settling dominates so about 80 percent of these white dwarfs have pure or nearly-pure hydrogen atmospheres. White dwarf stars with helium-rich atmospheres may appear as DO stars with He II lines ($`\mathrm{T}_{\mathrm{eff}}`$45 000K) or DB stars with He I lines (30 000K $`\mathrm{T}_{\mathrm{eff}}`$12 000K). The absence of non-DA white dwarfs with effective temperatures between 30 000K and 45 000K is known as the DB gap. The favoured explanation for the DB gap is that as the DO white dwarf cools to 45 000K, small amounts of hydrogen hidden in the atmosphere rise to the surface, eventually masking the underlying helium. A hydrogen layer mass of just $`10^{14}\mathrm{M}_{}`$ is enough to make the white dwarf look like a DA. Later, at $`\mathrm{28\hspace{0.17em}000}`$K, the convection zone in the underlying helium layer reaches the photosphere, mixing the hydrogen and helium so the white dwarf re-appears as a helium-rich, DB white dwarf. Those white dwarfs with thicker hydrogen layers remain as DA stars in this temperature range. In this scenario, $`20`$ percent of DA white dwarfs within the DB gap have thin hydrogen layers which are sufficient to mask the underlying helium.

The reason why $`20`$ percent of white dwarfs are born with very thin hydrogen envelopes is thought to be related to the born-again red giant phenomenon in which some stars form hot white dwarfs with thick hydrogen envelopes which undergo a late thermal pulse. This is caused by re-ignition of the helium shell which results in the stars rapidly expanding to red giant dimensions, e.g., Sakurai’s Object (Duerbeck 2000). The remaining hydrogen is quickly mixed into the helium layers and is burnt (Iben et al. 1983).

Those few white dwarfs which show a mixture of hydrogen and helium in their atmosphere are particularly interesting in the light of this debate. About 20 percent of DB stars show traces of hydrogen in their atmospheres and are classified as DBA (Shipman, Liebert & Green 1987). Much rarer are stars which show both strong hydrogen lines and strong helium lines. These are classified as DAB stars and only six of them have been identified. The prototype, GD 323, appears at the cool edge of the DB gap ($`\mathrm{T}_{\mathrm{eff}}`$= 28 750K) and appears to have a very thin hydrogen layer (10-17$`\mathrm{M}_{}`$) although no model can yet explain all the features of this peculiar object satisfactorily (Koester, Liebert & Saffer 1994). MCT 0128$``$3846 and MCT 0453$``$2933 appear to be binaries with a DA and a DB or DBA component (Wesemael et al., 1994). Holberg, Kidder & Wesemael (1990) detected a weak HeI 4471Å feature in G104$``$27 (WD 0612+177) but this was not confirmed by subsequent observations (Kidder et al. 1992). HS 0209+0832 is a DAB star in the middle of the DB gap which shows strong HeI lines whose strength varies from year-to-year (Heber et al. 1997) which may be explained by accretion onto the white dwarf from an inhomogeneous interstellar cloud (Wolff et al. 2000). Finally, PG 1115+166 was first reported as a DAB star by Burleigh et al. (2001). They note that the H$`\alpha `$ line shows a variable radial velocity but were unable to determine the orbital period. Bergeron & Liebert (2002) have analysed the optical spectrum of PG 1115+166 and found an excellent fit to the spectrum by assuming it is an unresolved DA+DB double degenerate binary. Their fit to the observed spectrum used a DA star with T$`{}_{\mathrm{eff}}{}^{}=22090`$K, $`\mathrm{log}g=8.12`$ combined with a helium-line DB star with T$`{}_{\mathrm{eff}}{}^{}=16210`$K, $`\mathrm{log}g=8.19`$. They estimate that both stars have masses of about 0.7$`\mathrm{M}_{}`$.

Those DAB stars which are binaries may be particularly interesting as they may be examples of stars which have undergone two mass transfer episodes, e.g., a common envelope phase (Iben & Livio 1993). The number of such close white dwarf pairs has increased rapidly in recent years. The majority of these binaries were identified as low mass white dwarfs ($`0.45\mathrm{M}_{}`$) from analysis of their Balmer lines. White dwarfs of such low mass are not formed in the standard picture of single star evolution within the lifetime of our Galaxy. It now appears that the majority of these binaries are the result of a common envelope phase in a star on the first giant branch (Marsh, Dhillon & Duck 1995). The properties of the binaries identified to-date, e.g., their orbital periods, masses and mass ratios, are strong tests of models for the formation of close binary white dwarfs and models for interacting binary stars in general (Nelemans et al. 2000, Nelemans et al. 2001). Identifying DAB stars may be a promising way to find more massive close binary white dwarf stars and so extend the range of masses over which observations can be used to test models of interacting binary stars.

In this paper we report radial velocity measurements of the Balmer lines and HeI lines of PG 1115+166. The Balmer lines show that this star is a binary with a period of 30.09d and the helium lines move in anti-phase to the hydrogen lines, so this is star is a DA+DB binary white dwarf. We also consider the formation of PG 1115+166.

## 2 Observations and Reductions

PG 1115+166 was first identified as a DAB star in February 1996 as part of a survey to find low mass white dwarfs (Moran 1999). It was re-observed several times thereafter as part of the programme at Southampton to find and characterize close binary white dwarfs. The Leicester group observed PG 1115+166 independently in April 1999 as a part of a programme to find white dwarfs with helium-rich atmospheres in the DB gap by observing hot white dwarfs undetected by EUV surveys and also identified it as a DAB star (Burleigh et al. 2001)

Observations of PG 1115+166 have been obtained with the intermediate dispersion spectrograph (IDS) on the 2.5m Isaac Newton Telescope (INT) and the dual-beam ISIS spectrograph on the 4.2m William Herschel Telescope (WHT) both on the Island of La Palma. A variety of gratings have been used to cover a variety of wavelengths at different times, as detailed in Table 1. The exposure times used were typically 10–30 minutes. Spectra of an arc lamp are taken before and after each target spectrum with the telescope tracking the star. None of the CCDs used showed any structure in unexposed images, so a constant bias level determined from a clipped-mean value in the over-scan region was subtracted from all the images. Sensitivity variations were removed using observations of a tungsten calibration lamp.

Extraction of the spectra from the images was performed automatically using optimal extraction to maximize the signal-to-noise of the resulting spectra (Marsh 1989). The arcs associated with each stellar spectrum were extracted using the profile determined for the stellar image to avoid possible systematic errors due to tilted arc lines. The wavelength scale was determined from a polynomial fit to measured arc line positions and the wavelength of the target spectra interpolated from the calibration established from the bracketing arc spectra. Uncertainties on every data point calculated from photon statistics are rigorously propagated through every stage of the data reduction. The normalized spectrum of PG 1115+166 is shown in Fig. 1 in two sections. The spectra used to produce this figure were the 1999 Dec WHT spectrum, the 1999 Apr WHT spectrum and the sum of the three WHT spectra from 2002 Jan. No attempt has been made to remove the telluric features near 6900Å and 7300Å.

## 3 Analysis

### 3.1 Radial velocities.

To measure the radial velocities we used least-squares fitting of a model line profile. This model line profile is the summation of one or more Gaussian profiles with different widths and depths but with a common central position which varies between spectra and is convolved with a Gaussian profile of an appropriate width to model the resolution of each spectrum.

One must be careful to allow for the asymmetry of the line due to pressure shifts, particularly when dealing with the wings of the higher order Balmer lines (Grabowski, Halenka & Madej 1987) or the HeI lines (Beauchamp, Wesemael & Bergeron 1997). To account for these asymmetries in the Balmer lines, we add a 3rd or 4th order polynomial in the model line profile. We also exclude data from the fitting process which is more than 1000 km s-1 from the rest wavelength of the lines measured. This reduces both the effect of the pressure shift and the possibility of blending with other lines. The radial velocity measurement is determined primarily by the centre of the line and the fitting process is identical for all lines, so any remaining pressure shifts would only be seen as an offset of a few km s-1 between the radial velocities measured for different Balmer lines. This is small compared to the typical uncertainties in the measured radial velocities.

The pressure shifts in the much weaker HeI lines are harder to deal with, so we simply report the radial velocities measured using a single Gaussian profile and a low-order polynomial for these lines and include a warning here that there are systematic shifts between the velocities measured for different HeI lines.

We first normalize the spectra using a polynomial fit to the continuum either side of the line of interest. We use a least-squares fit to one of the spectra to determine the model line profile. A least squares fit of this profile to each spectrum in which the position of the line and the polynomial coefficients are the only free parameters gives the measured radial velocities in Table 2. We tried three different initial values for the radial velocity in the least-squares minimization to ensure we found the optimum value of the radial velocity and excluded the results for any spectra where the optimum value could not be clearly identified.

The uncertainties on each radial velocity measurement are calculated by propagating the uncertainties on every data point in the spectra right through the data reduction and analysis. There are inevitably systematic errors in our radial velocity measurements, particularly since we are using spectra obtained with a variety of instruments. Maxted et al. (2000) found that these systematic differences are no more than $`1`$ km s-1 for the majority of the spectrographs used in this study. To allow for this, we include an additional uncertainty per radial velocity measurement equivalent to 1/40 of the resolution, typically about 1km s-1. This has been added in quadrature to the uncertainties given in Table 2. The radial velocity measurements of the Balmer lines are clearly variable. Of the HeI line measurements, only the measurements of the HeI 6678 lines are sufficiently numerous and accurate to show any sign of variability. There are too few accurate radial velocity measurements of the other HeI lines for this variability to be apparent. To quantify this we calculated their weighted mean radial velocity, which is the best estimate of the radial velocity assuming this quantity is constant. We then calculated the $`\chi ^2`$ statistic for this “model”, i.e., the goodness-of-fit of a constant to the observed radial velocities. We can then compare the observed value of $`\chi ^2`$ with the distribution of $`\chi ^2`$ for the appropriate number of degrees of freedom and find the probability of obtaining the observed value of $`\chi ^2`$ or higher from random fluctuations of constant value, $`p=0.023`$. This suggests that the HeI 6678 lines are likely to be variable, but the evidence from this test is not conclusive.

### 3.2 The orbital period

We used the 57 radial velocities measured from the Balmer lines to establish the orbital period. We used least-squares to fit sine waves at 500 000 trial frequencies uniformly distributed over the range 0 to 10 cycles per day and recorded the value of chi-squared for each trial frequency. Low values of chi-squared occur for orbital periods $`P31`$d with the lowest value of chi-squared being 68.3 for $`P=30.0873`$d. The next most likely orbital period is near the one cycle-per-day alias of this period, $`P=1.0316`$d, which has a chi-squared value of 87.3.

In order to test the robustness of our period determination, we randomly selected 50 of the radial velocity measurements and fitted sine waves at 100 000 trial frequencies uniformly distributed over the range 0 to 1.2 cycles per day and recorded the value of the period at which the minimum value of chi-squared occurs. We repeated this process 1000 times and found that in 973 trials this reduced dataset produced a minimum value of chi-squared at an orbital period near 30.09d. Most of the other trials favoured orbital periods near 1 day which, as we have seen, gives a poor fit to the complete data set.

As an alternative test, we used a Monte Carlo simulation to estimate how often a binary with an orbital period of 1.0316d will appear to have an orbital period near 30.09d given radial velocity meausurements with the same temporal sampling and accuracy as those presented here. We used a least-sqaures fit to our data to find the optimum period near 1.0316d and the corresponding semi-amplitude of a sine wave fit by least squares. We generated 1000 synthetic data sets with the same period and semi-amplitude as this fit and with the same temporal sampling as our data. We added noise to the synthetic data using a psuedo-random normal deviates multiplied by the uncertainties associated with each datum. For each of these synthetic data sets we fitted sine waves at 100 000 trial frequencies uniformly distributed over the range 0.01 to 1.2 cycles per day and recorded the value of the period at which the minimum value of chi-squared occurs. We also recorded the best values of chi-squared which occur within 10 percent of $`P=1.0316`$d and $`P=30.09`$d. Of the 1000 trials, 983 trials resulted in a minimum value of chi-squared near 1.0316d and only 17 gave a minimum value of chi-squared at a period near 30.09d. Furthermore, in all the trials for which $`P30.09`$d gave a lower value of chi-squared, the difference in the minimum chi-squared value for $`P1.0316`$d and $`P30.09`$d was usually only 1–3 and was always less than 8. In summary, this test suggests that there is a less than 1/1000 chance of a binary with $`P=1.0316`$d being incorrectly identified as a binary with $`P=30.09`$d from our data provided the least-squares fits of sine waves at the two periods result in a difference in chi-squared of 8 or more. The difference in chi-squared for our actual data is 19, so this condition is easily satisfied in this case.

These simple test both suggests that period determination is robust and that we can be confident that the orbital period is $`P30.09`$d.

### 3.3 Spectroscopic orbit

If PG 1115+166 is genuinely a binary DA+DB white dwarf we expect that the HeI lines should move in anti-phase to the Balmer lines. We investigated this scenario using a least-squares fit of a sine wave of the form $`v_i=\gamma _\mathrm{H}+K_\mathrm{H}\mathrm{sin}(2\pi (T_iT_0)/P`$ to the measured radial velocities, $`v_i`$, for all the Balmer lines and a sine wave of the form $`v_i=\gamma _{\mathrm{HeI}}+K_{\mathrm{HeI}}\mathrm{sin}(2\pi (T_iT_0)/P`$ to the 23 measured radial velocities for the HeI 6678 line. We used both data sets in a simultaneous fit to all these data, an approach which has the merit of using all the available data to determine the parameters $`T_0`$ and $`P`$. The results are given in Table 4 and are shown in Fig. 2. We have excluded data for the Balmer lines with uncertainties larger than 20 km/s and data for the HeI 6678 line with uncertainties larger than 40 km/s from Fig. 2 for clarity. The value of $`K_{\mathrm{HeI}}`$ derived is about 4 times greater than its uncertainty, so we can be fairly sure that the HeI 6678 line moves in anti-phase to the Balmer lines as expected for a DA+DB binary. We can certainly rule out the possibility that PG 1115+166 is a single DAB star since that would require the HeI lines and the Balmer lines to move in phase, which is clearly not the case. Also given in Table 4 are the results of a similar least-squares fit including the data for all the HeI lines. This reduces the uncertainty in the value of $`K_{\mathrm{HeI}}`$, but we prefer the results of fitting only the HeI 6678 line until the systematic errors introduced by pressure shifts in HeI lines is better understood. In any case, the results of the two least-squares fits are nearly identical.

The mass ratio implied by the parameters in Table 4 is $`q=M_{\mathrm{DA}}/M_{\mathrm{DB}}=0.84\pm 0.21`$. This is consistent with the conclusion of Bergeron & Liebert (2002) that the masses of the two stars are similar. We also find $`M_{\mathrm{DA}}\mathrm{sin}^3i=0.43\pm 0.15\mathrm{M}_{}`$, $`M_{\mathrm{DB}}\mathrm{sin}^3i=0.52\pm 0.12\mathrm{M}_{}`$ and $`a\mathrm{sin}i=40\pm 5\mathrm{R}_{}`$, where $`i`$ is inclination of the binary and $`a`$ is the separation of the stars. The masses measured by Bergeron & Liebert suggest that the inclination of the binary is approximately 60 and that the stars are separated by about 45$`\mathrm{R}_{}`$.

## 4 The formation of PG 1115+166

The basic mechanism for the formation of most white dwarfs is fairly well understood. They are, in general, remnants of intermediate mass stars (1–8$`\mathrm{M}_{}`$) which have evolved through the asymptotic giant branch (AGB) phase. Extreme mass loss towards the end of this phase exposes the core of the AGB star, which cools to form a white dwarf. Stars with larger initial masses are expected to form more massive white dwarfs and there is some observational evidence for this initial-final mass relation (IFMR; Weidemann 2000). White dwarfs with masses of 0.7$`\mathrm{M}_{}`$ are expected to form from stars with initial masses of about 3$`\mathrm{M}_{}`$.

Models of stars on the AGB are uncertain, but the radius of a star towards the end of this phase is several hundred solar radii (Girardi et al. 2000). This is much larger than the current separation of the stars so it appears that the orbit of PG 1115+166 is now substantially smaller than it was prior to the formation of the white dwarfs. Orbital shrinkage in long period binaries is explained by the “common envelope” scenario. In this scenario, the expanding red giant star comes into contact with its Roche lobe and begins to transfer mass to its companion star. This mass transfer is highly unstable, so a common envelope forms around the companion and the core of the red giant. The drag on the companion orbiting inside the common envelope leads to extensive mass loss and dramatic shrinkage of the orbit (Iben & Livio 1993). We expect the orbit of a post-CE binary to be circular, and have assumed this to be so in the case of PG 1115+166 throughout this paper. Much higher quality data than that presented here would be required to verify this assumption.

Common envelope evolution is poorly understood so it is useful to consider the limits imposed on the common envelope (CE) phase by the observed properties of PG 1115+166. Firstly, we note that published models for the population of white dwarf binaries do not predict any systems with the combination of high mass and long period seen in PG 1115+166. This is shown in Fig. 3 for the models of Iben, Tutukov and Yungelson (1997), but the same applies to the models of Nelemans (2001) and Han (1998). It is primarily the combination of high mass and long orbital period which is not explained by these models. A full exploration of all the possible formation scenarios for PG 1115+166 is beyond the scope of this paper, not least because there may have been up to three mass transfer episodes – two associated with the red giant phase prior to the formation of each white dwarf and a third associated with the born-again red giant phase of the DB white dwarf. These mass transfer episodes may not have been CE phases but may instead have involved stable mass transfer, i.e., an Algol-like phase. In this section we simply outline the observational constraints which can be applied to these scenarios and apply these constraints to the simplest case of a single CE phase.

Apart from the present day masses of the white dwarfs ($`0.7\mathrm{M}_{}`$) and their current separation ($`a_f45\mathrm{R}_{}`$), we also know the cooling ages of the stars are $`6\times 10^7`$y for the DA white dwarf and $`2.2\times 10^8`$y for the DB white dwarf (Bergeron & Liebert 2002). The cooling ages are the time elapsed since the AGB phase. If the progenitors of the white dwarfs were formed at the same time, this imposes a limit on the initial masses of the progenitor stars as follows. If we assume an initial mass for the progenitor of the DB star, $`M_1^{}`$, we can use an appropriate stellar model to find the lifetime of this star, $`\tau _1`$, i.e., the time taken since it formed for the star to reach the end of AGB phase. We then know that the progenitor of the DA star will take a further $`2.2\times 10^8\mathrm{y}6\times 10^7\mathrm{y}=1.6\times 10^8`$y to reach the AGB phase. We can then use grid of stellar models to find the initial mass, $`M_2^{}`$, of a star whose lifetime is $`\tau _2=\tau _1+1.6\times 10^8`$y, i.e, the mass of the progenitor of the DA star. The relationship between these initial masses is shown in Fig. 4 based on the models of Girardi et al. (2000) for Z=0.019 (solid lines) and Z=0.004 (dashed lines). Although the stars have not evolved as single stars, this argument regarding the cooling age difference between the stars remains valid because the AGB phase is much shorter than the lifetime of the star which precedes it, i.e., the stars spend most of their lives evolving as single stars. However, we should be aware that the IFMR may not apply to the white dwarfs in PG 1115+166 in the sense that a star may produce a white dwarf which is less massive than expected if it is involved in a CE phase which interrupts the growth of the future white dwarf forming in its core. Also, this constraint may not be valid if there has been significant accretion onto one of the stars. Similarly, if the born-again red giant phase of the DB star was short, the DB white dwarf may not have been re-heated to the temperatures expected for a white dwarf emerging from its first AGB phase. In this case the cooling age does not measure either the time since the AGB phase, nor the time since the born-again read giant phase, but is somewhere between these two timescales.

The observation that one of the stars is a DB white dwarf is also a constraint on possible formation scenarios for PG 1115+166. Any feasible scenario must allow for the DB white dwarf to undergo a born-again red giant phase or include an alternative explanation for how this star came to be hydrogen deficient. A stronger constraint, perhaps, is how this star avoids accreting hydrogen onto it surface despite appearing to be older than its companion which has, presumably, lost large quantities of hydrogen rich material.

With this information we consider the simplest case of a single CE phase in which the AGB star progenitor of the DA star with mass $`M_2`$ and radius $`R_2`$ interacts with the DB star when it is close to or exceeds its Roche lobe. Note that $`M_2`$ will be less than $`M_2^{}`$ if any mass loss occurs prior to the interaction. However, this mass loss cannot have been very extensive in this scenario because the DB white dwarf cannot have gained more than $`10^{14}\mathrm{M}_{}`$ so we assume $`M_2=M_2^{}`$. If the initial semi-major axis of the stars’ orbit is $`a_i`$ and the current separation is $`a_f`$, the change in the orbital binding energy is

$$\mathrm{\Delta }E_{\mathrm{orb}}=\frac{GM_2M_{\mathrm{DB}}}{2a_i}\frac{GM_{\mathrm{DA}}M_{\mathrm{DB}}}{2a_f}$$

Some fraction of this energy $`\alpha _{\mathrm{CE}}`$ is used to expel the envelope of the AGB star whose binding energy is given by

$$E_{\mathrm{env}}=\frac{GM_2(M_2M_{\mathrm{DA}})}{\lambda R_1}$$

So from the definition of $`\alpha _{\mathrm{CE}}`$ we have

$$\alpha _{\mathrm{CE}}\lambda \left(\frac{GM_{\mathrm{DA}}M_{\mathrm{DB}}}{2a_f}\frac{GM_2M_{\mathrm{DB}}}{2a_i}\right)=\frac{GM_2(M_2M_{\mathrm{DA}})}{R_1}$$

If the orbit prior to the interaction has an eccentricity $`e`$, then the $`R_1=r_La_i(1e)`$ at the time of interaction, where $`r_L`$ is the radius of the Roche lobe relative to the separation of the stars. This assumes that the interaction occurs at periastron, which is reasonable given that the orbital period ($`10`$y) is much shorter than the evolutionary timescale of an AGB star.

In Fig. 5 we show lines of constant $`\alpha _{\mathrm{CE}}\lambda `$ as a function of $`M_2`$ and $`R_2/(1e)`$ for $`M_{\mathrm{DA}}=M_{\mathrm{DB}}=0.7\mathrm{M}_{}`$ and $`a_f=45\mathrm{R}_{}`$. The maximum radius of an AGB star with an initial mass of $`3\mathrm{M}_{}`$ is a few hundred solar radii (Bloecker 1995), so unless the eccentricity of the orbit prior to the CE phase was quite extreme ($`e0.5`$) and/or the initial mass of the AGB star was lower than suggested by the IFMR ($`2\mathrm{M}_{}`$), we see that $`\alpha _{\mathrm{CE}}\lambda 1.0`$. This requires a large value of $`\lambda `$ because $`\alpha _{\mathrm{CE}}<1`$ in the absence of other sources of energy to eject the envelope. This strong limit on $`\lambda `$ is a consequence of the combination of long orbital period and high mass in PG 1115+166. Large values of $`\lambda `$ for AGB stars are predicted by Dewi & Tauris (2000) provided the calculation of the binding energy of the envelope includes the internal energy stored in the envelope, e.g., as ionized hydrogen. This implicitly assumes that this potential energy is efficiently converted into the kinetic energy of the envelope during the common envelope phase. Several issues arising from this scenario remain to be explored, e.g., how does the DB star avoid gaining more than $`10^{14}\mathrm{M}_{}`$ of hydrogen during the CE phase, why was there no CE phase when the progenitor of the DB star went through its AGB phase and where does the born-again red giant phase fit into this scenario? It is beyond the scope of this paper to tackle these questions, but they do demonstrate that our observations of PG 1115+166 provide a challenging test of models for binary star evolution and the formation of DB white dwarfs.

## 5 Conclusion

We have used radial velocity measurements of the Balmer lines and HeI lines of PG 1115+166 to show that this star is a DA+DB binary white dwarf with an orbital period of 30.09d. Published models for the formation of white dwarf – white dwarf binaries do not predict any white dwarfs with the characteristics of PG 1115+166. We consider the observational constraints on scenarios for the formation of PG 1115+166. We apply these to the simplest case of a single common-envelope phase. We find that this requires that some part of the internal energy in the envelope of the AGB star contributes to the ejection of the common envelope, but several questions regarding this scenario remain to be answered. We conclude that PG 1115+166 is a key object for testing models of binary star evolution. It may also be a key object for our understanding of the formation of DB white dwarfs.

## Acknowledgments

PFLM would like to thank Jacco van Loon for his helpful discussions regarding AGB stars and Phillip Podsiadlowski for sharing his thoughts on the formation of PG 1115+166. The Isaac Newton Telescope is operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias.