text:
# Quantum dots with two electrons: Singlet-triplet transitions

## I Introduction

The spin state of the many electron ground-state is determined by the interplay between the kinetic and the electrostatic interactions. According to Hund’s law the ground-state of electrons in a partially filled shell of an atom has the maximal possible spin, in order to minimize the electrostatic repulsion. On the other hand, Anderson’s super-exchange antiferromagnetic interaction, which favors zero total spin, arises from the reduction in the ground-state energy brought about by the hopping of electrons between adjacent ions. Another example is realized in gases made of homonuclear diatomic molecules, in which the distance between the nuclei determines the magnetic properties of the gas; in other words, the ground-state of the molecule can be either a triplet or a singlet.

The possibility to study these rules as function of controlled parameters, and in particular to observe deviations from them, has been opened up in recent years with the intensive experimental and theoretical investigations of small confined systems. Devices based on quantum dots formed in GaAs heterostructures allow to probe the electronic states in situations where there are only a few electrons in the system, as function of, e.g., the number of electrons in the sample, the spacing of the single-electron energy levels, the internal structure of the dots (e. g. the distance and the coupling between internal parts) or the coupling of the system to external leads. Thus, modulating the single-particle spectrum by a magnetic field, applied perpendicularly to the plane of the electron gas,  has produced a structure of the conductance peaks that has been interpreted in terms of singlet-triplet transitions of the last two electrons in the dot  (see also Ref. ). This has been observed on vertical quantum dots,  and also in lateral ones.  The ground-state spin of chaotic quantum dots has been studied by tracing the conductance peak spacing as function of a weak parallel magnetic field (which couples primarily to the spins).  A numerical investigation of such a configuration analyzed the influence of the exchange interaction on the peak structure.  In a similar fashion, the singlet-triplet transitions in such dots have been attriubted to avoided-crossings in the many-electron states, and the relation between those and kinks in the conductance pattern has been explored.

When the device consists of two or more coupled dots, another controllable parameter comes into the play: the interdot coupling. Experimentally,  it has been found that this coupling shows up in the conductance peak positions. Theoretically, the effect of the interdot distance on the states of a few electrons confined in a parabolic potential has been analyzed using mean-field methods,  the Kohn-Sham equation,  and numerical diagonalization. Different spin states of double-quantum-dots systems have been also studied using the numerical renormalization group method.

Here we present an exact analytical solution for the eigenvalues and eigenfunctions of two electrons on a double-dot system, which is coupled to two single-channel leads. The two electrons interact while they are on the dots, and we include in the calculation the direct Coulomb and the exchange interactions between the dots, and the on-site Hubbard interactions on the dots. We obtain the condition for the singlet-triplet transitions of the ground-state energies, and in particular examine the role of the delocalization effect of the interactions: It has been found in our previous works,  as well as in other studies,  that the interplay between the hybridization of the localized single-particle states on the dot with the propagating states on the leads, and the electron-electron interactions, may promote one of the electrons to the continuum. If the system had only one doubly bound state, then the above ‘promotion’ results in a ‘delocalization’ of the ground state. Once delocalized, the singlet and the triplet states become degenerate.

Our earlier work included only the case with on-site Coulomb repulsion. In that case, the bound ground state is a singlet. As the coupling of the dots with the leads was increased, the parameter range where this bound state exists was decreased, and the delocalization was accompanied by a transition from the singlet to a degenerate singlet-triplet ‘metallic’ state. In the present paper we add the effects of direct and exchange interactions. In the absence of the coupling to the leads, it is well known  that the ground state changes from a singlet to a triplet upon increasing the exchange interaction. Here we investigate what happens to that singlet-triplet transition in the presence of the coupling to the leads.

After discussing our model Hamiltonian in Sec. II, we present the exact two-electron solution in Sec. III. Section IV then discusses the singlet-triplet transition, mainly for a special choice of the parameters where it is easiest to explore the solution analytically. Our results are then summarized in Sec. V.

## II Hamiltonian

The Hamiltonian of the model reads

$`=_{\mathrm{sp}}+_\mathrm{c}.`$ (1)
The single-particle tight binding Hamiltonian is

$`_{\mathrm{sp}}={\displaystyle \underset{\mathrm{i}\sigma }{}}ϵ_\mathrm{i}c_{\mathrm{i}\sigma }^{}c_{\mathrm{i}\sigma }{\displaystyle \underset{\mathrm{ij}\sigma }{}}t_{\mathrm{ij}}c_{\mathrm{i}\sigma }^{}c_{\mathrm{j}\sigma },`$ (2)
where $`c_{\mathrm{i}\sigma }^{}`$ creates an electron with spin $`\sigma `$ on site i. The site energies $`ϵ_\mathrm{i}`$ are different from zero only on the ‘quantum dots’, which can be viewed as ‘impurities’ on the lattice. In what follows, we specifically consider two neighboring quantum dots, which will be denoted by $`\mathrm{}`$ and $`r`$, and study the symmetric case where $`ϵ_{\mathrm{},r}=ϵ_0`$. The hopping matrix elements $`t_{\mathrm{ij}}`$ are divided into three kinds: the hopping among the ‘dots’ is denoted by $`t_D`$ ($`t_\mathrm{}r`$), the hopping between a dot and a neighboring site on the ‘lead’ is denoted by $`t_0`$ (e. g. $`t_\mathrm{}\mathrm{i}`$ for i $`r`$), and the hopping between neighboring ‘lead’ sites is denoted by $`t=1`$, setting the units of energy.

The Coulomb interactions are assumed to exist only among electrons which sit on the dots. Generally, this interaction has the form

$`_\mathrm{c}={\displaystyle \underset{\mathrm{ijmn}}{}}\mathrm{\Gamma }_{\mathrm{ij}}^{\mathrm{mn}}{\displaystyle \underset{\sigma \sigma ^{}}{}}c_{\mathrm{i}\sigma }^{}c_{\mathrm{j}\sigma ^{}}^{}c_{\mathrm{m}\sigma ^{}}c_{\mathrm{n}\sigma },`$ (3)
where

$`\mathrm{\Gamma }_{\mathrm{ij}}^{\mathrm{mn}}={\displaystyle 𝑑𝐫𝑑𝐫^{}v(𝐫𝐫^{})\phi _\mathrm{i}^{}(𝐫)\phi _\mathrm{j}^{}(𝐫^{})\phi _\mathrm{m}(𝐫^{})\phi _\mathrm{n}(𝐫)},`$ (4)
while $`\phi _\mathrm{i}(𝐫)`$ is the Wannier wave function of an electron localized at site i. It is usually assumed that the dominant terms will be those in which the indices are equal pair-wise.  Here we follow the parametrization used in magnetic coupling studies,  in which the case $`\mathrm{i}=\mathrm{j}=\mathrm{m}=\mathrm{n}`$ is treated separately. Neglecting the coefficient with $`\mathrm{i}=\mathrm{j}\mathrm{m}=\mathrm{n}`$, which (when negative) leads to the superconductivity vertex, we are left with three possible parameters:

$`2\mathrm{\Gamma }_{\mathrm{ii}}^{\mathrm{ii}}`$ $`=`$ $`\mathrm{U},`$ (5)
$`2\mathrm{\Gamma }_{\mathrm{ij}}^{\mathrm{ji}}`$ $`=`$ $`\mathrm{V},2\mathrm{\Gamma }_{\mathrm{ij}}^{\mathrm{ij}}=\mathrm{K},\mathrm{i}\mathrm{j},`$ (6)
where $`\mathrm{U},\mathrm{V}`$ and $`\mathrm{K}`$ are the (intra-dot) Hubbard, (inter-dot) direct, and (inter-dot) exchange interactions, respectively. Then

$`_{\mathrm{int}}`$ $`=`$ $`\mathrm{U}{\displaystyle \underset{\mathrm{i}}{}}\widehat{n}_\mathrm{i}\widehat{n}_\mathrm{i}`$ (7)
$`+`$ $`{\displaystyle \frac{1}{4}}(2\mathrm{V}\mathrm{K}){\displaystyle \underset{\mathrm{i}\mathrm{j}}{}}\widehat{n}_\mathrm{i}\widehat{n}_\mathrm{j}\mathrm{K}{\displaystyle \underset{\mathrm{i}\mathrm{j}}{}}𝐒_\mathrm{i}𝐒_\mathrm{j},`$ (8)
where

$`\widehat{n}_{\mathrm{i}\sigma }`$ $`=`$ $`c_{\mathrm{i}\sigma }^{}c_{\mathrm{i}\sigma },\widehat{n}_\mathrm{i}={\displaystyle \underset{\sigma }{}}\widehat{n}_{\mathrm{i}\sigma },`$ (9)
$`𝐒_\mathrm{i}`$ $`=`$ $`{\displaystyle \underset{\sigma \sigma ^{}}{}}c_{\mathrm{i}\sigma }^{}\stackrel{}{\sigma }_{\sigma \sigma ^{}}c_{\mathrm{i}\sigma ^{}}`$ (10)
and $`\stackrel{}{\sigma }_{\sigma \sigma ^{}}`$ is the vector of Pauli’s spin matrices. It is important to note that with the above approximation for the Coulomb interaction vertex, the interaction Hamiltonian becomes spin-dependent.

A similar Hamiltonian has been used to study the effects of the spin states on the conducting properties of confined mesoscopic quantum dots.  In these studies, $`\mathrm{i}`$ represents some orbital state on the dot. In Ref.  it has been argued that when the Thouless conductance of the confined (mesoscopic) system is large, the Coulomb vertices can be parametrized in terms of two coupling constants, independent of the orbital indices. In terms of our parameters, these are the charging energy, $`\mathrm{E}_c=(2\mathrm{V}\mathrm{K})/4`$, and the exchange energy, $`\mathrm{J}=\mathrm{K}`$, while $`\mathrm{U}=\mathrm{V}+\mathrm{K}`$. It is straightforward to apply our results to that case.

In our special case of the double dot, the energies $`\mathrm{\Gamma }`$ in Eq. (6) differ from zero only for $`\mathrm{i},\mathrm{j}=\mathrm{},r`$, and thus

$`_{\mathrm{int}}`$ $`=`$ $`\mathrm{U}(\widehat{n}_{\mathrm{}}\widehat{n}_{\mathrm{}}+\widehat{n}_r\widehat{n}_r)`$ (11)
$`+`$ $`(\mathrm{V}{\displaystyle \frac{\mathrm{K}}{2}})\widehat{n}_{\mathrm{}}\widehat{n}_r2\mathrm{K}𝐒_{\mathrm{}}𝐒_r.`$ (12)
## III The two-electron ’molecule’

We will confine ourselves to a double-dot system (a ‘molecule’), containing two electrons. The system is modeled by two identical single-level impurities, each having an on-site single-particle energy level $`ϵ_{\mathrm{}}=ϵ_r=ϵ_0`$. The two impurities are coupled to one another by the interdot matrix element $`t_D`$. When the molecule is isolated ($`t_0=0`$), there are three degenerate triplet states, of energy $`2ϵ_0+\mathrm{V}\mathrm{K}`$, and three singlet states. One of the latter has the energy $`2ϵ_0+\mathrm{U}`$, while the other two energies are

$`\mathrm{E}=2ϵ_0+{\displaystyle \frac{\mathrm{U}+\mathrm{V}+\mathrm{K}}{2}}\pm \sqrt{4t_D^2+\left({\displaystyle \frac{\mathrm{U}\mathrm{V}\mathrm{K}}{2}}\right)^2}.`$ (13)
Examination of the interaction matrix elements, Eqs. (4) and (6), shows that $`\mathrm{U}>\mathrm{V}>0`$ and $`\mathrm{U}>\mathrm{K}`$. Moreover, the exchange interaction $`\mathrm{K}`$ involves the square of the overlap matrix element of the two impurities, and hence is of the order of $`t_D^2`$. Therefore one concludes that the lowest singlet energy is given by the minus sign in (13). The triplet states energy is below the lowest singlet one provided that

$`2t_D^2<\mathrm{K}(\mathrm{U}\mathrm{V}+\mathrm{K}).`$ (14)
One notes that when the direct and exchange Coulomb vertices V and K are disregarded, the energy of the singlet state is lowered by the interdot kinetic energy, producing the Anderson super-exchange antiferromagnetic coupling,  $`4t_D^2/\mathrm{U}`$. The triplet state becomes the ground-state once the exchange energy $`\mathrm{K}`$ wins over this kinetic energy.  Finally, with the choice of parameters in which $`\mathrm{U}=\mathrm{V}+\mathrm{K}`$ (see above), the condition becomes $`t_D^2<\mathrm{K}^2`$ (see also Ref. ). Both $`t_D`$ and K decay exponentially with the interdot distance. However, Eq. (4) indicates that K is roughly of order $`t_D^2`$, and thus decays faster. Therefore, one might expect a transition from the triplet to the singlet ground state as this distance increases.

The main purpose of the present paper is to study how the criterion for the singlet-triplet transitions, Eq. (14), is modified when the double quantum dot is connected to external leads, on which it is assumed that there are no electron-electron interactions. In other words, we study the changes in these transitions due to coupling the dot to a continuum of propagating states. A similar problem was discussed in Ref. , in the framework of the mean-field approximation, within a single-lead geometry. Here we first derive exactly the ground-state energy of two electrons, and then compare the singlet and the triplet ground-state energies.

For further simplification, we describe the external leads by single channel one dimensional chains, with nearest neighbor hopping $`t=1`$. In our previous work,  we have shown that the spectrum and the wave functions of two interacting electrons can be obtained in terms of the energy spectrum and the wave functions of the single-particle Hamiltonian. For our model, cf. Eq. (2), the latter can be divided into even and odd solutions, denoted by “e” and “o”. Of particular interest here will be the regions in the parameter plane $`\{ϵ_0\gamma \}`$, where $`\gamma =t_0^2`$, where the spectrum has bound states, see Fig. 1. For $`ϵ_0<\gamma 2+t_D`$, there appears one (even) bound (single-particle) state below the band of propagating states, of energy $`ϵ_\mathrm{B}^\mathrm{e}<2`$; For even smaller values of $`ϵ_0`$, such that $`ϵ_0<\gamma 2t_D`$ there appears a second (odd) bound state below the band, of energy $`ϵ_\mathrm{B}^\mathrm{o}`$, with

$`ϵ_\mathrm{B}^{\mathrm{e},\mathrm{o}}`$ $`=`$ $`{\displaystyle \frac{2\gamma }{1\gamma }}\left({\displaystyle \frac{ϵ_0^{\mathrm{e},\mathrm{o}}}{2}}\right)+{\displaystyle \frac{\gamma }{1\gamma }}\sqrt{\left({\displaystyle \frac{ϵ_0^{\mathrm{e},\mathrm{o}}}{2}}\right)^2+\gamma 1},`$ (15)
$`ϵ_0^\mathrm{e}`$ $`=`$ $`ϵ_0t_D,ϵ_0^\mathrm{o}=ϵ_0+t_D.`$ (16)
Similarly, for $`ϵ_0>2\gamma t_D`$ there appears the first bound state above the band, while for $`ϵ_0>2\gamma +t_D`$ there are two bound states above the band. Clearly, a necessary condition to have the two electrons bound in a triplet state is the existence of two distinct single-electron bound states, so that each electron occupies a different ‘state’. For simplicity, we shall confine ourselves to the case in which both occur below the continuum, i.e., to the lowest region II in Fig. 1.

It has been shown in Refs.  and  that when interactions between the two electrons are allowed, they may lead to the delocalization of one (or both) of the electrons from the doubly bound states. Thus, for example,  when V=K=0, in part of region I in Fig. 1 one finds that the doubly bound ground state is replaced by a ground state in which one of the electrons is shifted to the band. In these parameter regions there are no two-electron bound states, and the system may be called ‘metallic’. Indeed, such interaction-induced delocalization effects have been recently observed in quantum dots containing two potential minima.  This effect has also been found for two interacting electrons moving in a one-dimensional periodic structure,  and for spinless fermions on strongly disordered chains.  In contrast to region I, in region II there potentially exist two doubly bound states, i. e. the singlet and the triplet, and our goal is to study which of these states is the ground state in the presence of this interaction-induced delocalization.

To solve for the spectrum of the two interacting electrons we proceed as follows. Let us denote the eigenstates and the eigenvalues of the single-particle Hamiltonian (2) by $`\varphi _a(n)`$ and $`ϵ_a`$, respectively, where $`n`$ is the site index. Consider first the two-electron states with spin $`S_z=\pm 1`$,

$`|\mathrm{\Psi }^\sigma `$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{2}}}{\displaystyle \underset{ab}{}}\mathrm{X}_{ab}^\sigma c_{a\sigma }^{}c_{b\sigma }^{}|0`$ (17)
$`=`$ $`\sqrt{{\displaystyle \frac{1}{2}}}{\displaystyle \underset{ab}{}}\mathrm{X}_{ab}^\sigma {\displaystyle \underset{nn^{}}{}}\varphi _a^{}(n)\varphi _b^{}(n^{})c_{n\sigma }^{}c_{n^{}\sigma }^{}|0`$ (19)
$`\sqrt{{\displaystyle \frac{1}{2}}}{\displaystyle \underset{nn^{}}{}}\mathrm{X}_{nn^{}}^\sigma c_{n\sigma }^{}c_{n^{}\sigma }^{}|0,\sigma =\pm 1,`$
where the amplitudes $`\mathrm{X}_{nn^{}}^\sigma =\mathrm{X}_{n^{}n}^\sigma `$ are antisymmetric in the coordinates (the origin is half way between the dots), and $`_{nn^{}}|\mathrm{X}_{nn^{}}^\sigma |^2=1`$. The Schrödinger equation for the two electrons then yields

$`\mathrm{X}_{nn^{}}^\sigma `$ $`=`$ $`\mathrm{X}_r\mathrm{}^\sigma (\mathrm{V}\mathrm{K})`$ (21)
$`\times {\displaystyle \underset{ab}{}}{\displaystyle \frac{\varphi _a(r)\varphi _b(\mathrm{})\varphi _a(\mathrm{})\varphi _b(r)}{\mathrm{E}ϵ_aϵ_b}}\varphi _a^{}(n)\varphi _b^{}(n^{}).`$
Hence, the two-particle energies, E, are given by the solutions of the equation

$`{\displaystyle \frac{1}{\mathrm{V}\mathrm{K}}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{ab}{}}{\displaystyle \frac{|\varphi _a(\mathrm{})\varphi _b(r)\varphi _a(r)\varphi _b(\mathrm{})|^2}{\mathrm{E}ϵ_aϵ_b}}.`$ (22)
Making now use of the symmetry properties of $`\varphi _a`$, i.e., $`\varphi _a^\mathrm{e}(r)=\varphi _a^\mathrm{e}(\mathrm{})`$, and $`\varphi _a^\mathrm{o}(r)=\varphi _a^\mathrm{o}(\mathrm{})`$, we see that the only contributions to the sum come from the cases in which $`a`$ is even, and $`b`$ is odd, and vice versa. For the sake of brevity, in the following expressions we replace $`|\varphi _a^i(\mathrm{})|=|\varphi _a^i(r)|`$ ($`i=e,o`$) by $`|\varphi _a^i|`$. Consequently, we may write the result in the form

$`{\displaystyle \frac{1}{\mathrm{V}\mathrm{K}}}=G_{\mathrm{eo}}(\mathrm{E}),`$ (23)
where $`G_{\mathrm{eo}}`$ is a noninteracting two-electron Green’s function,

$`G_{\mathrm{eo}}(\mathrm{E})`$ $`=`$ $`4{\displaystyle \frac{|\varphi _\mathrm{B}^\mathrm{e}|^2|\varphi _\mathrm{B}^\mathrm{o}|^2}{\mathrm{E}ϵ_\mathrm{B}^\mathrm{e}ϵ_\mathrm{B}^\mathrm{o}}}+2{\displaystyle \underset{\stackrel{i,j=\mathrm{e},\mathrm{o}}{ij}}{}}{\displaystyle \underset{k}{}}{\displaystyle \frac{|\varphi _\mathrm{B}^i|^2|\varphi _k^j|^2}{\mathrm{E}ϵ_\mathrm{B}^iϵ_k}}`$ (24)
$`+`$ $`{\displaystyle \underset{kk^{}}{}}{\displaystyle \frac{|\varphi _k^\mathrm{e}|^2|\varphi _k^{}^\mathrm{o}|^2}{\mathrm{E}ϵ_kϵ_k^{}}}.`$ (25)
(In writing down this equation, we have assumed the existence of two single-particle bound states, as mentioned above.) Here the subscript B denotes a single-particle bound state, and $`k`$ refers to a band state, with $`ϵ_k=2\mathrm{cos}k`$. In the continuum of the band energies, there is no need to distinguish between $`ϵ_k^\mathrm{e}`$ and $`ϵ_k^\mathrm{o}`$. Also, the sum over all states $`k`$ is divided into the sum over the even propagating states, and the sum over the odd ones. We give in the Appendix the explicit expressions for the eigenstates required for the calculation of $`G`$.

Equation (23) is an implicit equation for the exact two-electrons eigenenergies E. We will postpone the discussion of these solutions, and consider now the two-electron states with $`S_z=0`$,

$`|\mathrm{\Psi }^0`$ $`=`$ $`{\displaystyle \underset{ab}{}}\mathrm{X}_{ab}^0c_a^{}c_b^{}|0`$ (26)
$`=`$ $`{\displaystyle \underset{ab}{}}\mathrm{X}_{ab}^0{\displaystyle \underset{nn^{}}{}}\varphi _a^{}(n)\varphi _b^{}(n^{})c_n^{}c_n^{}^{}|0`$ (28)
$`{\displaystyle \underset{nn^{}}{}}\mathrm{X}_{nn^{}}^0c_n^{}c_n^{}^{}|0,`$
with $`_{nn^{}}|\mathrm{X}_{nn^{}}^0|^2=1`$. For the triplet state, the amplitudes $`\mathrm{X}_{nn^{}}^0`$ are antisymmetric in the site indices. Then (inserting Eq. (28) in the Schrödinger equation) the energies are again given by Eqs. (23) and (25). For the singlet states, $`\mathrm{X}_{nn^{}}^0=\mathrm{X}_{n^{}n}^0`$ and the Schrödinger equation yields

$`\mathrm{X}_{nn^{}}^0={\displaystyle \underset{ab}{}}(\mathrm{U}{\displaystyle \underset{i=\mathrm{},r}{}}\mathrm{X}_{ii}^0{\displaystyle \frac{\varphi _a(i)\varphi _b(i)}{\mathrm{E}ϵ_aϵ_b}}+(\mathrm{V}+\mathrm{K})\mathrm{X}_\mathrm{}r^0`$ (29)
$`\times {\displaystyle \frac{\varphi _a(\mathrm{})\varphi _b(r)+\varphi _a(r)\varphi _b(\mathrm{})}{\mathrm{E}ϵ_aϵ_b}})\varphi ^{}_a(n)\varphi ^{}_b(n^{}).`$ (30)
We use this equation for $`n,n^{}=\mathrm{},\mathrm{}`$, $`n,n^{}=r,r`$, and $`n,n^{}=\mathrm{},r`$, and find two families of singlet solutions: (i) $`\mathrm{X}_{\mathrm{}\mathrm{}}^0=\mathrm{X}_{rr}^0`$, $`\mathrm{X}_\mathrm{}r^0=0`$, for which

$`{\displaystyle \frac{1}{\mathrm{U}}}=G_{\mathrm{eo}}(\mathrm{E});`$ (31)
(ii) $`\mathrm{X}_{\mathrm{}\mathrm{}}^0=\mathrm{X}_{rr}^0`$, with

$`1`$ $``$ $`(\mathrm{U}+\mathrm{V}+\mathrm{K})(G_{\mathrm{ee}}(\mathrm{E})+G_{\mathrm{oo}}(\mathrm{E}))`$ (32)
$`+`$ $`4\mathrm{U}(\mathrm{V}+\mathrm{K})G_{\mathrm{ee}}(\mathrm{E})G_{\mathrm{oo}}(\mathrm{E})=0.`$ (33)
\[Note that the last equation includes twice the number of solutions as Eqs. (23) and (31).\] Here, $`G_{\mathrm{ee},\mathrm{oo}}`$ are noninteracting, two-particle Green’s functions, which consist of the even and odd (with respect to interchanging the dots) solutions of the single-particle spectrum, respectively, with

$`G_{\mathrm{ee}(\mathrm{oo})}(\mathrm{E})`$ $`=`$ $`2{\displaystyle \frac{|\varphi _\mathrm{B}^{\mathrm{e}(\mathrm{o})}|^4}{\mathrm{E}2ϵ_\mathrm{B}^{\mathrm{e}(\mathrm{o})}}}+2{\displaystyle \underset{k}{}}{\displaystyle \frac{|\varphi _\mathrm{B}^{\mathrm{e}(\mathrm{o})}|^2|\varphi _k^{\mathrm{e}(\mathrm{o})}|^2}{\mathrm{E}ϵ_\mathrm{B}^{\mathrm{e}(\mathrm{o})}ϵ_k}}`$ (34)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{kk^{}}{}}{\displaystyle \frac{|\varphi _k^{e(o)}|^2|\varphi _k^{}^{\mathrm{e}(\mathrm{o})}|^2}{\mathrm{E}ϵ_kϵ_k^{}}}.`$ (35)
We next determine the ground-state energy, starting with the triplet states, whose energies are given by Eq. (23). We assume that we are in the lower region II of Fig. 1, where there are two bound states below the band. (Note that in the regions marked 0 and I in Fig. 1, where there is at most one single-particle bound state, the triplet ‘bound’ state will always lie in the continuum.) The function $`G_{\mathrm{eo}}(\mathrm{E})`$, Eq. (25), has the following behavior. As E approaches $`\mathrm{}`$, it goes to zero from below. At E$`=ϵ_\mathrm{B}^\mathrm{e}+ϵ_\mathrm{B}^\mathrm{o}`$, it diverges to $`\mathrm{}`$, jumps to $`+\mathrm{}`$ as E crosses that value, and then decreases, as E approaches the bottom of the two-electron continuum states, located at $`2+ϵ_\mathrm{B}^\mathrm{e}`$. As discussed in our earlier work , in the thermodynamic limit of infinite ‘leads’ $`G_{\mathrm{eo}}`$ has a finite value at this band threshold, due to the $`k`$-dependence of $`\varphi _k^{\mathrm{e}(\mathrm{o})}`$ at the impurities. The triplet ground-state energy is where $`G_{\mathrm{eo}}`$ crosses $`(\mathrm{V}\mathrm{K})^1`$. Hence, there will be a two-electron bound state only when $`(\mathrm{V}\mathrm{K})^1<G_{\mathrm{eo}}(2+ϵ_\mathrm{B}^\mathrm{e})`$. It follows that there are values of the direct and exchange Coulomb couplings such that the ground triplet state is not bound, but lies in the continuum.

We now turn to the singlet states, again assuming the existence of both $`ϵ_\mathrm{B}^\mathrm{e}`$ and $`ϵ_\mathrm{B}^\mathrm{o}`$. Consider first the solutions given by Eq. (31). Since U$`>`$V (and V–K), the lowest solution of this equation lies above the lowest solution of the triplet state, (which is given by the same function $`G_{\mathrm{eo}}`$). Hence, we need not consider anymore the states given by (31). To explore the other family of singlet solutions, it is convenient to re-write Eq. (33) in the form

$`{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{1}{G_{\mathrm{ee}}(\mathrm{E})}}+{\displaystyle \frac{1}{G_{\mathrm{oo}}(\mathrm{E})}}\right)={\displaystyle \frac{\mathrm{U}+\mathrm{V}+\mathrm{K}}{2}}`$ (36)
$`\pm \sqrt{{\displaystyle \frac{1}{16}}\left({\displaystyle \frac{1}{G_{\mathrm{ee}}(\mathrm{E})}}{\displaystyle \frac{1}{G_{\mathrm{oo}}(\mathrm{E})}}\right)^2+\left({\displaystyle \frac{\mathrm{U}\mathrm{V}\mathrm{K}}{2}}\right)^2}.`$ (37)
Since the behavior of $`G_{\mathrm{ee}}`$ and $`G_{\mathrm{oo}}`$ as function of E is similar to that of $`G_{\mathrm{eo}}`$ described above, it follows from Eq. (37) that the lowest singlet state energy obeys that equation with the minus sign.

## IV Singlet-triplet transitions

In order to decide when the lowest bound state of the two electrons is a singlet or a triplet, we need to (i) determine for which values of the Coulomb parameters Eqs. (23) and (37) have bound solutions; and (ii) to compare these two solutions, when they exist. Consider as an example the case in which there is only the on-site Hubbard interaction, that is, V=K=0. Then the triplet bound state has the energy $`\mathrm{E}_\mathrm{T}=ϵ_\mathrm{B}^\mathrm{e}+ϵ_\mathrm{B}^\mathrm{o}`$, cf. Eq. (23). The singlet energy, E<sub>S</sub>, in that case is given by the lowest solution of

$`{\displaystyle \frac{1}{\mathrm{U}}}=G_{\mathrm{ee}}(\mathrm{E})+G_{\mathrm{oo}}(\mathrm{E}).`$ (38)
Similarly to the behavior of $`G_{\mathrm{eo}}(\mathrm{E})`$, the right-hand-side of this equation starts at very small negative values when E tends to $`\mathrm{}`$. It then diverges to $`\mathrm{}`$ as E approaches $`2ϵ_\mathrm{B}^\mathrm{e}`$, jumps to $`+\mathrm{}`$ as E crosses that value, diverges again to $`\mathrm{}`$ as E$`2ϵ_\mathrm{B}^\mathrm{o}`$, then jumps to $`+\mathrm{}`$, and finally decreases towards a finite value as E approaches the bottom of the two-electron continuum. It follows that Eq. (38) has always a bound energy solution. Moreover, if $`G_{\mathrm{ee}}+G_{\mathrm{oo}}`$ is negative at E$`=\mathrm{E}_\mathrm{T}ϵ_\mathrm{B}^\mathrm{e}+ϵ_\mathrm{B}^\mathrm{o}`$, then that solution E<sub>S</sub> lies below E<sub>T</sub>, i.e., the ground-state is a singlet. This is indeed the case, as is shown in the Appendix \[Eq. (A17)\]. This is in accordance with the general rule, which states that in order for the ground-state to be a triplet, the exchange Coulomb energy has to overcome the kinetic energy.

For the sake of clarity of the presentation, we will carry the rest of the analysis to lowest order in the coupling to the leads, $`\gamma `$. In that case, it is possible to derive simple expressions for the two-particle Green’s functions, see Eqs. (A20). Using those equations, we find that the singlet energies are given by

$`\mathrm{E}2ϵ_0+\gamma (e^{\alpha ^\mathrm{e}}+e^{\alpha ^\mathrm{o}})={\displaystyle \frac{\mathrm{U}+\mathrm{V}+\mathrm{K}}{2}}`$ (39)
$``$ $`\sqrt{4t_D^2+\left({\displaystyle \frac{\mathrm{U}\mathrm{V}\mathrm{K}}{2}}\right)^2+4t_D\gamma (e^{\alpha ^\mathrm{e}}e^{\alpha ^\mathrm{o}})},`$ (40)
and the triplet energies are given by

$`\mathrm{E}2ϵ_0+\gamma (e^{\alpha ^\mathrm{e}}+e^{\alpha ^\mathrm{o}})=\mathrm{V}\mathrm{K},`$ (41)
where $`\alpha ^{\mathrm{e},\mathrm{o}}`$ is related to the corresponding $`\mathrm{E}`$ (with indices S or T) via

$`\mathrm{E}ϵ_\mathrm{B}^{\mathrm{e},\mathrm{o}}=2\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\alpha ^{\mathrm{e},\mathrm{o}}.`$ (42)
Let us examine the case in which in the absence of the coupling to the leads, the singlet and the triplet ground-state energies are equal, i.e., $`4t_D^2=2`$K(U–V+K), see Eq. (14). Then

$`\mathrm{E}_\mathrm{S}\mathrm{E}_\mathrm{T}`$ $`=`$ $`\gamma (e^{\alpha _\mathrm{T}^\mathrm{e}}+e^{\alpha _\mathrm{T}^\mathrm{o}}e^{\alpha _\mathrm{S}^\mathrm{e}}e^{\alpha _\mathrm{S}^\mathrm{o}})`$ (43)
$``$ $`{\displaystyle \frac{4t_D}{\mathrm{U}\mathrm{V}+3\mathrm{K}}}\gamma (e^{\alpha _\mathrm{S}^\mathrm{e}}e^{\alpha _\mathrm{S}^\mathrm{o}}).`$ (44)
The last term on the right-hand-side of this equation is negative. This follows from Eq. (42) and the fact that $`ϵ_\mathrm{B}^\mathrm{e}<ϵ_\mathrm{B}^\mathrm{o}`$. As for the first term, we use again Eq. (42), to write it in the form $`(\mathrm{E}_\mathrm{T}\mathrm{E}_\mathrm{S})\gamma ((e^{2\alpha _\mathrm{S}^\mathrm{e}}1)^1+(e^{2\alpha _\mathrm{S}^\mathrm{o}}1)^1)`$. Hence, E$`{}_{\mathrm{S}}{}^{}<\mathrm{E}_\mathrm{T}`$ and the singlet is prefered.

The above discussion shows that the coupling to the continuum of propagating states enhances the tendency of the two electrons to form a singlet state, in the situation where in the absence of that coupling, the singlet and the triplet states are degenerate. In order to investigate whether this tendency persists for other choices of parameters (and at the same time to keep the calculations tractable) we will now confine ourselves to the choice U=V+K. In this case, again to leading order in $`\gamma `$, the equation for the singlet energies (40) reads

$`f_\mathrm{S}(\mathrm{E})`$ $`=`$ $`\mathrm{V}+\mathrm{K},`$ (45)
$`f_\mathrm{S}(\mathrm{E})`$ $`=`$ $`\mathrm{E}2ϵ_0+2t_D+2\gamma e^{\alpha ^\mathrm{e}}.`$ (46)
Similarly, Eq. (41) for the triplet energies can be written as

$`f_\mathrm{T}(\mathrm{E})`$ $`=`$ $`\mathrm{V}\mathrm{K},`$ (47)
$`f_\mathrm{T}(\mathrm{E})`$ $`=`$ $`\mathrm{E}2ϵ_0+\gamma (e^{\alpha ^\mathrm{e}}+e^{\alpha ^o}).`$ (48)
Let us first determine for which parameters these equations yield bound, two-electron energies. To this end, we consider Eqs. (46) and (48) at the bottom of the two-electron continuum, E=$`2+ϵ_\mathrm{B}^\mathrm{e}`$. The first of these equations will have a bound state for $`f_\mathrm{S}(2+ϵ_\mathrm{B}^\mathrm{e})>\mathrm{V}+\mathrm{K}`$; the second will have such a solution when $`f_\mathrm{T}(2+ϵ_\mathrm{B}^\mathrm{e})>\mathrm{V}\mathrm{K}`$. These conditions are plotted in Fig. 2 as the thick lines there. A bound triplet state exists in regions I+III, below the heavy line of positive slope. A singlet bound state exists in I+II, below the heavy line of negative slope. In region IV, there are no bound states; both the triplet and the singlet states are in the continuum, and their energy is about the same. Crossing the line between regions II and IV (or III and IV) thus corresponds to the delocalization transition discussed above, from a singlet (triplet) bound ground state to a degenerate ‘metallic’ state. This transition is the most striking effect of the coupling to the leads, and we expect it to appear irrespective of the quantitative approximation used in Fig. 2.

In region I one has to compare the singlet energy with the triplet one. These two become degenerate along the diamond curve, whose equation is derived from (42), (46), and (48)

$`\mathrm{V}\mathrm{K}`$ $`=`$ $`{\displaystyle \frac{ϵ_\mathrm{B}^\mathrm{o}}{2}}(\mathrm{B}1)+{\displaystyle \frac{ϵ_\mathrm{B}^\mathrm{e}}{2}}(\mathrm{B}+1)2ϵ_0(1\gamma )\mathrm{AB},`$ (49)
$`\mathrm{A}`$ $`=`$ $`{\displaystyle \frac{2}{\gamma }}(\mathrm{K}t_D),`$ (50)
$`\mathrm{B}`$ $`=`$ $`\sqrt{1+4/[\mathrm{A}(\mathrm{A}+ϵ_\mathrm{B}^\mathrm{o}ϵ_\mathrm{B}^\mathrm{e})]}.`$ (51)
This line is almost vertical, with $`\mathrm{E}_\mathrm{S}<\mathrm{E}_\mathrm{T}`$ to its left, and $`\mathrm{E}_\mathrm{S}>\mathrm{E}_\mathrm{T}`$ to its right. The conclusion is that, as long as there exists a bound state of the two electrons, then this line moves slightly to the right as $`\gamma `$ is increased from zero (when the line was at K$`=t_D`$). However, the coupling to the continuum states delocalizes the electrons, making the two states degenerate over a large part of the parameter plane $`\{\mathrm{V}\mathrm{K}\}`$, i. e. region IV in Fig. 2. With all other parameters fixed, one might expect that increasing the distance between the dots causes a decrease in $`t_D`$, in V and in K, thus causing a shift towards to lower left side of Fig. 2, towards a bound singlet ground state.

## V discussion

We have derived analytical expressions for the spectrum of two interacting electrons on a simplified model for a double quantum-dot. When the dot is decoupled from the external leads, it is straightforward to obtain this spectrum, and discuss the criterion for its ground-state to be a singlet or a triplet. The question we have addressed is how this criterion is modified when the single-particle states which are localized on the dot are coupled to the continuum of extended states on the leads. A typical example of our results is shown in Fig. 2: as long as the the electron-electron interactions do not delocalize the ground-state, then the location for the singlet-triplet transition shifts continuously with the coupling to the leads, $`\gamma `$. In that case, one can still say that the singlet state is the ground one provided that the kinetic energy dominates over the exchange energy. In a way, the coupling to the leads enhances the kinetic energy, and therefore it slightly favors the singlet ground state. However, the coupling to the leads has a much more drastic effect: over a significant part of the parameter space (which consists of the Coulomb couplings, the single-particle energies on the dots, and the coupling to the leads), e. g. region IV in Fig. 2, the interplay between the coupling to the leads and the electron-electron interactions delocalizes one or both electrons. Then, the singlet and the triplet states are degenerate. In such a situation, the bound state (which may be either a singlet or a triplet) disappears and the ground state becomes a degenerate singlet-triplet ‘metallic’ state. We believe that this delocalization effect should be taken into consideration in the analyzes of experimental data related to this question.

###### Acknowledgements.

We enjoyed many discussions with Yoseph Imry. This project was supported by grants from the Israel Science Foundation and from the Israeli Ministry of Science and the French Ministry of Research and Technology, AFIRST.

## A The two-electron Green’s functions

We first list the eigenfunctions required for the calculation of the noninteracting, two-particle Green’s functions. These are needed only on the double quantum dot, that is, on the sites $`r`$ and $`\mathrm{}`$. Writing the even and odd bound-state energies, Eqs. (16), in the form

$`ϵ_\mathrm{B}^{\mathrm{e},\mathrm{o}}=2\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\kappa ^{\mathrm{e},\mathrm{o}},`$ (A1)
the wave functions on the dot sites is

$`|\varphi _\mathrm{B}^{\mathrm{e},\mathrm{o}}|^2={\displaystyle \frac{1}{2}}\left(1+{\displaystyle \frac{\gamma }{e^{2\kappa ^{\mathrm{e},\mathrm{o}}}1}}\right)^1.`$ (A2)
The band states have been calculated assuming periodic boundary conditions for a system of N sites (these include the leads). Then

$`|\varphi _k^{\mathrm{e},\mathrm{o}}|^2={\displaystyle \frac{2}{\mathrm{N}}}{\displaystyle \frac{\gamma \mathrm{sin}^2k}{\gamma ^2\mathrm{sin}^2k+(ϵ_kϵ_0^{\mathrm{e},\mathrm{o}}+\gamma \mathrm{cos}k)^2}},`$ (A3)
with $`ϵ_k=2\mathrm{cos}k`$, for both the even and odd states. In calculating the sums over $`k`$ in the two-particle Green’s functions, we shall use the continuum limit, dividing the $`k`$-integrations on the even (odd) functions by 2.

We next derive $`G_{\mathrm{ee}}`$, $`G_{\mathrm{oo}}`$, and $`G_{\mathrm{eo}}`$. It is convenient first to calculate the function

$`Q^{\mathrm{e},\mathrm{o}}(\omega )={\displaystyle \underset{k}{}}{\displaystyle \frac{|\varphi _k^{\mathrm{e},\mathrm{o}}|^2}{\omega ϵ_k}},\omega <2.`$ (A4)
Writing

$`\omega =2\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\alpha ,`$ (A5)
and using Eq. (A3), we find

$`Q^{\mathrm{e},\mathrm{o}}`$ $`=`$ $`\gamma {\displaystyle \frac{e^{\kappa ^{\mathrm{e},\mathrm{o}}}}{e^{2\kappa ^{\mathrm{e},\mathrm{o}}}1+\gamma }}`$ (A7)
$`\times {\displaystyle \frac{e^{\alpha +\kappa ^{\mathrm{e},\mathrm{o}}}}{(e^{\alpha +\kappa ^{\mathrm{e},\mathrm{o}}}1+\gamma )(e^{\alpha +\kappa ^{\mathrm{e},\mathrm{o}}}1)}},`$
where we have also used

$`{\displaystyle \underset{k}{}}|\varphi _k^{\mathrm{e},\mathrm{o}}|^2=12|\varphi _\mathrm{B}^{\mathrm{e},\mathrm{o}}|^2.`$ (A8)
Exploiting this result, we now introduce the function $`F`$,

$`F^{\mathrm{e},\mathrm{o}}(\omega )`$ $`=`$ $`|\varphi _\mathrm{B}^{\mathrm{e},\mathrm{o}}|^2+(\omega ϵ_\mathrm{B}^{\mathrm{e},\mathrm{o}}){\displaystyle \frac{1}{2}}Q^{\mathrm{e},\mathrm{o}}(\omega )`$ (A9)
$`=`$ $`{\displaystyle \frac{1}{2}}\times {\displaystyle \frac{1}{1+(\gamma /(e^{\alpha +\kappa ^{\mathrm{e},\mathrm{o}}}1))}}.`$ (A10)
Note that

$`F^\mathrm{e}(ϵ_\mathrm{B}^\mathrm{e})`$ $`=`$ $`|\varphi _\mathrm{B}^\mathrm{e}|^2,F^\mathrm{o}(ϵ_\mathrm{B}^\mathrm{o})=|\varphi _\mathrm{B}^\mathrm{o}|^2,`$ (A11)
$`F^\mathrm{e}(ϵ_\mathrm{B}^\mathrm{o})=F^\mathrm{o}(ϵ_\mathrm{B}^\mathrm{e})`$ $`=`$ $`{\displaystyle \frac{1}{2}}\times {\displaystyle \frac{1}{1+(\gamma /(e^{\kappa ^\mathrm{e}+\kappa ^\mathrm{o}}1))}}.`$ (A12)
It is now straightforward to show, using Eqs. (35), (25), and (A10) that

$`G_{\mathrm{ee}}(\mathrm{E})=2{\displaystyle \frac{(F^\mathrm{e}(\mathrm{E}ϵ_\mathrm{B}^\mathrm{e}))^2}{\mathrm{E}2ϵ_\mathrm{B}^\mathrm{e}}}`$ (A13)
$`+`$ $`{\displaystyle \underset{kk^{}}{}}{\displaystyle \frac{|\varphi _k^\mathrm{e}|^2|\varphi _k^{}^\mathrm{e}|^2(ϵ_\mathrm{B}^\mathrm{e}ϵ_k)(ϵ_\mathrm{B}^\mathrm{e}ϵ_k^{})/2}{(\mathrm{E}ϵ_kϵ_k^{})(\mathrm{E}ϵ_\mathrm{B}^\mathrm{e}ϵ_k)(\mathrm{E}ϵ_\mathrm{B}^\mathrm{e}ϵ_k^{})}},`$ (A14)
with an analogous result for $`G_{\mathrm{oo}}`$, with e replaced by o, and

$`G_{\mathrm{eo}}(\mathrm{E})=4{\displaystyle \frac{F^\mathrm{e}(\mathrm{E}ϵ_\mathrm{B}^\mathrm{o})F^\mathrm{o}(\mathrm{E}ϵ_\mathrm{B}^\mathrm{e})}{\mathrm{E}ϵ_\mathrm{B}^\mathrm{e}ϵ_\mathrm{B}^\mathrm{o}}}`$ (A15)
$`+`$ $`{\displaystyle \underset{kk^{}}{}}{\displaystyle \frac{|\varphi _k^\mathrm{o}|^2|\varphi _k^{}^\mathrm{e}|^2(ϵ_\mathrm{B}^\mathrm{o}ϵ_k)(ϵ_\mathrm{B}^\mathrm{e}ϵ_k^{})}{(\mathrm{E}ϵ_kϵ_k^{})(\mathrm{E}ϵ_\mathrm{B}^\mathrm{e}ϵ_k)(\mathrm{E}ϵ_\mathrm{B}^\mathrm{o}ϵ_k^{})}}.`$ (A16)
For energies E below the bottom of the two-electron continuum, i.e., E$`<2+ϵ_\mathrm{B}^\mathrm{e}`$, the double sum on $`k`$ and $`k^{}`$ in (A14) is negative. Using Eq. (A12) for the first terms in the equations for $`G_{\mathrm{ee}}`$ and $`G_{\mathrm{oo}}`$, it follows that

$`G_{\mathrm{ee}}(ϵ_\mathrm{B}^\mathrm{e}+ϵ_\mathrm{B}^\mathrm{o})+G_{\mathrm{oo}}(ϵ_\mathrm{B}^\mathrm{e}+ϵ_\mathrm{B}^\mathrm{o})<0.`$ (A17)
This result is used to show that when the only Coulomb coupling is the Hubbard U, the singlet is always the ground state.

Up to this point, the results were given for general $`\gamma `$. To lowest order in the coupling to the leads, $`\gamma `$, we may discard the double sums in (A14) and (A16). Then, using (A10), we find

$`{\displaystyle \frac{1}{G_{\mathrm{eo}}(\mathrm{E})}}`$ $`\mathrm{E}2ϵ_0+\gamma (e^{\alpha ^\mathrm{e}}+e^{\alpha ^\mathrm{o}}),`$ (A18)
$`{\displaystyle \frac{1}{2G_{\mathrm{ee}}(\mathrm{E})}}`$ $`\mathrm{E}2ϵ_0+2t_D+2\gamma e^{\alpha ^\mathrm{e}},`$ (A19)
$`{\displaystyle \frac{1}{2G_{\mathrm{oo}}(\mathrm{E})}}`$ $`\mathrm{E}2ϵ_02t_D+2\gamma e^{\alpha ^\mathrm{o}},`$ (A20)
where we have used Eq. (42).


format: text

text:
# Dust-enshrouded Asymptotic Giant Branch Stars in the Solar Neighbourhood

## 1 Introduction

Stars of low to intermediate main-sequence mass, up to about 6-8$`M_{}`$, will eventually evolve up the Asymptotic Giant Branch (AGB) (Iben & Renzini 1983). At the tip of the AGB, where they combine a low effective temperatures ($`T_{eff}\stackrel{<}{}3000K`$) with high luminosities ($`L\stackrel{>}{}3000L_{}`$), the stars pulsate and are subject to rapid mass loss. Typical examples of such stars are the Miras and the OH/IR111A small fraction of OH/IR stars are supergiants, i.e. their progenitors were massive stars. stars (e.g. Feast & Whitelock 1987; Habing 1996). Their light output varies on timescales from around one hundred days to well over a thousand days. They also lose mass rapidly ($`10^7\stackrel{<}{}\dot{M}\stackrel{<}{}10^4M_{}yr^1`$). They follow various period-luminosity relations and are at or near the maximum luminosity they will ever achieve during their evolution. This makes them very useful as distance indicators and as probes of galactic structure. They are also important contributors of dust into the interstellar medium, which may be crucial for the formation of new stars and planets.

Jura & Kleinmann (1989 henceforth JK89) identified a group of mass losing AGB stars within about 1 kpc of the sun. In 1989 a monitoring programme of stars in the JK89 sample visible from the southern hemisphere was initiated from SAAO. This entailed obtaining near-infrared photometry at $`JHKL`$ over a time interval of several years, from which pulsation periods could be determined. The aim of this study was to investigate the physical properties of the sample, in particular variability, colours, mass-loss rates and kinematics, using SAAO photometry together with infrared photometry and kinematic data from the literature.

Finally these data are used to estimate the average progenitor mass and lifetime of stars in this high mass-loss phase.

## 2 The Sample

JK89 aimed to develop a more quantitative understanding of mass loss from local AGB stars. They identified AGB stars undergoing heavy mass loss ($`\dot{M}>10^6M_{}yr^1`$), henceforth referred to as dust-enshrouded, within about 1 kpc of the sun, and within $`33^{}<\delta (1950)<82^{}`$. Their selection was based on infrared photometry from the Two Micron Sky Survey, IRAS and the AFGL survey from which they identify 63 sources. They give reasons for thinking this sample was reasonably complete.

A few comments should be made at this point. JK89 assume a luminosity of $`10^4L_{}`$ for all of these AGB stars. It will be seen later that, if the period luminosity relation (PL) is correct, this leads to the inclusion of more luminous OH/IR stars outside the 1 kpc sphere and up to $`2`$ kpc away. They note that some red supergiants will satisfy their selection criteria; those they largely excluded, but included the star 3 Puppis, an A supergiant, with the idea that it might be undergoing the transition from AGB star to planetary nebula. However, Plets, Waelkens & Trams (1995) convincingly showed that 3 Puppis is a genuine massive supergiant, and we therefore exclude it from this study. JK89 also included the carbon-rich planetary nebulae NGC 7027 (e.g. Justtanont et al. 2000) in their sample; this object is excluded from the present study. The remaining 61 sources are listed in Table 1 (this includes 3 M supergiants, $`175132313`$, $`181351641`$ and $`182041344`$ (Ukita & Goldsmith 1984; Winfrey et al. 1994) of which we only became aware part the way through this project). Near-infrared data are tabulated for two of them, but none are used in the analysis.

Chemical types (O and C for oxygen- and carbon-rich, respectively) for all the sources were taken from Loup et al. (1993), and spectral types obtained from the SIMBAD data base for 51 of the sources are also listed in Table 1. Although $`04307+6210`$ is listed as M6 in SIMBAD, Groenewegen (1994) suggests this is incorrect and the star is actually carbon rich; therefore no spectral type is assigned to it here. The supergiants are marked (I).

The final sample thus consists of 27 carbon-rich stars and 31 oxygen-rich stars (excluding the 3 supergiants) of which 2 are S-type. Almost all of the stars are of very late spectral type (M, C, N, R and S), with two exceptions: the Red Rectangle and the Egg Nebula. Both are well-known carbon-rich bipolar reflection nebulae (see, e.g. Sahai et al. 1998; Waters et al. 1998). The Red Rectangle is a wide binary system with an oxygen-rich dust disk (Waters et al. 1998). These two objects are generally thought to be post-AGB stars on their way to become Planetary Nebulae.

The Benson et al. (1990) catalogue lists sources which have been examined for OH, SiO and $`\mathrm{H}_2\mathrm{O}`$ maser emission. Of the 28 oxygen-rich stars in our sample which have been examined for OH and $`\mathrm{H}_2\mathrm{O}`$, 20 were detected in OH and 20 in $`\mathrm{H}_2\mathrm{O}`$. While 19 of the 21 sources searched for SiO were detected. Interestingly, 4 of the 11 carbon-rich stars which were searched for SiO also had positive detections. We can thus estimate that 70 percent of the oxygen-rich stars have OH and $`\mathrm{H}_2\mathrm{O}`$ maser emission and about 90 percent have SiO emission.

From the number of references listed in Tables 1, CW Leo stands out as the best studied source in this sample. About 50 molecular species have been detected in its envelope (Wallerstein & Knapp 1998) which is more than for any other source. The full sample consists mainly of cool variable stars with a few possible post-AGB objects.

## 3 Observations

### 3.1 Near-infrared Photometry

Near-infrared observations were obtained in the $`JHKL`$ bands (1.2, 1.65, 2.2 and 3.45$`\mu `$m) for 42 of the 45 southern sources, defined henceforth as stars with $`\delta <+30^{}`$. Some of these 42 stars were selected exclusively from JK89, while the others were already part of other monitoring programmes. The remaining three southern stars ($`06500+0829`$, $`175132313`$ and RAFGL 1406) were not observed because they were in crowded fields and confusion rendered observations with a single channel photometer impossible (note that Winfrey et al. (1994) suggest that $`175132313`$ is actually a supergiant and it is therefore not referred to again).

The near-infrared observations for the 42 southern sources, accurate to $`\pm 0.1`$ mag and better (0.03 mag in $`JHK`$ and 0.05 mag in $`L`$ in the majority of cases), are listed in Table 2. The magnitudes from the 0.75 m reflector are on the SAAO system as defined by Carter (1990). The $`HKL`$ bands for the 1.9 m reflector are assumed to be identical to those of the 0.75 m reflector, while the $`J`$-magnitude is converted to the SAAO system.

In some instances, no observations are listed because the sources were too faint to be reliably measured through the relevant filter with a particular telescope. The times of the observations are given as Julian dates from which 244 0000 days has been subtracted. In addition to the data from Table 2, published SAAO observations from Whitelock et al. (1994, 1997a) and Lloyd Evans (1997) have been used in the analysis, together with several unpublished observations from Lloyd Evans (private communication) with a sum total of 848 measurements.

Mean magnitudes calculated from the average of observed values at maximum and minimum light are listed in Table 3 and compared with the Fourier mean in appendix A. These two averages compare fairly well, differences were less than 0.3 mag for the $`J`$-band and less than about 0.1 mag for the $`HKL`$ bands. Table 3 also lists mean $`K`$-magnitudes calculated as above from values listed in the literature, for the 20 remaining stars with no near-infrared SAAO photometry.

### 3.2 Mid- and Far-infrared Observations

#### 3.2.1 Photometry

Infrared fluxes at 12, 25 and $`60\mu `$m were taken from the IRAS Point Source Catalog (PSC) and colour corrected as specified in the IRAS Explanatory Supplement (ES). The results are listed in Table 4 together with the IRAS flux-quality designations and the IRAS variability index (V.I.), which gives the probability that the source is variable in units of 0.1 (see IRAS ES). Table 4 also contains $`11\mu `$m magnitudes from the *Revised Airforce Geophysics Laboratory, RAFGL*, survey.

#### 3.2.2 Low Resolution Spectra (LRS)

The IRAS Low Resolution spectral classifications, from Olnon & Raimond (1986) or Loup et al. (1993), are listed in Table 4. Note that $`181351641`$ does not have a published classification so we examined the spectrum in the IRAS Low Resolution Spectra electronic database at the University of Calgary. It shows a weak $`10\mu `$m emission feature, but no $`18\mu `$m feature, which is present in the spectra of all the other oxygen-rich stars. It is actually a supergiant as discussed below (section 4).

## 4 Variability

Using the data in Table 2, the southern sources were classified as variable in the near-infrared if the difference between the maximum and minimum magnitudes at $`K`$ and/or at $`L`$ was greater than 0.1 mag; thirty seven of the stars are variable and 2 non-variable. The only two decisively non-variable stars, $`181351641`$ and $`182041344`$, are the supergiants (see section 2). There were insufficient observations to determine if the other 6 stars are variable or not. Table 5 lists the variability information.

Period determinations, using Fourier transforms of the $`K`$ light curve, were made for all variables with measurements on at least 10 epochs. Plausible, although not definitive, periods were also established for some sources with as few as 8 measurements. For the very red source $`174113154`$, a period could only be estimated from the $`L`$ data. Periods, listed in Table 5, were thus determined for a total of 20 stars, and phased $`K`$ light curves for these are shown in Fig. 1. Where possible the Fourier transform analysis was carried out at the other wavelengths and in all cases the periods found were within 2 percent of those listed in Table 5. No period determination was possible for $`191260708`$ (W Aql), although it had been observed on 12 epochs. This is because 6 of the epochs were between 1968 and 1982 when this carbon star was undergoing a dust obscuration episode (Mattei & Foster 1997). Such an event would disrupt the apparent periodicity.

The average of periods found in the literature are listed in Table 5, with references, for the 47 stars concerned. Where there are periods in common none of them differs significantly from those determined here. In the following analysis our own determinations, where available, are used in preference to others, except for $`173603012`$ and $`183480526`$ where our determination is uncertain due to the small number of measurements available.

### 4.1 Phased Light-curves

Using the periods found for these 20 stars, the relative data were phased, arbitrarily assuming zero phase at JD 2440000.00, where JD is the heliocentric Julian date. These phased light-curves are plotted in Fig. 1 together with least squares fits to the data of the form:

$$M=\frac{\mathrm{\Delta }M}{2}\mathrm{sin}(\varphi \varphi _0)+M_{av}$$
(1)
where M is the magnitude in a particular band, $`\mathrm{\Delta }M`$ the peak to peak amplitude, $`\varphi `$ related to the phase, $`p`$, by $`\varphi =2\pi p`$ and $`\varphi _0`$ the phase angle at zero phase. $`M_{av}`$ is the mean magnitude.

One immediate observation from most of the phased light-curves is that at a given phase there is a range of magnitudes, rather than one sharply defined magnitude. This is not, in general, due to observational error, but is contributed to by various effects. First, Miras and OH/IR stars are erratic pulsators in the sense that their light curves do not exactly repeat from one cycle to the next at visual or infrared wavelengths (Koen & Lombard 1995; Whitelock, Marang & Feast 2000) Secondly, some stars, particularly carbon stars go through faint phases from time to time due to variable dust obscuration, e.g. $`022702619`$ (R For, see Whitelock et al 1997a). Thirdly, a very small number of stars show temporal or long term changes in luminosity, amplitude and/or period, possibly due to evolutionary effects (Whitelock 1999 and references therein).

Many of the light curves in Fig. 1 show clear departures from simple sinusoids and would be better fitted by the inclusion of one or two harmonics, in particular, $`094292148`$ and $`12447+0425`$. The most strikingly odd curve is that of $`104912059`$ (V Hya) which is discussed below.

#### 4.1.1 V Hya

V Hya is a well studied carbon star and semi-regular variable. It has a bipolar circumstellar outflow, and models invoke rotation of the central star or the presence of a binary companion to explain bipolar outflows (Barnbaum, Morris & Kahane 1995).

Fig. 2a illustrates the light curve for V Hya using $`K`$ magnitudes measured from SAAO, including material from Lloyd Evans (1997 and private communication) together with the data from Table 2. In addition to the regular pulsation V Hya experienced two faint phases around JD 244 3500 and 244 9500. A Fourier analysis of the $`K`$ data when the star is bright, between JD 244 6000 and 244 8250 provide a period of 530 days, identical to the value given in the GCVS. Fig. 2b shows the data phased at this period and illustrates the double peak caused by harmonics in the light curve.

Knapp et al. (1999) analysed 3260 AAVSO visual magnitude estimates covering almost 35 years. They confirm and refine the very long period found by Mayall (1965) at $`6160\pm 400`$ days. The two faint phases seen in Fig. 2a are consistent with this long period. The amplitudes of these 6160 day variations are at least 2.4, 2.1, 1.7 and 0.7 mag, at $`JHKL`$, respectively, while $`\mathrm{\Delta }V3.5`$ mag (Knapp et al. 1999).

Barnbaum et al. (1995) found evidence for rapid rotation in V Hya, using spectral line broadening analysis. They suggest that the long term variation can be attributed to a coupling between the radial pulsation, with period of 530 days, and the rotation of the star with a very similar period. It is not clear how this interaction would give rise to the very long period without strongly modulating the shorter period, or why it should result in the observed reddening of the star during the faint phases.

Lloyd Evans (1997) suggested that the long period may be the result of obscuration events like those observed in R For. Although the colours would be consistent with this, the regularity of the long period (Knapp et al. 1999) is quite unlike the effect observed in R For and this explanation therefore seems unlikely.

Knapp et al. (1999) on the other hand, suggest that the long-period dimming of this star is due to a thick dust cloud orbiting the star and attached to a binary companion. Indeed, there is evidence for large numbers of late-type binaries in the LMC (Wood et al. 1999), so this may not be rare occurrence. The colour variations of V Hya are less extreme than those expected from optically thin interstellar dust, but they could be explained by large particles.

### 4.2 Period Distributions

The period distributions for the oxygen- and carbon-rich stars are plotted in Fig. 3. Most of the stars in the sample have periods between 500 and 700 days and almost all have periods between 300 and 900 days. There are only three extreme OH/IR stars with periods in excess of 1000 days. Fig. 3 can be compared with Whitelock et al. (1994) fig. 8 which shows the period distribution of IRAS selected Miras in the South Galactic Cap. Most of the Cap sample have $`150<P<450`$ days, i.e. shorter than the sample under discussion, but longer than optically selected Miras. Longer period Miras have more massive progenitors and are therefore lower scale heights. It is not surprising to find more of them in the solar neighbourhood than in the Galactic Cap.

A comparison of the period distributions of the oxygen- and carbon-rich stars is instructive. There are no carbon-rich stars with periods longer than 1000 days. From Fig. 3 it can be seen that the maximum of the carbon-star period distribution seems to be slightly shifted to longer periods than that of the oxygen-rich stars. A Kolmogorov-Smirnov (KS) test, when applied to these two distributions for periods less than 1000 days, yields a probability that the two distributions come from the same parent population of $`0.40`$. This probability ($`>0.05`$) implies that there is no significant difference in the period distributions for the carbon- and oxygen-rich stars.

The KS test cannot be used if the very long period stars are included since it is not very sensitive to differences in the tails of the distributions, as is well known. In this case the F-test and then the appropriate version of the student t-test is applied. The F-test gives a probability that the variances of the two distributions are the same of $`2.5\times 10^5`$, i.e. very unlikely, which is obvious from Fig. 3. The appropriate version of the t-test, which takes differences in variances into account, gives a probability that the two distributions have the same mean of $`0.24`$. This probability ($`>0.05`$) implies that there is no statistical evidence for any difference in the mean of the two distributions.

### 4.3 Amplitudes

Table 5 lists the stars for which periods were determined together with their amplitudes, from equation 1, at $`JHKL`$. These amplitudes are plotted against period in Fig. 4, where a clear correlation can be seen at all wavelengths.

The amplitude clearly decreases with wavelength (see also Feast et al. 1982). Quantitatively, the average amplitudes, relative to that at $`J`$ are, for $`H`$, $`K`$ and $`L`$ respectively, $`0.84\pm 0.1`$, $`0.69\pm 0.14`$ and $`0.59\pm 0.14`$. The linear correlation coefficient between the average relative amplitude and wavelength is $`0.93`$ with a correlation probability of $`Prob=0.07`$ (using the relative amplitude at $`J`$ of $`1`$). This probability is considered as marginally significant. The KS test gives the probability for the null hypothesis that the light amplitude distributions for the carbon- and oxygen-rich stars (with $`P<1000`$ days) came from the same parent population at $`J`$, $`H`$, $`K`$ and $`L`$ of 0.25, 0.92, 0.51 and 0.52 respectively. Hence there is no statistical reason to suspect that these carbon- and oxygen-rich stars have significantly different light-amplitudes in the near-infrared.

## 5 Colours

### 5.1 Average Colours

Average near-infrared colours, i.e. consisting of $`JHKL`$ magnitudes only, were calculated for each star using the magnitudes at the dates of minimum and maximum $`K`$-light. This average is compared to a Fourier-fit colour in Appendix A. The two values are quite comparable with mean differences for all the near-infrared colours not exceeding 0.1 mag. Note that due to the very small phase coverage of $`174113154`$ in the $`K`$-band, the derived $`KL`$ colour of this star may not be a good representative of the average $`KL`$ colour.

The IRAS colours were obtained using the data in Table 4, where the colour-corrected IRAS fluxes were converted to magnitudes using the zero points listed in the IRAS-ES. Colours consisting of $`JHKL`$ and IRAS photometry, were calculated by combining average magnitudes in the $`JHKL`$ bands (Table 3) with the IRAS magnitudes.

For the northern sources and the three non-IRAS sources the average $`K`$ magnitudes were taken from Table 3, to calculate the $`K[12]`$ colour. For the non-IRAS (carbon-rich) sources the average $`12\mu `$m magnitudes were estimated from the average $`[11][12]=0.1\pm 0.3`$ colour for the carbon-rich sources and the 11$`\mu `$m magnitudes listed in Table 4. All the derived infrared colours are given in Table 6.

### 5.2 Colour-Colour Diagrams

Infrared two colour diagrams are illustrated in Figs. 5 to 8. There is a clear separation of the oxygen- and carbon-rich stars in these diagrams with the carbon stars falling to the lower right of the oxygen-rich stars in the near-infrared two colour diagrams (Figs. 5 and 6) and to the left of them in the combined diagram, Fig. 7. The $`JHKL`$ colours of these AGB stars are more extreme than those of the well studied Miras (e.g. Feast et al. 1982), or even of IRAS selected Miras in the Galactic Cap (Whitelock et al. 1994).

#### 5.2.1 The $`K[12]`$ versus $`[12][25]`$ Diagram

Of particular interest is the $`K[12]`$ vs. $`[12][25]`$ diagram (Fig. 7). The colours of the oxygen- and carbon-rich stars separated rather clearly in this diagram. This has been noted previously (Le Bertre et al. 1994) and shown to be the result of the difference in the ratios of the near-infrared to mid-infrared opacities for the carbon- and oxygen-rich dust.

The single carbon star which falls among the oxygen-rich stars is $`061761026`$ (the Red Rectangle), a well studied, but peculiar star. It has an IRAS LRS classification of $`80`$ indicating the presence of the $`11.3\mu `$m line in emission in contrast to the other carbon stars which show $`11.3\mu `$m absorption, hence its position in Fig. 7. It is thought to be a pre-planetary nebula with an oxygen-rich circumstellar disk (see section 2.1).

The dashed line in Fig. 7 has been positioned so as to separate the carbon- and oxygen-rich stars as far as possible, it is described by:

$$(K[12])=11.5([12][25])3.8$$
(2)
Thus if we define $`\delta _{C/O}`$ as follows:

$$\delta _{C/O}=11.5([12][25])(K[12])3.8,$$
(3)
The two chemical types can be separated by a single parameter, $`\delta _{C/O}`$, such that for carbon stars $`\delta _{C/O}<0`$ and for most oxygen-rich stars $`\delta _{C/O}>0`$. This parameter might be useful for separating these objects over the colour range discussed. On this sample, obviously the one used to define the parameter, only 6 percent of stars were misclassified.

The two S-type stars were omitted from the above discussion, we note that one of them falls with the carbon- and the other with the oxygen-rich stars.

#### 5.2.2 IRAS Two Colour Diagram

Van der Veen & Habing (1988) demonstrated that the position of an AGB star in an IRAS two-colour diagram, of the type shown in Fig. 8, was indicative of its chemical type and the thickness of its dust shell. Combined with information on variability and IRAS spectral-type the diagram provides a useful diagnostic. The colour indices for this plot ($`C_{12/25}`$, $`C_{25/60}`$) are defined by:

$$C_{12/25}=2.5\mathrm{log}(\frac{F_{25\mu m}}{F_{12\mu m}}),$$

$$C_{25/60}=2.5\mathrm{log}(\frac{F_{60\mu m}}{F_{25\mu m}}),$$

where the IRAS fluxes are *not colour corrected*.

Following van der Veen & Habing the diagram is divided into regions, each of which is dominated by a certain type of object. Stars with low mass-loss rates will fall in the lower left of the diagram and increasingly higher mass-loss rates are found as one moves to the upper right.

Most of the oxygen-rich stars under discussion occupy region IIIa which van der Veen & Habing found to contain variable stars with moderate oxygen-rich dust-shells. Four stars are located in regions IIIb and IV, indicating thick to very thick oxygen-rich dust-shells and high mass-loss rates. Four oxygen-rich stars are located in region VII, an area dominated by carbon stars, but where van der Veen & Habing also found oxygen-rich stars.

Most of the LRS spectra for oxygen-rich stars in this sample are classified as 2n, indicating that the $`10\mu `$m silicate feature is in emission. The two oxygen-rich variables with the reddest $`[12][25]`$ colour, $`174113154`$ and $`183480526`$, have LRS classifications of 3n, i.e. their spectra show the $`10\mu `$m silicate feature in absorption. In these stars the dust-shell has become so cool and optically thick that the silicate feature has gone into absorption.

Most of the carbon-stars fall in region VII, which is characterised by variable stars with well developed carbon-rich dust-shells. All of these stars have an LRS classification of 4n, indicating that the SiC $`11\mu `$m feature is in emission. The Red Rectangle, $`061761036`$, (discussed above) is located in region IIIa, as is $`23166+1655`$. The latter object has an LRS classification of O2, probably due to self absorption of the SiC feature in the very think shell (Loup et al. 1993).

## 6 Bolometric Magnitudes and Distances

### 6.1 Apparent Bolometric Magnitudes

Given that these sources were selected on the basis of their high fluxes and therefore small distances, together with the fact that the bulk of their energy is emitted at infrared wavelengths it is reasonable to assume that the effects of interstellar extinction are negligible. The apparent bolometric magnitude, $`m_{bol}`$, were calculated by integrating under a spline curve (Hill 1982), fitted to the $`J`$, $`H`$, $`K`$, $`L`$, 12- and 25-$`\mu m`$ fluxes as a function of frequency. The end-points were dealt with by extrapolating a line, which joins the $`J`$ flux and the point lying halfway between the $`H`$ and $`K`$ fluxes, to zero flux. Zero flux at zero frequency was assumed. The derived values of $`m_{bol}`$ are listed in Table 7 for the sources with sufficient photometry.

For the very red sources $`23166+1655`$ and $`183480526`$, which have no measurements in the $`J`$ and/or $`H`$ bands, zero flux at $`J`$ and $`H`$ was assumed when making the spline fit.

To estimate an upper-limit to the amount of flux not accounted for by ignoring the tail in the flux distribution, a Wien approximation, $`A\nu ^3e^{b\nu }`$ (where A and b are constants), was assumed for flux at frequencies higher than $`J`$ with the same slope as the spline at $`J`$. In more than 75 percent of the cases the difference between the fluxes calculated in this way and the simple spline fit was less than 1 percent. $`204400105`$, which has the smallest $`K[12]=3.16`$, shows the maximum difference of 5 percent. Neglecting $`\mathrm{H}_2\mathrm{O}`$ absorption between the $`J`$,$`H`$,$`K`$ and $`L`$ bands (Robertson & Feast 1981) will have a negligible effect on the bolometric magnitude estimated for stars whose flux originates largely from the circumstellar dust-shell.

The bolometric corrections at 12$`\mu m`$, $`BC_{12}`$, is plotted against the $`K[12]`$ colour in Fig. 9, together with curves given by:

$$BC_{12}=15.4e^{0.70(K[12])}+8.10,$$
(4)
and

$$BC_{12}=7.59e^{0.43(K[12])}+8.41,$$
(5)
which are least-square fits of this form to the oxygen- and carbon-rich star data respectively. The rms difference between the data and these curves is 0.04 and 0.06 mag respectively. For the remaining sources $`m_{bol}`$ was calculated from  and $`K[12]`$ using equations 4 and 5.

The apparent luminosities derived here are typically rather smaller than those estimated by JK89, because of the colour corrections applied to the IRAS photometry.

### 6.2 Absolute Bolometric Magnitudes and Luminosities

For the 47 stars with periods the absolute bolometric magnitudes, $`M_{bol}`$, was calculated from the period-luminosity (PL) relation, equation 6, derived for LMC variables (Feast et al. 1989) assuming an LMC distance modulus of 18.60 mag (Whitelock et al. 1997b) and the results given in Table 7.

$$M_{bol}=2.34\mathrm{log}P+1.26$$
(6)
The average luminosity, for all stars with periods less than 1000 days is $`(0.98\pm 0.17)\times 10^4L_{}`$, or $`M_{bol}=5.23\pm 0.19`$. For stars with unknown periods this mean value is assumed; which seems reasonable since the carbon stars seem to have periods less than 1000 days and the remaining oxygen-rich stars have $`[12][25]`$ colours suggesting that they too have periods less than a 1000 days. This value of the mean luminosity suggests that the $`10^4L_{}`$ often assumed for the luminosity of the tip of the AGB is consistent with the application of the PL relation to this type of AGB star. Note this value may not be applicable to the Red Rectangle and the Egg Nebula, since they may have left the AGB already.

It is not really clear at this stage how well the various PL relations fit long period variables or if different relations should be used for the oxygen- and carbon-rich stars. Provided this caveat is born in mind when the results are interpreted it seems reasonable to use the luminosities derived as described above.

Note, in particular, that the PL relation used above was derived for stars with periods less than 420 days and it may not be appropriate to use it for the present sample. There certainly are long period AGB ($`P>420`$ days) stars in the LMC that lie above the PL relation (e.g. Feast et al. 1989). However, recent work on dust-enshrouded stars in the LMC suggest that these AGB stars are over-luminous because they are undergoing hot bottom burning and that more normal long period ($`P400`$ day) carbon- and oxygen-rich AGB stars fall very close to the PL derived for shorter period oxygen-rich stars (Whitelock & Feast 2000).

### 6.3 Distances

Distances to all the stars in the sample (excluding the Red Rectangle and the Egg Nebula) were calculated from $`m_{bol}`$ and $`M_{bol}`$ and are listed in Table 7. These distances are typically larger than those estimated by JK89. As discussed above the apparent luminosities used here are smaller than JK89’s, while the absolute luminosities are on average the same, except for the 3 stars with $`P>1000`$ where they are brighter. Thus we find sources with distances up to 2 kpc.

The other plausible way of estimating distances to OH/IR stars is the phase lag method, although it is notoriously difficult and time consuming. Van Langevelde, van der Heiden & Schooneveld (1990) measured phase lag distances for two of the stars under discussion. For $`01037+1219`$ (WX Psc) they get $`0.74\pm 0.15`$ kpc, identical to the PL distance given in Table 7. For $`183480526`$ (OH26.5+0.6) they give two values, $`1.44\pm 0.27`$ and $`1.30\pm 0.35`$ kpc, both somewhat closer than the 1.9 kpc listed in Table 7. Given the uncertainties in phase lag and the PL method this difference is not problematic. West (1998) has measured a phase lag distance for $`174113154`$ (OH357.3–1.3) of $`1.2\pm 0.4`$ kpc, in agreement with the distance of 0.99 kpc listed in Table 7.

## 7 Mass-loss rates

Mass-loss rates, $`\dot{M}`$, were $`estimated`$ using the expression given by Jura (1987):

$$\dot{M}=1.7\times 10^7v_{15}r_{kpc}^2L_4^{1/2}F_{\nu ,60}\overline{\lambda }_{10},$$
(7)
where $`v_{15}`$ is the outflow velocity of the gas in units of 15 $`\mathrm{km}\mathrm{s}^1`$, $`r_{kpc}`$ the distance to the star in kpc, $`L_4`$ the luminosity in units of $`10^4L_{}`$, $`F_{\nu ,60}`$ the flux measured at $`60\mu `$m in Jansky and $`\overline{\lambda }_{10}`$ the mean wavelength of light emerging from the star and its circumstellar dust-shell in units of $`10\mu `$m.

Equation 7 assumes all the grains have an emissivity of $`150\mathrm{cm}^2\mathrm{g}^1`$ at $`60\mu `$m and a dust to gas ratio of $`4.5\times 10^3`$ (Jura 1986). The mean wavelength of light emerging from the star and its circumstellar dust-shell, $`\overline{\lambda }_{10}`$, is given by:

$$\overline{\lambda }_{10}=\frac{_0^{\mathrm{}}\lambda F_\lambda 𝑑\lambda }{_0^{\mathrm{}}F_\lambda 𝑑\lambda }$$
(8)
where $`F_\lambda `$ the monochromatic flux density in units of energy received per unit area per unit time per unit wavelength interval.

For the stars for which the function $`F_\nu `$ could be estimated (section 6.1), $`\overline{\lambda }_{10}`$ was calculated from equation 8 and is listed in Table 7. The relationship between $`\overline{\lambda }_{10}`$ and $`K[12]`$ is illustrated in Fig. 10 and can be expressed as:

$$\overline{\lambda }_{10}=c+b(K[12])+a(K[12])^2,$$
(9)
where a, b and c are constants.

Least square fits of equation 9 to the data results in the following values (a,b,c) for oxygen-rich (–0.0110,0.284,–0.392) and carbon-rich (–0.0137,0.338,–0.513) stars, respectively, with an rms error of $`0.05\mu `$m for both chemical types. The discrepant carbon star is the Red Rectangle which was omitted from the solution.

Equation 9 with the appropriate parameters is used to estimate $`\overline{\lambda }_{10}`$ for the remaining stars in the sample; the results are listed in Table 7

The outflow velocities, $`v_{15}`$, were taken to be the expansion velocities of the shell obtained from CO($`10`$) line measurements. Appropriate values from the literature are listed in Table 7 with references. For the three stars without measurements an average outflow velocity, for stars with $`v_{out}<25\mathrm{km}\mathrm{s}^1`$, of $`17\pm 4\mathrm{km}\mathrm{s}^1`$ was used. Jones et al. (1983, henceforth JHG) have shown there to be a correlation between outflow velocity and luminosity for a sample of OH/IR stars (see fig. 1 in JHG). The outflow velocity versus luminosity for stars with periods in our sample is illustrated in Fig. 11. It shows no obvious correlation, which should not be entirely unexpected, considering the large spread in outflow velocity at a given luminosity (up to $`10kms^1`$) in fig. 1 of JHG and the small range in $`M_{bol}`$ in Fig. 11 (compared to that in fig. 1 of JHG).

The mass-loss rates for 58 of the stars in the sample were calculated from equation 7 and are listed in Table 7 together with mass-loss rates, $`\dot{M}_{model}`$, taken from the compilation by Loup et al. (1993). Where $`\dot{M}_{model}`$ was determined from detailed models of the CO emission. A comparison of these values is illustrated in Fig. 11 which shows that they agree within a factor of 3. This is reasonable given the uncertainty in the various assumptions used to derive equation 7. $`174113154`$ does not appear in Fig. 11, it has an extremely high mass-loss rate, $`\dot{M}>10^4M_{}yr^1`$, from equation 7, much higher than the value, $`0.87\times 10^5M_{}yr^1`$, given in Loup et al. (1993). Kastner (1992) applies a revised model of circumstellar CO emission, and suggests that it is losing mass at a rate greater than $`\dot{M}>10^5M_{}yr^1`$, consistent with the value derived here.

Fig. 12 shows the relationship between mass-loss rate and $`K[12]`$ ($`\dot{M}`$ correlates with other colours, but none as strongly as this). This relation for Miras has been discussed by Whitelock et al. (1994) and Le Bertre & Winters (1998). A least squares fit can be made to the data in Fig. 12 assuming the form:

$$\mathrm{log}(\dot{M})=\frac{a}{(K[12]+b)}+c,$$
(10)
where a, b and c are constants.

This results in values (a,b,c) for carbon- (–2.89, –1.79, –4.32) and oxygen-rich (–21.34, 3.00, –2.40) stars respectively with rms uncertainties of 0.14 and 0.16. This expression can then be used to estimate mass-loss rate for the RAFGL 1406, which appear in brackets in Table 7.

Fig. 13 shows that there is a correlation between mass-loss rate and the period of pulsation. Some dependence of mass-loss rate on period might be introduced by using the PL relation to calculate the quantities $`L_4`$ and $`r_{kpc}`$ for equation 7. If we assume, for the sake of argument, that these sources are all the same luminosity, then use of the PL will result in a mass-loss/period dependence of the form $`\dot{M}P^{0.47}`$, i.e. a much weaker effect than is actually seen. A correlation between period and mass-loss rate was also found by Whitelock et al. (1994).

The average mass-loss rate for stars in this sample is $`\dot{M}1.5\times 10^5M_{}yr^1`$. Stars with $`P>1000`$ day (i.e. $`[12][25]>1.7`$) have $`\dot{M}1.3\times 10^4M_{}yr^1`$ and stars with $`P<1000`$ day (i.e. $`[12][25]<1.7`$) have $`\dot{M}0.87\times 10^5M_{}yr^1`$. It will be seen later that stars with $`P>1000`$ day probably have main-sequence masses greater that $`4M_{}`$, while those with $`P<1000`$ day have main-sequence masses between 1 and $`2M_{}`$. Thus intermediate mass stars lose mass an order of magnitude faster than do low mass-stars. The average mass-loss rates of carbon- and oxygen-rich (including S-type) stars are $`\dot{M}(1.0\pm 1.0)\times 10^5M_{}yr^1`$ and $`\dot{M}(2.0\pm 1.2)\times 10^5M_{}yr^1`$ respectively. This slight difference may be a consequence of assuming the same grain emissivity and gas to dust ratio, in equation 7, for the two chemical types. Within the limits of the assumptions used here, the mass-loss rates of the carbon and oxygen-rich stars are plausibly identical.

One of the stars in this sample, RAFGL 1406 (IY Hya), is a known super-lithium-rich carbon star (i.e. it has $`\mathrm{log}\frac{n(Li)}{n(H)}>8`$) and it is estimated that about 2 percent of carbon-stars are of this type (Abia et al. 1993a). Furthermore there is evidence that these stars are currently important contributors to lithium in the Galaxy (Abia et al. 1995 and 1998).

Abia et al. (1993b) estimate a lithium production rate of $`<\dot{m}_7>=4.5\times 10^{14}M_{}yr^1`$ per carbon star with a return rate into the interstellar medium $`P_7=1.9\times 10^9M_{}pc^2Gyr^1`$. The sample of stars they studied have a typical mass-loss rate of $`2\times 10^7M_{}yr^1`$, i.e. more than a order of magnitude less than the average mass-loss rate of carbon-rich DE-AGB stars. Note that Abia et al. (1993b) derive a average progenitor mass for the carbon stars in their sample of 1.2-1.6$`M_{}`$, which is similar to the average progenitor mass estimated for this sample (see later). Assuming a super-lithium-rich carbon star to have a typical lithium abundance of $`\mathrm{log}\frac{n(Li)}{n(H)}=7`$ (see fig. 3 in Abia et al. 1993a) we estimate a lithium production rate from carbon-rich DE-AGB stars of $`1.3\times 10^{13}M_{}yr^1`$ per carbon star, with the fraction of super-lithium-rich carbon stars in this sample taken to be 0.02. The total return rate of lithium by these stars into the interstellar medium is then about $`1.0\times 10^9M_{}pc^2Gyr^1`$. This indicates that carbon-rich DE-AGB stars may contribute up to a third of the lithium in the Galaxy produced from carbon-rich AGB stars. However, granting the underlying assumptions and uncertainty in the numbers above, further investigation into the problem is necessary (see e.g. Abia et al. 1993b).

## 8 Scale-height, Progenitor Mass and Lifetime

### 8.1 Scale-height

Assuming the usually exponential distribution of stellar densities the number of stars per height interval, $`\rho `$, varies with distance from the galactic plane, $`z`$, as:

$$\rho =\rho _0e^{|z|/H},$$
(11)
where $`\rho _0`$ is the density of stars at $`z=0`$ and $`H`$ is the scale-height. For such a population the average height, $`\overline{|z|}`$, is equal to the scale-height, $`H`$, provided the sample is reasonably complete (as JK89 suggest it is).

The sample of stars discussed here were necessarily selected from a limited region of the sky by JK98. Taking this into account we use only stars with $`14^{}<l<230^{}`$ and find that $`H=235`$ pc.

Alternatively, the fraction of stars, $`f_{|z|}`$, within a given $`|z|`$ of the galactic plane is given by:

$$f_{|z|}=1e^{|z|/H}.$$
(12)
$`f_{|z|}`$ is plotted in Fig. 14 for each 100 pc interval for the sample of stars with $`14^{}<l<230^{}`$. A curve of the form given by equation 12, fitted by least squares, is also shown. This fits the data very well (rms difference 0.03) and provides a scale-height of $`H=236\pm 10`$ pc, comparable to the estimate made above. This value will be used below as representative of the scale height of the sample.

### 8.2 Progenitor Mass

The scale height derived above can be used to estimate a typical progenitor mass for the AGB population. A comparison with Miller & Scalo (1979 - their table 1) indicates an average main-sequence mass $`M_i1.2M_{}`$ for the group under discussion, although it is difficult to estimate the uncertainty on this value.

An alternative approach uses the data in table 3 of Miller & Scalo (1979) which gives the number of main-sequence stars pc$`{}_{}{}^{2}(\mathrm{log}(M))_{}^{1}`$, $`\varphi _i`$, as a function of $`\mathrm{log}(M)`$. A maximum likelihood straight line fit to these data over the range $`0.14<\mathrm{log}(M)<0.92`$ has a slope $`\alpha =4.20\pm 0.32`$, a zero-point $`\beta =1.20\pm 0.14`$ and fits the data fairly well. Let $`M_{lower}`$ and $`M_{upper}`$ be the lower and upper progenitor mass limits respectively of the stars in this sample. Then it is straight forward to show that the average mass, $`\overline{M}`$, in $`M_{}`$ is given by:

$$\overline{M}=\frac{\alpha }{\alpha +1}\frac{M_{upper}^{\alpha +1}M_{lower}^{\alpha +1}}{M_{upper}^\alpha M_{lower}^\alpha }$$
(13)
To get an idea of the values of $`M_{lower}`$ and $`M_{upper}`$ the relation between $`M_{lower}`$, $`M_{upper}`$ and the $`average`$ scale-height $`\overline{H}`$ was investigated. $`\overline{H}`$, for a given mass range $`[M_{lower},M_{upper}]`$ is simply given by:

$$\overline{H}=\frac{_{M_{lower}}^{M_{upper}}H(M)𝑑N(M)}{_{M_{lower}}^{M_{upper}}𝑑N(M)}$$
(14)
where $`H(M)`$ is the scale-height as a function of mass $`M`$ (table 1 from Miller & Scalo), $`dN(M)=\frac{\varphi _i}{Mln(10)}dM`$ the number of stars$`\mathrm{pc}^2`$ in the infinitesimally small interval $`dM`$ centred around $`M`$.

Assuming an average scale-height, $`\overline{H}_0`$, then for a given upper mass limit $`M_{upper}`$ equation 14 can be solved numerically for $`M_{lower}`$, since the functions $`H(M)`$ and $`\varphi _i`$ are known. This was done for $`\overline{H}_0=H=236`$ pc. The resulting relation between $`M_{lower}`$ and $`M_{upper}`$ is plotted in Fig. 15, from which it is clear that for $`M_{upper}>3M_{}`$, $`M_{lower}`$ has a constant value of $`0.96M_{}`$.

Now since we expect to find dust envelope (DE) AGB main-sequence progenitors up to about $`8M_{}`$ and since the value of $`\overline{M}`$, from equation 13, is relatively insensitive to the value of $`M_{upper}`$ if $`\frac{M_{upper}}{M_{lower}}>3`$, values of $`M_{lower}=0.93M_{}`$ and $`M_{upper}=8M_{}`$ will be assumed. It is important to note that the value of $`\overline{M}`$ is essentially proportional to the value of $`M_{lower}`$ and almost independent of $`M_{upper}`$ if $`\frac{M_{upper}}{M_{lower}}>3`$. Also when calculating the surface density of progenitors, it will be seen that this is very sensitive to the value of $`M_{lower}`$. Changing the assumed average scale-height to 216 pc or 256 pc (i.e. effectively changing the assumed scale-height by about $`2\sigma `$) gives values of $`M_{lower}`$, assuming $`M_{upper}=8M_{}`$, of $`0.99M_{}`$ and $`0.86M_{}`$ respectively. Thus an error in $`M_{lower}`$ of $`0.07M_{}`$ will be assumed. This then gives an average mass from equation 13 of $`(1.22\pm 0.18)M_{}`$.

From the discussion above, an average scale-height of 236 pc implies a progenitor mass range starting from about $`0.96M_{}`$ up to possibly $`8M_{}`$. Hence using equation 18 (below), one can estimate the ratio of the number of stars in the mass range $`2.0M_{}<M_i<8.0M_{}`$ to the number of stars in the mass range $`1.0M_{}<M_i<2.0M_{}`$. This ratio is only 0.06. The three stars with periods greater than 1000 days have an average height above the galactic plane of $`20\pm 2`$ pc. This indicates that the progenitors of these stars were probably more massive stars, with $`M_i>4.0M_{}`$, as can be seen from Miller & Scalo’s table 1. Interestingly the ratio of stars with $`P>1000`$ day to the rest of the sample (which are expected to have $`P<1000`$ day from their colours) is 0.05, comparable to the ratio calculated above. The ratio, from equation 18, of the number of stars in the mass range $`0.93M_{}<M_i<1.0M_{}`$ to the number stars in the mass range $`1.0M_{}<M_i<2.0M_{}`$ is about 0.4. Thus the majority of the stars in this sample have main-sequence masses in the range $`1.0M_{}<M_i<2.0M_{}`$.

Limited discussion in the literature also supports this conclusion. Thus Kahane et al. (2000) estimate the main-sequence progenitor of $`09452+1330`$, one of the stars in our sample, to be of low mass ($`M_i2M_{}`$) by comparing the predictions from AGB stellar models of solar composition and various initial mass with the observed chlorine isotopic ratio in its circumstellar shell. Also Nishida et al. (2000) estimate the progenitor mass of three carbon-rich Miras, with $`P>450`$ day, in Magellanic Cloud clusters to be $`1.5`$-$`1.6M_{}`$ by main-sequence isochrone fitting.

The velocity dispersion of the stars with $`P<1000`$ provides another indication of initial mass. Fifty three sources have velocities, relative to the local standard of rest, tabulated by Loup et al. (1993), they show a dispersion $`\sigma =24\pm 2\mathrm{km}\mathrm{s}^1`$. After correcting for galactic rotation, using the constants from Feast & Whitelock (1997), this reduces to $`\sigma =19\pm 2\mathrm{km}\mathrm{s}^1`$. Table 8 lists the velocity dispersions (calculated from table 2 of Delhaye 1965) and initial masses (from Allen 1973) for main-sequence stars between A0 and G5. The velocity dispersion thus indicates that the mean sample has $`M_i<2M_{}`$, consistent with numbers given above.

### 8.3 Surface Density

Fig. 16 shows the surface distribution of the stars in the galactic plane, calculated from the distances derived previously. As can be seen, most of the stars (90 percent) are within $`1.4`$ kpc of the sun. There are very few in quadrant $`IV`$ due to positional selection effects in the original sample. The surface density, $`\sigma _i`$ in region $`i`$, $`i=I,II,III`$, was estimated from:

$$\sigma _i=\frac{N_i}{\frac{1}{4}\pi (1.0)^2}\mathrm{stars}\mathrm{kpc}^2$$
(15)
where $`N_i`$ is the number of stars in a particular region $`i`$.

The average surface density, $`\sigma _{DEAGB}`$ of DE-AGB stars in the solar neighbourhood was then estimated from:

$$\sigma _{DEAGB}=\frac{2\sigma _I+\sigma _{II}+\sigma _{III}}{4}$$
(16)
where $`\sigma _I\sigma _{IV}`$, assuming symmetry about $`l=0`$.
Using equations 15 and 16 we find $`\sigma _{DEAGB}=15\pm 4`$ $`\mathrm{stars}\mathrm{kpc}^2`$.

### 8.4 Lifetime

Assuming the average lifetime of stars in any given phase of evolution is proportional to the number of stars in that phase, it follows that:

$$\overline{\tau }_{DEAGB}=\frac{\sigma _{DEAGB}}{\sigma _{MS}}\overline{\tau }_{MS}$$
(17)
where $`\overline{\tau }_{DEAGB}`$ is the average lifetime of DE-AGB stars, $`\overline{\tau }_{MS}`$ and $`\sigma _{MS}`$ the average lifetime and surface density, respectively, of the progenitors of DE-AGB stars on the main-sequence.

From the linear relation between $`\mathrm{log}(\varphi _i)`$ and $`\mathrm{log}(M)`$ it is straight forward to show that $`\sigma _{MS}`$ in the mass range $`[M_{lower},M_{upper}]`$, is simply given by:

$$\sigma _{MS}=\frac{10^\beta }{\alpha ln(10)}[M_{upper}^\alpha M_{lower}^\alpha ]$$
(18)
where $`\sigma _{MS}`$ is in units of $`\mathrm{stars}\mathrm{pc}^2`$.
It is clear from equation 18 and the value of $`\alpha `$ that $`\sigma _{MS}`$ is very sensitive to the value of $`M_{lower}`$ and $`M_{}^{\alpha }{}_{lower}{}^{}`$ if $`\frac{M_{upper}}{M_{lower}}>3.0`$. Using the upper and lower mass limits assumed above, gives $`\sigma _{MS}=2.22\pm 0.76`$ $`\mathrm{stars}\mathrm{pc}^2`$.

Miller & Scalo’s table 1 also gives $`\mathrm{log}(\tau _{MS})`$ as a function of $`\mathrm{log}(M)`$. $`\overline{\tau }_{MS}`$ is simply given by:

$$\overline{\tau }_{MS}=\frac{_{M_{lower}}^{M_{upper}}\tau _{MS}(M)𝑑N(M)}{_{M_{lower}}^{M_{upper}}𝑑N(M)}$$
(19)
Equation 19 was integrated numerically with the above mass limits, giving $`\overline{\tau }_{MS}=(5.5\pm 1.5)\times 10^9`$ years. The errors above were estimated by changing the average scale-height (and hence $`M_{lower}`$) to 216 and 256 respectively.

Using equation 17 an average lifetime for the stars in this sample was estimated, giving $`\overline{\tau }_{DEAGB}=(3.7\pm 1.9)\times 10^4`$ years. This value is in excellent agreement with lifetimes of similar stars at the tip of the AGB estimated by various other investigators (see Habing 1996 p. 183).

## 9 Evolutionary Status and Conclusions

The stars under discussion were selected to isolate high mass-loss, thick dust-shell AGB stars within the Solar Neighbourhood. This study indicates that the selection criteria employed were fairly successful in this respect. The $`JHKL`$ photometry, the IRAS two-colour diagrams and other published data prove that the majority of these sources are cool variable stars with moderate to high mass-loss rates.

As shown above most of these stars will have had a main-sequence progenitor with $`1M2`$ $`M_{}`$. While the average progenitor mass is $`1.3M_{}`$. Using this average value and an empirical “initial/final mass”-relation (Weidemann 1987) an average final white dwarf mass can be estimated at $`0.56M_{}`$.

Assuming the core masses of these stars are close to the final white dwarf mass, it is found that the thermal interflash period ($`t_{if}`$) will be, $`t_{if}10^5`$ years, using the core-mass interflash period relation given by Boothroyd & Sackmann (1988) (their equation 6 which is valid for solar composition). A more recent core-mass interflash period relation given by equation 11a of Wagenhuber & Groenewegen (1998) also suggests a thermal interflash period of $`t_{if}10^5`$ years; this is provided there is no hot bottom burning and that the star has already undergone several helium shell flashes (see Wagenhuber & Groenewegen). The first assumptions is reasonable as stars under 2 $`M_{}`$ do not experience hot bottom burning. The second is also an essential assumption for using the final white dwarf mass for the current core mass and obviously the carbon stars, at least, have undergone several flashes.

This interflash period, $`t_{if}10^5`$, is more than a factor of two larger than the average lifetime of these DE-AGB stars, $`\overline{\tau }_{DEAGB}=(3.7\pm 1.9)\times 10^4`$ years, deduced above. Hence we can conclude that these DE-AGB stars will experience at most one thermal pulse before leaving the TP-AGB phase.

The average mass-loss rate, $`\dot{M}_{ave}`$ for stars in the sample was found to be $`1.5\times 10^5M_{}yr^1`$. Thus the total mass lost, $`M_{lost}`$, by an average star in this phase is $`M_{lost}=\dot{M}_{ave}\overline{\tau }_{DEAGB}0.5M_{}`$. Which is very close the average initial-final mass difference for these stars. However, this value for $`M_{lost}`$ is very uncertain due to the uncertainty in the mass-loss rates. Nevertheless, this indicates that a significant fraction of the envelope mass is lost during this high mass-loss dust-enshrouded TP-AGB phase.

The total rate of mass returned into the interstellar medium, $`\dot{M}_{ISM}`$, by stars in this phase is $`\dot{M}_{ISM}=\dot{M}_{ave}\sigma _{DEAGB}2.3\times 10^4M_{}kpc^2yr^1`$. Jura & Kleinmann (1990) estimate that M supergiants in the solar neighbourhood return mass to the interstellar medium at only a rate between $`1\times 10^5`$ and $`3\times 10^5M_{}kpc^2yr^1`$. Thus AGB stars must be important contributors of material to the interstellar medium in the solar neighbourhood.

### Acknowledgements

We would like to thank Chris Koen, Michael Feast and Peter Wood for helpful discussions and advice, Luis Balona for his Fourier analysis program and Tom Lloyd Evans for allowing us to use observations in advance of publications. We are also grateful to, Robin Catchpole and Dave Laney for making near-infrared observations on our behalf. This paper made use of the SIMBAD database and observations made from the South African Astronomical Observatory.

## Appendix
Average Magnitudes and Colours

### Average Magnitudes

For variable stars, such as those discussed here, it is often useful to know the average magnitude at a particular wavelength. The average magnitude as discussed in section 3, was taken as the mean of the minimum and maximum magnitude measured for the particular star. This average, $`A_M`$, is compared in Table A1 with the average magnitude, $`\overline{M}`$, derived from the fit given by equation 4.1 for the 19 stars with sufficient data. The average absolute difference, $`\mathrm{\Delta }_M`$, between these two averages, for each near-infrared band, with the sample standard deviation, $`\sigma _\mathrm{\Delta }`$, on this average is also listed in Table A1. From these values for $`\mathrm{\Delta }_M`$ it can be seen that $`A_M`$ represents the average magnitude fairly well. Obviously in instances where there are insufficient observations to cover a complete cycle, $`A_M`$ will be a poor representation of the mean.

### Average Colours

The stars in this sample not only vary in light output, but also in colour. In section 5 the average colour of a star is taken as the mean of the value calculated at each Julian Date, $`A`$. This average is then compared to the average colour, $`A_f`$, derived from the Fourier fits given by equation 4.1. $`A_f`$ is simply given by:

$$A_f=<M_1M_2>=<M_1><M_2>=M_{1,av}M_{2,av}$$
(20)
where $`M_1`$ and $`M_2`$ are the magnitudes in the bands making up the colour. $`M_{1,av}`$ and $`M_{2,av}`$ are the average magnitudes in these bands obtained from the Fourier fits given by equation 4.1.

Both the averages, for each near-infrared colour, are listed in Table A2 for the 18 stars with sufficient data. Also listed in the Table are the average absolute difference, $`\mathrm{\Delta }_A`$ and sample standard deviation, $`\sigma _\mathrm{\Delta }`$, on this average. The two values are quite comparable.