text:
# Hydrogen model atmospheres for white dwarf stars

## 1 Introduction

White dwarf (WD) stars represent the most common end-product of the stellar evolution. On average, WDs have a mass of $`0.6M_{}`$, a central density of $`10^6`$ gr cm-3 and a core composed of carbon and oxygen resulting from helium burning during the preceding phases of the evolution (D’Antona & Mazzitelli 1990). Due to the high degree of electron degeneracy in the interior of WDs, the outer nondegenerate layers, and specially their atmospheres, control the loss of energy into outer space and therefore, the cooling of these objects. From a technical view point, detailed atmosphere models provide an adequate surface boundary condition for the computation of evolutionary stellar models. Besides, the atmosphere models are necessary to determine the emergent flux from the star, which is required to interpret the luminosities and colours of observed WDs.

Current interest in WDs is motivated by two reasons. These stars are candidates for the Galactic dark matter and can also be used to estimate the age of the star system they belong to. The recent detection of gravitational microlensing events in the direction of the Magellanic Glouds (Alcock et al. 1997, Alcock et al. 1999, Renault et al. 1997) suggests that very low mass stars, WDs among others, might constitute a significant fraction of the mass of the Galatic halo. The luminosity function necessary to study this interpretation requires the calculation of cooling tracks using accurate atmosphere models.

The spectral evolution of a WD during its cooling can be used to calculate the age of the star. The first detailed calculation of WD evolution was performed by Iben & Tutukov (1984), for the cooling of a carbon-oxygen WD from the planetary nucleus stage through the stage of complete internal crystallization. Recent works (Hansen 1999; Chabrier 1999) show that the time scales involved, particularly in an advanced stage of the cooling process, depend strongly on the atmosphere composition (among other factors such as the stellar mass), which through the opacity regulates the energy loss from the stellar interior. Thus, for instance, the strong opacity from hydrogen molecules causes a slow cooling of the WDs with hydrogen-rich atmosphere, resulting these objects potentially observable for times comparable to the Hubble time (Richer et al. 2000). The technique based on the observation of a low-luminosity cutoff has been applied to estimate the age of the local Galatic disk (Winget et al. 1987; Wood 1992) and the age of clusters (von Hippel, Gilmore & Jones 1995; von Hippel & Gilmore 2000).

The majority of WDs have pure hydrogen atmospheres due to gravitational settling, which is efficient at removing helium and heavier elements from atmosphere toward inner layers (Fontaine & Michaud 1979; Muchmore 1984; Iben & MacDonald 1985; Althaus & Benvenuto 2000). Convective mixing could lead to a helium-rich atmosphere in cool WD when the mass of the superficial hydrogen layer is small enough (Koester 1976; Vauclair & Reisse 1977; D’Antona & Mazzitelli 1979). In particular, Bergeron, Ruiz & Leggett (1997) have reported evidence for convective mixing in cool WD atmospheres. However, hydrogen-rich atmospheres are observationally most abundant than those of helium or other with metals, in part due to the long cooling times. Furthermore, these atmospheres present a structure which is insensitive to small admixtures of helium (Hansen 1999).

In this context, recent model atmosphere calculations (Hansen 1998, 1999; Saumon & Jacobson 1999) suggest that WDs of very low luminosity, and therefore very old, become bluer and not redder as they age (and their effective temperature decreases), resulting their detection easier. This being the case, an important number of these objects should be observed with the new generation of telescopes.

Another important question of current interest is the evolution of WDs with very low mass having helium cores (He WD). Such stars, due to the low stellar mass, would be formed from some binary systems. Indeed, numerous binary configurations of this type have been observed lately (Marsh 1995; Moran, Marsh & Bragaglia 1997; Edmonds et al. 1999). In particular, low mass WDs have been detected in binary systems composed of a WD star and a pulsar. A particularly interesting binary system with these characteristics is the millisecond pulsar J1012+5307 and its WD companion, because it is possible to determine the properties of the pulsar from the analysis of the WD (Van Kerwijk, Bergeron & Kulkarni 1996). The study of these binaries would allow us to determine an upper limit to the neutron-star mass (and in this way, to constrain the still uncertain properties of the equation of state of the cool nuclear matter) and on its age (with implications for the evolution of its magnetic field). Obviously, to obtain accurate results from the study of these systems it is necessary to have suitable evolution models and, in particular, detailed model atmospheres for WDs.

The primary aim of this paper is to present results of model calculations from an independent atmosphere code. This will be applied to determine accurate cooling tracks of WDs, by providing appropriate outer boundary conditions to evolution models. Our model atmospheres with pure hydrogen compositions cover extensive ranges of effective temperature ($`2500T_{eff}60000`$ K) and surface gravity ($`6.5\mathrm{log}g9.0`$). Although several sets of models available in the literature (cf. Bergeron, Ruiz & Leggett 1997; Saumon & Jacobson 1999; Hansen 1999; Chabrier 1999) have already investigated different parts of this ($`T_{eff},g)`$ region, we offer here an overall analysis. We have considered it important to extend the study to the very low gravities (log $`g=6.5`$) which are found in early stages of the WD formation, and also in the cool stages of He WDs with very low mass, which have received little or no attention so far. In particular, we have also computed model atmospheres following the cooling tracks of He WDs with hydrogen envelopes and several masses.

The results presented here are based on ideal equations of state. The incorporation of density effects in our code using the occupation probability formalism of Hummer & Mihalas (1988) is in progress. However, nonideal effects can be neglected for hydrogen WD models with T$`{}_{eff}{}^{}4000`$ K (Bergeron, Saumon & Wesemael 1995a). They increase for lower Teff and are significant for T$`{}_{eff}{}^{}2500`$ K (Saumon & Jacobson 1999).

In Section 2 we derive explicit analytical expressions for the thermodynamic quantities required in the calculations. The method used for the model atmosphere computations is presented in Section 3, together with the opacities and convective theory adopted. In Section 4 we display the atmosphere structure of WDs obtained from our code. Comparisons with observations and other recent model results, through the analysis of colour diagrams, are also presented there. Our conclusions are given in Section 5.

## 2  Thermodynamic quantities

The computation of convective flux of energy requires some properties of the gas such as the adiabatic gradient and the specific heat. To determine them, we consider a mixture of H2, H, H+, e-, and radiation occuping a volume $`V`$, in thermodynamic equilibrium at temperature $`T`$. In the present calculations, species with low abundances can be ignored (Saumon & Chabrier 1992).

Let $`N`$ be the number of hydrogen nuclei or protons. Then, the number conservation and charge neutrality constrains demand that the particle density is

$`{\displaystyle \frac{xN}{V}},`$ (1)
with

$`x=1+x_1x_2,`$ (2)
where $`x_1=n_{H^+}/n`$ and $`x_2=n_{H_2}/n`$ (with $`n=N/V`$) are the number concentration of H+ and H2. We use the ideal state equation to express the gas pressure. The total pressure, including the radiation contribution, is then

$`P={\displaystyle \frac{xNkT}{V}}+{\displaystyle \frac{1}{3}}aT^4.`$ (3)
We adopt a simple approach for the total energy in the system

$`U={\displaystyle \frac{3}{2}}xNkT+\left[x_1\chi _1+\left({\displaystyle \frac{1}{2}}x_2\right)\chi _2\right]N`$ (4)
$`+\left[kT+kT_v\left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{e^\theta 1}}\right)\right]x_2N+aT^4V,`$ (5)
where $`\theta =T_v/T`$. The first term on the right-hand side is the translational energy, which assumes the particles to behave as classical point particles. The second term gives the potential energy of the hydrogen nuclei, which increases with the dissociation of H2 molecules and the ionization of H atoms according to the corresponding binding energies $`\chi _1`$ and $`\chi _2`$. For simplicity, we ignore in this term the excited states of these species. The zero point of the energy scale is taken to be the ground state of H2. The third term represents the rotation-vibration energy of H2, where we use a vibration temperature $`T_v=6100`$ K. Finally, the last term takes into account the energy in the radiation field.

Expressions for the adiabatic gradient $`_A=\mathrm{ln}T/\mathrm{ln}P|_S`$, and heat capacity at constant pressure $`C_P=U/T|_P+PV/T|_P`$, can be obtain from differentiation of Equations (3) and (4). We have followed a similar procedure to that of Krishna-Swamy (1961; see also Mihalas 1965) and, therefore, we only show here the results obtained including molecular species in the gas mixture. The expresion for the adiabatic gradient is

$`_A={\displaystyle \frac{a_1}{a_2}},`$ (6)
with

$`a_1=1+4\gamma b_1{\displaystyle \frac{x_1}{\mathrm{ln}P_g}}|_T+b_2{\displaystyle \frac{x_2}{\mathrm{ln}P_g}}|_T`$ (7)
and

$`a_2={\displaystyle \frac{5}{2}}+16\gamma +\left[1+{\displaystyle \frac{\theta ^2e^\theta }{(1e^\theta )^2}}\right]{\displaystyle \frac{x_2}{x}}+b_1{\displaystyle \frac{x_1}{\mathrm{ln}T}}|_{P_g}`$ (8)
$`b_2{\displaystyle \frac{x_2}{\mathrm{ln}T}}|_{P_g},`$ (9)
where

$`b_1=\left({\displaystyle \frac{5}{2}}+{\displaystyle \frac{\chi _1}{kT}}+4\gamma \right){\displaystyle \frac{1}{x}},`$ (10)
$`b_2=\left\{{\displaystyle \frac{5}{2}}+{\displaystyle \frac{\chi _2}{kT}}+4\gamma \left[1+\theta \left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{e^\theta }{1e^\theta }}\right)\right]\right\}{\displaystyle \frac{1}{x}}.`$ (11)
For the heat capacity we obtain

$`{\displaystyle \frac{C_P}{kN}}=\left[{\displaystyle \frac{3}{2}}+12\gamma +\left(1+4\gamma \right)^2\right]x`$ (12)
$`+\left[1+{\displaystyle \frac{\theta ^2e^\theta }{(1e^\theta )^2}}\right]x_2+\left({\displaystyle \frac{3}{2}}+{\displaystyle \frac{\chi _1}{kT}}\right){\displaystyle \frac{x_1}{\mathrm{ln}T}}|_P`$ (13)
$`\left\{{\displaystyle \frac{3}{2}}+{\displaystyle \frac{\chi _2}{kT}}\left[1+\theta \left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{e^\theta }{1e^\theta }}\right)\right]\right\}{\displaystyle \frac{x_2}{\mathrm{ln}T}}|_P.`$ (14)
The partial derivatives of the concentrations $`x_1`$ and $`x_2`$, have been calculated considering that, except in the regime of pressure dissociation, H+ and H2 do not have simultaneously high abundances. Then, using Saha equations for the first and third reactions in Table 1, result

$`{\displaystyle \frac{x_1}{\mathrm{ln}P_g}}|_T={\displaystyle \frac{x_1}{2}}\left(1x_1^{\mathrm{\hspace{0.33em}\hspace{0.33em}2}}\right),`$ (15)
$`{\displaystyle \frac{x_1}{\mathrm{ln}T}}|_{P_g}={\displaystyle \frac{x_1}{2}}\left(1x_1^{\mathrm{\hspace{0.33em}\hspace{0.33em}2}}\right)\left({\displaystyle \frac{5}{2}}+{\displaystyle \frac{\chi _1}{kT}}\right),`$ (16)
$`{\displaystyle \frac{x_1}{\mathrm{ln}T}}|_P={\displaystyle \frac{x_1}{2}}\left(1x_1^{\mathrm{\hspace{0.33em}\hspace{0.33em}2}}\right)\left({\displaystyle \frac{5}{2}}+{\displaystyle \frac{\chi _1}{kT}}+4\gamma \right),`$ (17)
$`{\displaystyle \frac{x_2}{\mathrm{ln}P_g}}|_T=\left({\displaystyle \frac{1}{x_2}}+{\displaystyle \frac{4}{12x_2}}{\displaystyle \frac{1}{1x_2}}\right)^1,`$ (18)
$`{\displaystyle \frac{x_2}{\mathrm{ln}T}}|_{P_g}=\left({\displaystyle \frac{5}{2}}+{\displaystyle \frac{\chi _2}{kT}}\right){\displaystyle \frac{x_2}{\mathrm{ln}P_g}}|_T,`$ (19)
$`{\displaystyle \frac{x_2}{\mathrm{ln}T}}|_P=\left({\displaystyle \frac{5}{2}}+{\displaystyle \frac{\chi _2}{kT}}+4\gamma \right){\displaystyle \frac{x_2}{\mathrm{ln}P_g}}|_T,`$ (20)
where the dependence of the partition functions with the temperature and pressure has been neglected.

Since the adiabatic temperature gradient forms the basis of the Schwarzschild criterion for convective instability, this gradient plays a central role in the computation of atmosphere models. Therefore, we check our expression (6-8) by comparing with results from the equation of state (EOS) for hydrogen developed by Saumon & Chabrier (summarized in Saumon, Chabrier & Van Horn 1995, SCVH). This EOS, considered the best available at present for hydrogen, is based on the free-energy minimization technique (for details, see Graboske, Harwood & Rogers 1969; Fontaine, Graboske & Van Horn 1977; Hummer & Mihalas 1988) and includes a sophisticated treatment of nonideal effects.

Comparisons between values for the adiabatic gradient obtained from SCVH and our EOS are shown in Fig. 1. The curves correspond to models of WD atmosphere with $`T_{eff}=2500`$, 4000, 6000 and 10000 K, and log $`g=8`$. Both present and SCVH calculations agree quite well except for high pressures of the model $`T_{eff}=2500`$ K, where non-ideal effects can be present. For the hottest two models, discrepancies at low pressure occur beyond the regime where SCVH EOS can be applied ($`4<\mathrm{log}P`$(dyn cm-2)$`<19`$). In particular, the upper layers of the $`T_{eff}=6000`$ K atmosphere are dominated by atomic hydrogen (see Section 4.1), and a tendency to a value 0.4 for the temperature gradient (corresponding to an ideal monoatomic gas) is reproduced by our calculations. A partial molecular formation (up to $`n_{H_2}/n=0.12`$) takes place as one proceeds inward in this atmosphere, yielding a dramatic drop in the value of $`_A`$.

Other familiar features can be identified in this figure. For instance, $`T_{eff}=2500`$ K model represents an atmosphere formed mainly by H2. At low temperatures present in its superficial layers, the H2 vibrational modes remain frozen and then, the adiabatic gradient takes a value close to $`2/7=0.286`$, i.e. that due to a fluid of rigid rotating molecules with two degrees of rotational freedom and three degrees of translational freedom. The adiabatic gradient decreases toward deep layers where the molecular gas is partially dissociated.

The low values reached by $`_A`$ in the $`T_{eff}=4000`$ K model are due to a marked sensitivity of the molecular abundance to changes in temperature and pressure. While the behaviour of $`_A`$ in the hottest model ($`T_{eff}=10000`$ K) is produced by different ionization degree at different depth.

In conclusion, the relative simplicity of our EOS allows us to obtain the analytic expressions for the heat capacity and the adiabatic gradient necessary to compute atmosphere models. This has the advantage of avoiding potential errors due to numerical inaccuracies that can be present in a sophisticated EOS, which requires interpolation from tables and numerical evaluation of derivatives.

## 3  Model atmosphere

Cooling of WDs is a very slow process so that the atmosphere structure responds to quasi-static boundary conditions. Furthermore, the region where the emergent radiation is formed has a thickness much less than the stellar radius. These properties allow us to choose a static model atmosphere with the approximation of plane-parallel layers and constant surface gravity throughout the atmosphere.

On the other hand, because the gas densities are very high in the layers where the emergent radiation is formed, local thermodynamic equilibrium (LTE) is a good approximation for the matter state and we adopt that.

We take the atmophere gas to be composed of H, H2, e-, H-, H+, H$`{}_{}{}^{+}{}_{2}{}^{}`$ and H$`{}_{}{}^{+}{}_{3}{}^{}`$. The main constituents are H2, H, H+ and e-, depending on the physical conditions. The other species are usually present as traces, but they can play an important role in the absorption and emission of radiation. In fact, the H- ion can be a dominant source of opacity in cool atmospheres, whereas H$`{}_{}{}^{+}{}_{2}{}^{}`$ and H$`{}_{}{}^{+}{}_{3}{}^{}`$ are important because their presence affects the concentration of H- and e-.

To determine the equilibrium concentrations we consider the reactions given in Table 1. The values of energy are those compiled by Lenzuni & Saumon (1992). We use the partition functions calculated by Irwin (1981) for H and H2, by Sauval & Tatum (1984) for H$`{}_{}{}^{+}{}_{2}{}^{}`$, and Neale & Tennyson (1995) for H$`{}_{}{}^{+}{}_{3}{}^{}`$. We have checked our chemical equilibrium against results from an independent code (Saumon 1999; based on an equation of ideal state, except in a sharp cuttoff in the H partition function depending on the distance to the nearest neighbor). Though some differences (at most 12%) appear between both calculations (mainly due to different ways of computing the H2 partition function), the agreement is satisfactory in general.

### 3.1  Opacities

We have consider the most important opacity sources in a hydrogen atmosphere. Bound-free and free-free absorption coeficients of H have been taken from Mihalas (1978), with Gaunt factors from Menzel & Pekeris (1935). Isotropic Thomson scattering is assumed for the opacity of free electrons. We adopted absorption coeficients of the negative hydrogen ion according to John (1988) for bound-free transitions, and analytic fits by Gray (1992) based on Bell & Berrington (1987) results for free-free transitions.

Collisions between H2 molecules yield collision-induced absorption (CIA), consisting of rototranslational and rotovibrational bands with peaks located in far-red and infrared spectral regions. We use H2-H2 CIA cross sections calculated by Borysow and coworkers, including the most recent results (Borysow, Jorgensen & Zheng 1997). For increasingly large densities, absorption from collisions involving more than two H2 molecules becomes important. Three-body collisions are no longer negligible for densities $`\rho >0.02`$ g cm-3 and they dominate the opacity above $`\rho >0.2`$ g cm-3 (Lenzuni & Saumon 1992). Ignoring terms of higher order, the total frequency-integrated CIA coefficient can be expressed as (Lenzuni, Chernoff & Salpeter 1991)

$`k=k_2\rho ^2+k_3\rho ^3,`$ (21)
where $`k_2`$ and $`k_3`$ account for absorption by binary and ternary collisions. We assume the ratio $`k_3/k_2=0.05/\rho `$ meassured experimentally by Hare & Welsh (1958).

Finally, cross section of Rayleigh scattering from molecular hydrogen is taken from Dalgarno & Willians (1966). We also consider other opacity sources such as Rayleigh scattering from H, bound-free and free-free absorption of H$`{}_{}{}^{+}{}_{2}{}^{}`$, and free-free opacity of H$`{}_{}{}^{}{}_{2}{}^{}`$. For all of them, we use coefficient expressions from Kurucz (1970).

### 3.2 Convective Transport

Energy transport by convection is treated within the mixing-length theory (e.g., Cox & Giuli 1968; Mihalas 1978) in most model atmosphere calculations of WDs. Although this phenomenological formalism involves free parameters, Bergeron, Wesemael & Fontaine (1992) have found that model results with $`T_{eff}`$ below $``$ 8000 K are insensitive to the assumed parameterization. In addition, they found that the predicted emergent fluxes show dependence on the efficiency of convection in the range $`T_{eff}`$ $``$ 8000-15000 K. However, it is difficult to decide which is the most appropriate parameterization for WDs. Preliminary analyses of the effects of convective efficiency on atmospheric parameters of ZZ Ceti 111WDs of DA class (i.e. with Balmer lines and no He I or metals present, and therefore, having H-rich atmospheres) showing non-radial pulsations. and other lukewarm DA stars, suggest that ML1 and ML3 parameterizations (following the nomenclature of Fontaine, Villeneuve & Wilson 1981) could be inadequate (Wesemael et al. 1991, Bergeron et al. 1992). In this context, the ML2 version of the mixing-length theory has been considered a viable alternative. The results presented here use this version of convection.

### 3.3  Numerical Method

The model atmosphere may be calculated using an iteration procedure over linearized structure equations in terms of pertubations of the temperature and density (or pressure). However, we adopt a linearization method for the temperature alone (Gustafsson 1971), the density and pressure structure being obtained from a given temperature distribution by integration of the hydrostatic equilibrium equation. This procedure yieds similar results to a complete linearization method, but has a much shorter computing time (Saumon et al. 1994; SBLHB). Additional economy of execution time is obtained from formulating the radiative transfer in terms of moment equations, with closure using variable Eddington factors (Auer & Mihalas 1970).

To calculate a model atmosphere, the equations of radiative transfer and constant flux condition are linearized as functions of the temperature. The equations are discretized and organized following the Rybicki scheme (for details see Gustafsson 1971 and Gustafsson & Nissen 1972). The calculation of a model starts with a given temperature distribution from which the hydrostatic equilibrium equation is integrated, and the opacities and thermodynamic quantities derived. Then, transfer and energy constrain equations are solved and we obtain a temperature correction. The whole process is iterated to convergence. The Eddington factors are evaluated in each step from a Feautrier solution of the radiative transfer, using known values of the source function.

### 3.4  Model convergence properties

We remark on some points about the performance of the calculations. Although the correction procedure adopted here has proved effective and stable, certainly it is not as stable as when the energy transfer is purely radiative. Numerical difficulties can appear to solve the coupled equations of transfer and conservation energy when both convective and radiative modes are present (Gustafsson 1971; Grenfell 1974; Koester 1980; Bergeron et al. 1991; SBLHB). These difficulties are not exclusive to partial or complete linearization methods, since they have also been reported by Mihalas (1965) using a modification of the temperature-correction method proposed by Avrett and Krook (1963).

We found that numerical instability can occur mainly in atmospheres of WDs with low Teff and high gravity. In these models, the convection becomes extremely efficient and, therefore, the temperature gradient is nearly adiabatic ($`_A<10^3`$ in convective deep layers). Then, the temperature correction oscillates yielding a large convective flux in one iteration, and zero convection in that following for the same layer, according to the Schwarzschild instability criterion. Among the different ways that have been proposed to stabilize the iterations (see above references), we found convinient a similar method to the one suggested by Koester (1980). This consists of setting $`F_c=FF_{rad}`$ instead of $`F_c=0`$ when $`<_A`$ for a layer which would be convective in the converged model.

On the other hand, a good initial temperature distribution decreases the number of iterations necessary and increases the range of convergence of the model in the ($`T_{eff}`$, log $`g`$) plane. For purely radiative models, one can start with an approximate solution of the gray atmosphere problem, e.g. the Eddington distribution (Mihalas 1978). However, when the convection is present, this class of distributions usually leads to an enormous convective flux which may have destabilizing consequences on the temperature correction method. A more suitable starting model is obtained when in the convective zone where $`F_c>F`$, the temperature gradient is changed and forced to yield $`F_c=F`$ (Nordlund 1974). This procedure yields a suitable initial temperature stratification in convective models, from which to apply the rigorous temperature correction scheme.

Finally, as in the report of SBLHB, we found that a discontinuity develops in the temperature distribution of converged models with T$`{}_{eff}{}^{}`$ 4400-5500 K at $`\mathrm{log}g7.5`$, extending to lower effective temperatures when the gravity decrease. Such as discussed further by Bergeron et al. (1995a), this effect is due to competition between opacities of H2 molecule and hydrogen negative ion. This convergence problem can be solved using an appropriate form of the radiative equilibrium equation (SBLHB), but unfortunately needs a large number of iterations.

## 4 Results and discussion

We have calculated a grid of models covering an extensive range of effective temperatures $`2500T_{eff}60000`$ K and surface gravities $`6\mathrm{log}g9`$. In these computations we used 76 depth points to solve the hydrostatic equilibrium and energy tranfer equations in the range $`6\mathrm{log}\tau 1.5`$. The continuum radiation field was computed at 87 to 107 frequency points depending on the effective temperature of the model. The points were also chosen taking into account the location of photometric bandpasses used in the colour computations (see Section 4.3).

The temperature-correction procedure was repeated until the temperature changes were everywhere $`\mathrm{\Delta }T/T<0.1\%`$, for most models. In hot models, $`T_{eff}14000`$ K, which are radiative or weakly convective, less than 10 iterations were sufficient to reach the convergence, yielding deviations from flux constancy less than 0.1 per cent everywhere. For lukewarm models, $`9000T_{eff}14000`$ K, convergence of $`\mathrm{\Delta }F/F<1\%`$ is obtained with about 10 iterations. As the convection becomes more efficient, $`T_{eff}9000`$ K, the number of iterations increases to about 40. In the range of temperatures where the onset of the H2 formation takes place, we need to arise the iteration number to approximately 150, due to the change of the boundary conditions in energy balance equations (see SBLHB).

### 4.1  Atmospheric structure

Strongly convective flux transport develops in the inner layers of cool and lukewarm WD atmospheres. This occurs preferentially in zones where hydrogen is partially ionized or dissociated, which yields a significant drop in the value of the adiabatic temperature gradient favouring the convective instability. Fig. 2 shows the fraction of flux carried by convection in atmospheres with surface gravities log $`g=6.5`$ and 8.0. It is also shown for log $`g=6.5`$ models the regions where H2, H and H+ in turn dominate the gas. Convection is absent from the atmosphere in the hotter early stages of a WD. According to our calculations, a superficial convection zone appears as the star cools below $`T_{eff}12000`$ K for WDs with low surface gravities (log $`g=6.5`$) and below $`T_{eff}16000`$ K for more massive stars (log $`g=8`$). The convective zone extends over the region where part of the emergent continuum radiation is formed ($`0.1\tau 1`$).

For cooler models, approximately $`T_{eff}<10000`$ K for log $`g=6.5`$ and $`T_{eff}<11000`$ K for log $`g=8.0`$, the efficiency of convection increases and a considerable fraction of the flux (more than 90 per cent) is carried by convection in the innermost layers ($`\tau 10`$). Radiation always dominates in the optically thin layers $`\tau 1`$. However, thanks to high densities and molecule formation, a little convective flux persists up to very transparent layers, $`\tau 10^3`$, in atmospheres around $`T_{eff}4000`$ K. The computed models predict a reduction of convection for very cool atmospheres $`T_{eff}2500`$ K, but this result may be affected by non-ideal contributions to the EOS. Slight pressure dissociation of H2 can be expected at large depths in these models, affecting the profile of the convective flux as consequence of changes in the adiabatic gradient (Saumon & Jacobson 1999).

For log $`g=8.0`$ models, the convection structure displayed in Fig. 2 is very similar to that showed by Bergeron et al. (1991, Fig. 2) in the common region ($`T_{eff}5000`$ K). As a difference, our results show as the molecular formation (not included in that work) strongly affects this structure below $`T_{eff}=6000`$ K, increasing the convection toward superficial layers.

Results about the run of several physical variables are shown in Figs 3-5 for selected models with different effective temperatures at log $`g=6.5`$ and 8.0. For the common models at log $`g=8`$, a good agreement is obtained from our $`T(P)`$ profiles compared to those by Bergeron, Wesemael & Fontaine (1995a, Fig. 2, $`T_{eff}=4000`$ \- 10000 K) and Saumon & Jacobson (1999, Fig. 1, $`T_{eff}=2500`$ and 3000 K). The only differences are found in the transparent layers of hot models ($`T_{eff}=8000`$ \- 10000 K), where the line blanketing (considered by Bergeron et al.) yields there slight changes in the temperature distribution.

Our sets of models with log $`g=6.5`$ and 8.0 present similar aspects except that higher pressures and densities are found in each model atmosphere with increasing $`g`$ (Figs 3-5). In terms of $`T_{eff}`$, it is possible to distinguish roughly three classes of atmospheric structures. One corresponds to hotter models, $`T_{eff}12000`$ K at log $`g=6.5`$ and $`T_{eff}14000`$ K at log $`g=8.0`$, which are purely radiative or with a superficial zone weakly convective. They show rather smooth temperature and density gradients without additional significant details. Throughout these atmospheres, most material is ionized.

As $`T_{eff}`$ decreases, a new atmospheric structure takes place. The models with $`T_{eff}`$ 4500-12000 K at log $`g=6.5`$ and $`T_{eff}`$ 4500-14000 K at log $`g=8.0`$, cover the range from the onset of convection with little flux transported by it to the models where the molecular formation becomes important. In the present group of models, the relationship between temperature and pressure exhibits (Fig. 3) a clear transition from a very small gradient in radiative outer layers to a steep gradient in the top of the convective zone. In contrast, the run of density with pressure becomes relatively shallow in the transition region between radiative and convective zones (Fig. 5). In this sense, the model $`(T_{eff},\mathrm{log}g)=(12000,6.5)`$ is interesting, because it shows a slight inversion of the density in approximately $`\mathrm{log}P=4.7`$. It can be interpreted as result of a too fast temperature rise with increasing $`P`$ (see Fig. 3), which is compensated by a decrease of the density ($`P\rho T`$).

Finally, the coolest models studied, $`2500T_{eff}4500`$ K, show some proper features. They are characterized by molecular formation which, as mentioned earlier, increases the convective efficiency and extends the convective zone toward very transparent regions (Figs 2 and 4), reaching a maximum at $`T_{eff}=4000`$ K. This effect has also been found in very low mass stars with log $`g=5`$ (SBLHB). At such low $`T_{eff}`$, our models show a uniform temperature in the upper radiative layers and a rather sudden rise with depth in still transparent regions (Fig 4). Deep inside the convective zone results a low temperature gradient, which is very close to the adiabatic value.

Fig. 5 shows that very cool atmospheres extend to extremely high density and pressure regimes (note that the top and bottom layers of all models are located at log $`\tau =6`$ and 1.5, respectively). Non-ideal effects would actually be very small in models with relatively low gravity log $`g=6.5`$, but it can become important in the deep layers (below $`\rho 10^2`$) of very cool models with log $`g=8.0`$. Nevertheless, the $`T(P)`$ structures of the atmosphere models are away from the region where Saumon & Chabrier (1992) predict the plasma phase transition of hydrogen (Fig. 3).

Compared to hot and lukewarm models analysed above (except at $`T_{eff}`$ around of the onset of convection), the temperature distributions of very cool atmospheres ($`T_{eff}4500`$ K) are particularly sensitive to the value of the gravity. The surface temperature decreases and the gradient in the ($`\tau ,T`$) plane increases as the gravity is raised (Figs 3 and 4). This is a consequence of the CIA opacity (Section 3.1) which, being a collisional process, becomes stronger in atmospheres with high gravity (therefore more dense) increasing the cooling of upper layers.

### 4.2  Emergent flux distribution

To illustrate the behaviour of the WD spectra in the ($`T_{eff},g`$) plane, we display in Figs 6a-6f the emergent flux corresponding to some selected $`T_{eff}`$ values for log $`g=6.5`$, 8.0 and 9.0. The importance of different opacity sources can be noted in Fig. 7, where individual contributions to the total monochromatic opacity are shown for some relevant physical conditions. These correspond to a Rosseland optical depth $`\tau =1`$ in atmospheres with log $`g=8`$ and $`T_{eff}=60000`$, 10000, 6000 and 2500 K (from 7a to 7d, respectively). Fig. 7d shows, with heavy lines, the opacity contributions at $`\tau =0.1`$ for the 2500 K model.

Emergent fluxes for the three highest values of $`T_{eff}`$ clearly show the effect of the atomic hydrogen absorption, which is a dominant source at such high temperatures (Figs 6a-c and 7a-b). Among the spectra present in Fig. 6, only the $`T_{eff}=60000`$ K models show an important radiation flux emerging from the Lyman continuum.

At $`T_{eff}=15000`$ K and log $`g=9`$, an important convective zone (with a maximum of 40% of the flux transported by convection) develops between log $`\tau 1.5`$ and 1.3. This strongly affects the Balmer continuum decreasing Lyman and Balmer jumps, with respect to lower gravity models (Fig. 6b). This sensitivity offers the possibility of an observational test of the convection theory even with relatively low-quality data.

Bound-free and free-free process of H- and atomic Rayleigh scattering start to take place as the effective temperature diminishes to 10000 K. The most prominent feature in the spectra of these atmospheres is the Balmer jump, which markedly weakens for large gravities. We note that these changes in the Balmer jump can also be reproduced by the presence of helium in the atmosphere and an increasing of its abundance (Wegner & Schulz 1981; Bergeron et al. 1991).

A large variety of opacity sources are found in the atmospheres of lukewarm WDs (Fig. 7c), the H- opacity being the most important. Since H- bound-free opacity is dominant and weakly frequency dependent, the emergent flux of $`T_{eff}=6000`$ K model is relatively featureless and resembles more a blackbody distribution (Fig. 6d). Furthermore, it is relatively insensitive to $`g`$.

Fig. 7d emphasizes the importance of the molecular opacity at $`\tau =1`$ in a ($`T_{eff},\mathrm{log}g)=(2500,8.0)`$ model. As one proceeds outward in the atmosphere, the opacity becomes completely dominated by CIA process and Rayleigh scattering (heavy lines in Fig. 7d, corresponding to $`\tau =0.1`$), due to a marked decrease of the abundance of all species but the H2 molecule. Besides, as the temperature decreases, separate vibrational bands are clearly identifiable in the CIA opacity. Note that the decrease of the density (from $`\tau =1`$ to $`\tau =0.1`$) does not affect the molecular Rayleigh scattering contribution to the total opacity per gram, because the H2 remains as the most abundant component.

Synthetic spectra obtained for $`T_{eff}=2500`$ and 4000 K (Figs 6e-f) show a well-known result: the effect of CIA opacity is largest for WDs of low $`T_{eff}`$ and high gravity, corresponding to a large abundance of H2 and high density, and, therefore, a large number density of H2-H2 pairs. The effect of the four vibrational bands of H2 (Fig. 7d) is clearly visible in the spectrum of models with $`T_{eff}=2500`$ K. This strong absorption over infrared spectral regions coincides with the peak emission of a blackbody at the same $`T_{eff}`$. As a consequence, the stellar flux is forced to emerge at shorter wavelengths, around the minimum in the opacity caused by the H2 CIA and the Rayleigh scattering.

### 4.3  Colour indices

From the emergent flux provided by the atmospheric model, we can derive magnitudes and colours. We have calculated broadband colour indices using the Johnson-Cousins BVRI and Johnson-Glass JHK response functions of Bessell (1990) and Bessell & Brett (1988), respectively. The magnitude $`m`$ corresponding to a bandpass with tranmission function $`S_\lambda ^m`$ was computed from

$`m=2.5\mathrm{log}\left[4\pi \left({\displaystyle \frac{R}{D}}\right)^2{\displaystyle \frac{_0^{\mathrm{}}H_\lambda S_\lambda ^m𝑑\lambda }{_0^{\mathrm{}}S_\lambda ^m𝑑\lambda }}\right]+c_m,`$ (22)
where $`H_\lambda `$ is the monochromatic Eddington flux from the model atmosphere for a WD of radius $`R`$ and distance $`D`$ from Earth. We adopted the calibration constants $`c_m`$ given in Bergeron, Ruiz & Leggett (1997, BRL), yielding magnitudes on the Carnegie image tube system. A set of colour indices is obtained from the difference between magnitudes, resulting independent of the radius and distance of the star. Fig. 8 shows the location of the filter transmitions over two illustrative emergent fluxes. We have also calculated the bolometric correction at the visual using the expression derived by Bergeron, Wesemael & Beauchamp (1995b, Eq.)

$`BC=2.5\mathrm{log}{\displaystyle _0^{\mathrm{}}}H_\lambda S_\lambda ^V𝑑\lambda 10\mathrm{log}T_{eff}+15.6165.`$ (23)
In view of increasing attention on very low mass WDs in the last few years, we report the colours and bolometric corrections for a subsample of our grid, corresponding to low gravity (log $`g=6.5`$) and cool ($`T_{eff}=4000`$-$`8000`$ K) models (Table 2). We also show the results for a canonical value $`\mathrm{log}g=8.0`$.

Fig. 9 shows the ($`VR,RI`$) two-colour diagram for all models calculated. The solid lines connect models with same gravity and log $`g=6.5`$ (0.5) 9.0 (from top to botton). Also displayed is a subset of 44 DA stars (filled circles) and 45 non-DA stars (open circles) observed by BRL and having complete BVRI-JHK photometric data 222The sample contains WDs from the catalog of McCook & Sion (1987) reclassified DA by BRL.. By comparison, we include the hydrogen atmosphere sequence at log $`g=8`$ of Bergeron (2000; diamonds), and that of Saumon & Jacobson (1999; triangles) covering $`T_{eff}4000`$ K. The cooling sequence of Richer et al. (2000; dashed line) for a WD with carbon-oxygen core, hydrogen atmosphere and mass 0.6 $`M_{}`$ is also shown.

Our models predict that all stars with hydrogen atmospheres above $`T_{eff}4500`$ K, fall on a line of almost constant surface gravity with little dispersion. The DA stars plotted in the figure remain very close to the synthetic colour sequences over the range $`3500T_{eff}8500`$ K. In contrast, the non-DA stars reveal more dispersion, although most of them seem to follow the hydrogen sequences. Only a few non-DA stars (marked on the figure) show colour indices strongly affected by other non-hydrogen sources of opacity. Among them, the very low luminosity WD LHS 3250 has a DC classification (continuous spectrum without lines deeper than 5%) and its atmospheric composition is not known (Harris et al. 1999). BPM 27606 is a DQ star (ie. with carbon features) whose spectrum presents very strong C2 Swan bands (BRL). LP $`79029`$ is also a DQ star.

Fig. 9 shows that the theoretical sequence at log $`g=8.0`$ of Bergeron (2000) is well reproduced by our models (this occurs also for other log $`g`$ values, not shown here). For the coolest atmospheres ($`T_{eff}<3500`$ K), our sequence at log $`g=8`$ becomes bluer than that of Saumon & Jacobson (1999). Nevertheless, the three independent models found an equivalent location for the turnoff in $`RI`$. This turnoff is a consequence of the onset of the H2 formation (which starts the CIA and molecular Rayleigh opacities), and therefore it varies with the surface gravity. The largest values of $`RI`$ correspond, at instance, to $`T_{eff}=3500`$ K at log $`g=6.5`$ and $`T_{eff}=4000`$ K at log $`g=9.0`$. In agreement with Saumon & Jacobson, we found that several colour indices show the presence of a turnover at low effective temperatures, with $`BV`$ being an exception (see below).

In addition, the cooling curve of Richer et al. (2000) follows the other sequences of models and the sample observed for $`T_{eff}>4000`$ K, but crosses several sequences with different $`g`$ for lower $`T_{eff}`$ where a constant or smoothly increasing gravity should be expected as the star cools. We also note a certain irregular behaviour where our grid predicts that colours are insensitive to changes in the gravity ($`T_{eff}>4500`$ K).

The ($`BV,VK`$) two-colour diagram is displayed in Fig. 10. There is a close agreement between Bergeron’s and our models at log $`g=8.0`$, except for the hot models. A similar turnover in $`VK`$ colour (reported earlier by Mould & Liebert 1978) is predicted by both sets of models. Nevertheless, significant departures of our sequence from that of Bergeron can be observed below $`BV0.3`$ ($`T_{eff}>7500`$ K). These discrepancies presumably have their origin in the hydrogen line opacity and the occupation probability formalism (Hummer & Mihalas 1988) used by his code and not included in our calculations.

The diagram also reveals a certain difficulty of the theory in reproducing the observed values in the $`BV`$ colour at low effective temperatures. In fact, there is an increasing difference between observed and theoretical colour indices above $`BV0.7`$. BRL have interpreted this in terms of a pseudo continuum opacity at the Lyman edge (not taken into account in the atmosphere codes), which affects the flux in the $`B`$ magnitude (for details see their paper).

Fig. 11 illustrates that the sample of observed stars in the ($`VI,VK`$) diagram is well reproduced by our model grid, which is also in good agreement with the theoretical sequence of Bergeron. It is also observed that a marked sensitivity of $`VK`$ and $`VI`$ colours to the effective temperature and surface gravity in models beyond the turnover. The reason of this is that $`I`$ and $`K`$ bandpasses lie just where the H2 CIA has strong contributions from two major vibrational transitions, while the $`V`$ bandpass receives a flux excess that is forced to emerge in the optical spectral range. One would expect that this sensitivity in the WDs with hydrogen atmospheres could be very useful for cosmochronological purposes.

### 4.4 White dwarfs with low mass

Helium core white dwarf (He WD) stars are considered the final product of the evolution of some close binary systems. The mass of these objects should be smaller than that required for degenerate helium ignition, $`0.5`$ $`M_{}`$ (Mazzitelli 1989), thus, resulting very low mass WDs. They have been detected in large surveys (Bragaglia et al. 1990; Bergeron, Saffer & Liebert 1992; Bragaglia, Renzini & Bergeron 1995; Saffer, Livio & Yungelson 1998), concluding that about 10% of presently known WDs are He WDs.

We have used our atmosphere code to study the changes in colours and absolute visual magnitude of He WDs during the cooling. The calculations were carried out by using the cooling sequences of Benvenuto & Althaus (1998) for helium-core composition and hydrogen envelopes. The chosen evolutionary models cover stellar masses $`M=0.15`$, 0.30 and 0.45 $`M_{}`$ with $`M_H/M=4\times 10^3`$, $`2\times 10^3`$ and $`4\times 10^4`$, respectively, where $`M_H`$ is the mass of the H envelope. We wish to point out that recent studies (Driebe et al. 1998; Sarna, Antipova & Muslimov 1998; Sarna, Antipova & Ergma 1999) predict the formation of He WDs with massive hydrogen envelopes. Since the H-burning increases with $`M_H`$, thus implying longer evolutionary times, the ages estimated in this paper (using models with rather thin hydrogen envelopes) should be considered as lower limit values.

We summarized here the most important results. Interesting observations can be inferred from the analysis of Fig. 12, which shows evolutionary tracks along with an observational sample of cool stars in ($`JH,VK`$) and ($`JK,RI`$) two-colour diagrams. The cooling tracks in both diagrams show an almost linear behaviour following the observed photometry up to a predicted effective temperature of approximately 4500 K. Although there is an important scatter of observed colours, it is clear that $`JH`$ and $`JK`$ show a shifting to the blue in the cool extreme of the sample, for DA and non-DA stars. This effect is very well represented by the theoretical cooling sequences, which gets dramatically bluer in these colour indices due to the onset of the H2 collision-induced absorption for $`T_{eff}4500`$ K. Note how the theoretical sequences follow the red edge in $`VK`$ and $`RI`$ of the observational sample, with the less massive models on the right.

Fig. 12 seems to provide observational evidence of the presence of a strong pressure-induced opacity in most of the cool atmospheres. It also suggests that very cool WDs identified as non-DA, could have an important abundance of hydrogen in their atmospheres (perhaps beyond the spectroscopic limit of detection due to the low temperatures). A good example is LHS 1126 (see Fig. 12), the peculiar emergent flux has been interpreted in terms of H2-He CIA in a mixed H/He atmosphere with $`\mathrm{log}N`$(He)$`/\mathrm{log}N`$(H)$`=0.8`$ (Bergeron et al. 1994).

Another feature exhibited in this figure is the prediction of a second turnover in the $`JH`$ colour (around $`VK=0.3`$) for an advanced stage of the cooling. In this part of the ($`JH,VK`$) two-colour diagram, the tracks of stars with different mass coincide, but the cooling times differ significantly. For example, the lowest value of $`JH`$ is reached by a 0.15 $`M_{}`$ He WD at 7 Gyr, while a 0.45 $`M_{}`$ star needs 16 Gyr more, due to its lower radius and, therefore, to a lower rate of cooling.

In the ($`VI,M_V`$) colour-absolute magnitude diagram of Fig. 13, the tracks corresponding to masses 0.15, 0.30 and 0.45 $`M_{}`$ are shown from right to left. The diagram shows that when a 0.15 $`M_{}`$ He WD is older than 5 Gyr, it gets bluer in the $`VI`$ colour. There, the absolute magnitude changes slowly and the star remains brighter than $`M_V=16.5`$ up to 10 Gyr. Similar behaviour is found for more massive models, but a longer time (around 20 Gyr at $`M=0.45`$ $`M_{}`$) is spent by these stars before the turnover in $`VI`$ takes place, as a consequence of an abundant H2 formation in their atmospheres.

A number of observed WDs appear redward in $`VI`$ of the 0.45-$`M_{}`$ sequence. Since the stellar evolution theory predicts that single WDs with very low mass cannot be formed within the lifetime of the Galaxy, those stars in the diagram correspond to unresolved (or barely resolved) double degenerate systems which appear overluminous, or single low mass WDs forming a binary system with a much less luminous companion.

We have used our code to derive atmospheric parameters, effective temperature and surface gravity, for five of these stars using photometric data of BRL. Combining these with evolutionary results we have also derived masses and ages for the stars. The effective temperature can be constrained from atmosphere fits to the observed colour indices. However, the surface gravity is poorly determined by this procedure for hot stars ($`T_{eff}4500`$ K), because the colours are rather insensitive to changes of $`g`$ in such models. Nevertheless, it is possible to derive the gravity and the mass of a star when its absolute magnitude is known. Therefore, we determine the effective temperature, the surface gravity and the mass of each star, fitting simultaneously its absolute visual magnitude and colours with a least-squares method. The $`BV`$ colour was not included in this procedure because there is an opacity source affecting the $`B`$ filter, which is not considered in the model (see Section 4.3). The mass-radius relations employed in the fitting and the evolution ages were taken from models kindly provided by Benvenuto & Althaus (2000, private communication). All interior models used here have hydrogen envelopes.

Treating the observed objects as single WDs, we obtain the results summarized in Table 3. Successive columns give the WD identifier, name, effective temperature, surface gravity, estimated mass, mean deviation of the magnitude and colour fitting, calculated and observed absolute visual magnitude, and age. The quoted errors have been determined as the changes necessary to produce a mean incertainty of 2$`\sigma `$ provided by the fitting procedure.

The effective temperatures determined here are on average 65 K hotter than those previously estimated by BRL. These differences are not significant against to the mean errors of both calibrations (120 K here and 110 K in the BRL work). As discussed above, the determination of the $`T_{eff}`$ is mainly provided by the broad-band fit, whereas the temperature uncertainties can be significantly constrained by the absolute magnitude fit.

Our $`\mathrm{log}g`$ values are 0.1 larger than the results of BRL, and we also obtain systematically larger mass for these stars (by $`0.06`$ $`M_{}`$), in part due to the use of different mass-radius relations. Furthermore, using the helium core models of Benvenuto & Althaus, we estimated ages about twice than the BRL determinations, who used carbon core models of Wood333Evolutionary models of helium core WDs are available only recently.. These differences have an origin mainly in that the heat content stored per gram in a helium core is larger than that of carbon (or carbon-oxygen) core for fixed temperature. Consequently, He WDs remain bright for comparatively longer times, while the C and C-O WDs cool more rapidly.

Among the five stars analysed here, only L$`8702`$ is confirmed as binary. This star is a double-lined spectroscopic binary composed of a detached pair of DA WDs (Saffer, Liebert & Olszewski 1988). LHS 239 is a component of a resolved binary (e.g. Poveda, Herrera, Allen, Cordero & Lavalley 1994). Its companion is LHS 240, a WD classified DC9 in the McCook & Sion catalog (1987). For LHS 239, BRL have found a discrepancy between the effective temperature derived from its energy distribution (colour-fit) and that obtained from spectroscopy (fit to H$`\alpha `$ line profile). They consider that it could be explained if this object is itself an unresolved double system. The remain stars in the sample (G$`1412`$, G$`1878`$ and G$`1287`$) are suspected double degenerate binary (BRL) because their derived masses fall below 0.4 $`M_{}`$ and the current theories attribute to such WDs a binary origin.

Table 4 shows the results of interpreting this group of stars as unresolved binaries containing two identical WDs. The combined mass estimated for these systems is below the Chandrasekhar limit, although the individual stars are too massive to have helium cores, except perhaps in the case of L$`8702`$. It is worth noting that combinations of two different atmospheres (not analysed here) could allow mixed systems of He WD + CO WD. The effective temperatures in Table 4 are equivalent to those deduced from a single WD fitting due to the lack of sensitivity to gravity variations of colours, except for LHS 239.

From a detailed spectroscopic analysis Bergeron, Wesemael, Liebert & Fontaine (1989) have deduced that L$`8702`$ is composed of a pair of DA WDs with comparable effective temperatures ($`T_1=7470`$ K and $`T_2=6920`$ K) and masses ($`M_1=0.47`$ $`M_{}`$ and $`M_2=0.52`$ $`M_{}`$). For these determinations, they adopted the zero-temperature mass-radius relation of Hamada-Salpeter (1961). Assuming identical components for the system and using both helium core and carbon-oxygen core models recently available, we confirm that their individual masses are within the core-helium ignition limit ($`0.450.50`$ $`M_{}`$, see Table 4). Our result for the effective temperature ($`T_{eff}=7200`$ K) is intermediate between those deduced by Bergeron et al..

We note that LHS 239 cannot be modelled satisfactorily by a binary with identical components (the mean deviation of the fit, $`\sigma =0.13`$, is much larger than the uncertainty of the observed magnitude, $`\sigma _{obs}=0.03`$), while a good fit is obtained considering this star as a single object ($`\sigma =0.05`$). These results suggest that the photometry observed for LHS 239 corresponds to a single He-core WD. We consider that another origin, different to unresolved binary possibility, would be need to explain the discrepancy in the calibrations of $`T_{eff}`$ obtained by BRL (see above).

The photometric fits for G$`1412`$ and G$`1878`$ do not show any clear preference between the single WD or unresolved binary interpretations. In particular, we have not found a very satisfactory hydrogen atmosphere fit to the observed colour indices for these stars, perhaps due to a significant presence of helium or metals in their atmospheres. For G$`1287`$, a good fit is obtained for both single and double star possibilities. Unfortunatly, for the deduced effective temperature (5600 K), the energy distribution is gravity-insensitive and the fits do not highlight the binary nature of this star. Additional spectroscopic data and line calculations could improve the energy distribution fits and the identification of the components of these and other (suspected) unresolved binaries.

## 5 Conclusions

We have computed a grid of hydrogen model atmospheres for white dwarfs covering an extensive region in the ($`T_{eff},g`$) plane. For this, we have employed a partial method of linearization of equations for stellar atmospheres in which the energy is transported by radiation and convection. In view of numerical instabilities usually found to solve the resulting set of equations, we use a suitable initial temperature distribution which improves the convergence properties taking into account the low values of temperature gradient expected in deep layers of convective models.

From a relatively simple thermodynamic model, we have derived analytical expressions for the adiabatic temperature gradient and heat capacities of a mixture of hydrogen species. They result appropriate and easy to apply in the calculation of the convective flux. In particular, for the density regimes present in most of the atmospheres studied, our expression of the adiabatic gradient successfully reproduces the values of the ideal equation of state developed by Saumon & Chabrier.

We have analysed in detail the structure and emergent flux of hydrogen white dwarf atmospheres over a wide range of effective temperatures and gravities. This analysis complements the most recent studies, to provide an overall picture of these atmospheres and the changes to them expected during the cooling. In particular, we find strong modifications of the emergent flux for $`T_{eff}4500`$ K produced by the pressure-induced absorption of molecular hydrogen. The results obtained from our code have been compared to those recently published. Our calculations confirm the importance of infrared photometry for the study of cool white dwarfs, and its use in the determination of ages of the Galaxy and globular clusters.

Spectral evolution of helium core white dwarfs has briefly been studied using evolutionary sequences recently available for these objects. Colour-colour diagrams show that the theoretical cooling sequences follow the low mass boundary at the cool extreme of an observational sample of white dwarfs. It is very interesting to note that on a timescale less than 9 Gyr, very low mass He WDs ($`0.15`$ $`M_{}`$) show strongly non-blackbody colours with absolute visual magnitude that are still below 16.5, while massive He WDs ($`0.45`$ $`M_{}`$) remains brighter than $`M_V=17`$ for 25 Gyr.

Combining atmosphere calculations with interior model results, we determined the possible fundamental parameters of some low mass white dwarfs. Because, in general, very low mass WDs are suspected to be close binaries, the estimated parameters depend on the possible components of each system. With the assumption that both components are identical, the analysed stars present preferentially carbon-oxygen core compositions but with combined masses below the Chandrasekhar limit.

We will present self-consistent cooling calculations in the near future, by incorporating our atmosphere code as an appropriate surface boundary condition for evolutionary models of white dwarfs.

Acknowledgments

I am deeply grateful to O. Benvenuto and L. Althaus for encouraging me to perform this study. I thank them for many useful and stimulating discussions, and their thorough reading of the manuscript. L. Althaus has provided me with invaluable help for this investigation. His suggestions and contributions resulted in substantial improvements to the atmosphere code and the results presented here. I am also grateful to D. Saumon, who offered significant comments and provided me with his abundance and colour calculations for comparison. It is a pleasure to thank P. Bergeron for kindly making available his colour calculations and A. Borysow for providing the updated cross-sections for the collision-induced absorption by H2. Many thanks are due to N. Cepeda and M. Gomez for kindly revising the English text. I wish to thank also the referee for his useful remarks which improved the earlier version of this paper.

* Figure 1.- Comparisons of values for the adiabatic temperature gradient $`_A=\mathrm{ln}T/\mathrm{ln}P|_S`$, obtained from SCVH EOS (dotted lines) and our expression (Eqs. (6-8), solid lines) calculated for different atmosphere models at log $`g=8`$ (with effective temperatures indicated on the plot).
* Figure 2.- The structure of the convection zones for models with log $`g=6.5`$ and $`8.0`$. Curves show the layers where convection carries 1, 10 (20) 90 and 95 per cent (from bottom to top) of the total flux at log $`g=6.5`$ (solid lines) and 8.0 (dashed lines). The dotted curves show from left to right 95, 50 and 5% of H2 abundance, and 5, 50 and 95% of H+ abundance for models at log $`g=6.5`$.
* Figure 3.- Temperature-pressure stratifications for atmosphere models with log $`g=6.5`$ (solid lines) and 8.0 (dotted lines). Models shown have $`T_{eff}=2500`$ (250) 4000, 4500, 5000 (1000) 12000, 15000 (5000) 30000 (10000) 60000 K (from bottom to top). For models at log $`g=6.5`$, the triangles indicate the $`\tau =2/3`$ level, and the circles show the regions where convection carries 1% of the total flux near the transition between radiative and convective zones. The plasma phase transition predicted by Saumon & Chabrier (1992) is represented by the heavy curve labeled “PPT”. This figure extends the results in Fig. 2 of Bergeron et al. (1995a).
* Figure 4.- Runs of temperature with Rosseland mean optical depth for models presented in Fig. 3. Solid circles indicate the top of the convective zone of models at log $`g=6.5`$.
* Figure 5.- Density profiles for models as those referred to in Fig. 4. Visible extremes of each curve are located at log $`\tau =6`$ (botton) and 1.5 (top).
* Figure 6.- Emergent flux distributions for various effective temperatures (indicated on the plot) and log $`g=6.5`$ (solid lines), 8.0 (dashed lines) and 9.0 (dash-dotted lines). Those are compared with a blackbody spectrum at the same $`T_{eff}`$ (dotted lines).
* Figure 7.- In an obvious notation, individual contributions to the opacity at a depth $`\tau =1`$ of the models with log $`g=8.0`$ and $`T_{eff}=60000`$, 10000, 6000 and 2500 K, from $`(a)`$ to $`(d)`$. These correspond to $`(\mathrm{log}T,\mathrm{log}\rho )=`$ ($`4.76,6.53`$), ($`4.03,5.75`$), ($`3.80,4.15`$), ($`3.53,1.13`$), respectively. Heavy lines in (d) indicate the most important opacity contributions at $`\tau =0.1`$ of the 2500 K model, where $`(\mathrm{log}T,\mathrm{log}\rho )=(3.09,2.46)`$.
* Figure 8.- Normalized profiles of emergent fluxes for $`T_{eff}=2500`$ (dashed line) and 10000 K (solid line) at log $`g=8`$, and the blackbody profile at 2500 K (dash-dotted line). Dots indicate the filter transmissions of $`BVRI`$ Johnson-Cousins and $`JHK`$ Johnson-Glass systems.
* Figure 9.- The ($`VR,RI`$) two-colour diagram shows hydrogen model atmospheres of white dwarfs. Heavy solid lines represent our models at log $`g=6.5`$ (0.5) 9.0 from top to bottom, while thin lines connect those with the same $`T_{eff}`$ (labeled). It is also shown the log $`g=8.0`$ sequences of Saumon & Jacobson (1999; triangles) and Bergeron (2000; diamonds); and the cooling track for a 0.6-$`M_{}`$ carbon-oxygen WD calculated by Richer et al. (1999; dashed line). Photometric data of DA (filled circles) and non-DA (open circles) stars are taken from BRL, except the colour indices of LHS 3250 which are taken from Harris et al. (1999).
* Figure 10.- The ($`BV,VK`$) colour-colour diagram shows the photometric data of BRL (circles), the log $`g=8.0`$ sequence of Bergeron (diamonds) and our model grid (lines). Notation is as Fig. 9.
* Figure 11.- Same as Fig. 10 for the ($`VI,VK`$) diagram.
* Figure 12.- Evolutionary tracks of He WDs with masses 0.15, 0.30 and 0.45 $`M_{}`$ (solid lines from right to left) displayed in ($`JH,VK)`$ and ($`JK,RI`$) two-colour diagrams. A subset of cool WDs observed by BRL is represented with filled (DA) and open circles (non-DA WDs). Ages in Gyr are labeled for 0.15 and 0.45-$`M_{}`$ models. Dashed lines indicate the cooling track of a black body (independently of its size!).
* Figure 13.- The colour-absolute magnitude diagram shows the cooling tracks of WDs with helium cores, hydrogen atmospheres and masses 0.15, 0.30 and 0.45 $`M_{}`$ (solid lines from right to left). Ages are indicated in units of 109 years. Dotted lines show the same cooling tracks assuming stellar radiation as black body. The observed photometry is taken from BRL (DA WDs: filled circles; non-DA WDs: open circles); Monet et al. (1992; diamonds) and Harris et al. (1999; LHS 3250).


format: text

text:
# Optimal conditions for observing Josephson oscillations in a double-well Bose-gas condensate

## I Introduction

In this paper we consider a Bose-Einstein condensate (BEC) at zero temperature trapped in a double-well potential, which is created by superimposing a harmonic trap and a Gaussian barrier. With a finite barrier height that is greater than the condensate chemical potential, the Gross-Pitaevskii dynamics of the condensate can be well described as a superposition of right and left localized condensates, the so-called two-mode ansatz Raghavan1999b . The weak link allows Josephson tunneling of condensate atoms between wells when there is a chemical potential difference (or a difference in phases) between the wells.

A double well system was implemented in the experiment described in Ref. Andrews1997b , where the Gaussian barrier was created by focusing a blue-detuned far-off-resonant light sheet into the center of a magnetic trap. The trap had a “cigar” geometry, that is $`\lambda 1`$, where we have defined the aspect ratio of the radial and axial trap frequencies $`\lambda \omega _\rho /\omega _z`$, which turns out to play an important role in our calculations. In the experiment of Ref. Andrews1997b , an *infinite* barrier along the $`z`$-axis was used to separate two independent condensates of about five million 23Na atoms and then removed to allow the condensates to interfere. In this paper we focus on Josephson tunneling through a *finite* barrier, which was not reported in Ref. Andrews1997b . To date, there have been no reports of experimental observations of Josephson tunneling in a double-well trap. Josephson tunneling of a condensate in a one dimensional optical lattice was reported in Ref. Anderson1998a . Each well in the array held a condensate of about one thousand 87Rb atoms. The chemical potential offset due to gravity caused the condensates to tunnel out of the wells. A crucial characteristic of the experiment is that the wells had an extremely “pancake” geometry ($`\lambda 1`$) with $`\lambda 0.002`$ Anderson1998a ; AndersonPrivate , in contrast to the “cigar” geometry of the trap used in Ref. Andrews1997b . A related system of two hyperfine states coupled by applied electromagnetic fields was reported in Refs. Hall1998b ; Matthews1999a ; Matthews1999b  for the case of strong coupling. Tunneling in a three-component condensate was reported in Ref. Miesner1999a .

The double-well tunneling system for a Bose condensate has been treated extensively in the recent theoretical literature Javanainen1986a ; Milburn1997a ; Zapata1998a ; Raghavan1999b ; Salasnich1999a . Nonlinear Josephson equations describing the oscillations of the relative population $`\eta (t)=(N_1(t)N_2(t))/N_c`$ and phase $`\varphi (t)=\varphi _1(t)\varphi _2(t)`$ between the wells have been thoroughly explored Raghavan1999b ; Marino1999a , and this analysis will be the basis of our paper. Here, the indices $`1`$ and $`2`$ label the left and right wells and $`N_c=N_1+N_2`$ is the fixed total condensate population, so that $`1\eta (t)1`$. A crucial property of this system is that, due to the nonlinear terms in the Josephson equations arising from the mean field interaction, there is a maximum amplitude of oscillation $`\eta _{\mathrm{max}}<1`$ less than unity. If the initial population difference $`\eta (0)\eta _{\mathrm{max}}`$, then very little tunneling occurs, a phenomenon referred to as macroscopic quantum “self-trapping” Raghavan1999b . The calculations of Ref. Salasnich1999a  were carried out using realistic parameters corresponding to the experiment of Ref. Andrews1997b , but used a finite barrier height to allow tunneling to occur. An extremely small value of $`\eta _{\mathrm{max}}0.01`$ was found, indicating that Josephson oscillations in that system would be very difficult to observe due to experimental error in measuring $`\eta (t)`$ (which we assume to be of the order of a few percent).

Because of the importance of having clear experimental demonstrations of quantum phase oscillations, we feel it is useful to explore a wider range of physical parameters in order to find optimal conditions for observing Josephson oscillations in a double-well. This is the main purpose of the present paper. The key quantities that determine $`\eta _{\mathrm{max}}`$ are the condensate density, the s-wave scattering length $`a`$, the atomic mass $`m`$, and the coupling strength between wells. In the absence of interactions ($`a=0`$), the problem reduces to a simple linear two-level system, so that $`\eta _{\mathrm{max}}=1`$. The effect of the mean field interaction can be reduced by reducing the density, scattering length, or the mass. In this paper, we do not vary $`m`$ or $`a`$, but rather, we consider two types of atoms, 23Na and 87Rb, each with a fixed mass and scattering length. We calculate $`\eta _{\mathrm{max}}`$ for a wide range of trap parameters by varying the trap geometry determined by $`\omega _z`$ and $`\omega _\rho `$, the barrier geometry, and the number of condensate atoms $`N_c`$. We find $`\eta _{\mathrm{max}}`$ is optimized by making the trap more “pancake” shaped ($`\lambda 1`$), which has the effect of maximizing the weak link area, while maintaining a sufficiently low condensate density. This result is consistent with the fact that Josephson tunneling was clearly observed in Ref. Anderson1998a  for an array of wells possessing an extreme “pancake” geometry. This behavior is in contrast to that for a trap with a “cigar” geometry, which minimizes the weak link area and compresses the condensate in two dimensions rather than one, resulting in a greater density.

## II Derivation of Model

At zero temperature, the condensate dynamics is well described by the Gross-Pitaevskii (GP) equation

$$i\mathrm{}\frac{\mathrm{\Phi }(𝐫,t)}{t}=\left[H_0+gn_c(𝐫,t)\right]\mathrm{\Phi }(𝐫,t),$$
(1)
where $`n_c(𝐫,t)=|\mathrm{\Phi }(𝐫,t)|^2`$ is the condensate density and

$$H_0=\frac{\mathrm{}^2}{2m}^2+U_{\mathrm{ext}}(𝐫).$$
(2)
Here, the external potential $`U_{\mathrm{ext}}(𝐫)=U_H(𝐫)+U_B(𝐫)`$ is created by superimposing a harmonic trap

$$U_H(\rho ,z)=m\omega _z^2(\lambda ^2\rho ^2+z^2)/2,$$
(3)
and a Gaussian barrier along the $`z`$ axis.

$$U_B(z)=U_0\mathrm{exp}(z^2/2\sigma ^2).$$
(4)
In our calculations, we consider a range of values for the trap frequencies $`\omega _z`$ and $`\omega _\rho `$, the barrier height $`U_0`$ and width $`\sigma `$, and the condensate population $`N_c`$.

### II.1 Two-mode ansatz

For a barrier height that is sufficiently larger than the condensate chemical potential

$$U_0>\mu _c,$$
(5)
there exist two nearly degenerate stationary solutions of the GP equation

$$\left[H_0+gn_c(𝐫)\right]\mathrm{\Phi }_\alpha (𝐫)=\mu _\alpha \mathrm{\Phi }_\alpha (𝐫),$$
(6)
where $`n_c(𝐫)=N_c|\mathrm{\Phi }_\alpha (𝐫)|^2`$. The index $`\alpha =\{S,A\}`$ indicates the reflection symmetry along the $`z`$-axis, symmetric or antisymmetric, respectively. These solutions satisfy the orthonormality condition $`𝑑𝐫\mathrm{\Phi }_\alpha (𝐫)\mathrm{\Phi }_\beta (𝐫)=\delta _{\alpha \beta }`$. Left and right localized states can be formed from these eigenstates $`\mathrm{\Phi }_1(𝐫)=\left[\mathrm{\Phi }_S(𝐫)+\mathrm{\Phi }_A(𝐫)\right]/\sqrt{2}`$ and $`\mathrm{\Phi }_2(𝐫)=\left[\mathrm{\Phi }_S(𝐫)\mathrm{\Phi }_A(𝐫)\right]/\sqrt{2}`$, respectively, which are also orthonormal.

In the two-mode approximation, the condensate dynamics is described by the following ansatz solution Raghavan1999b

$$\mathrm{\Phi }(𝐫,t)=\psi _1(t)\mathrm{\Phi }_1(𝐫)+\psi _2(t)\mathrm{\Phi }_2(𝐫).$$
(7)
Here, the complex coefficients $`\psi _i(t)=\sqrt{N_i(t)}\mathrm{exp}[i\varphi _i(t)]`$ are spatially uniform and contain all of the time dependence, while the two states $`\mathrm{\Phi }_1(𝐫)`$ and $`\mathrm{\Phi }_2(𝐫)`$ are localized in the left and right wells, respectively, and contain all of the position dependence. In a more accurate theory, one could account for the slow time evolution of the states $`\mathrm{\Phi }_i(𝐫)`$ due to the mean field interaction since the population $`N_i(t)`$ in each well evolves in time. If the amplitude of oscillation is relatively small, this can be shown to have a negligible effect Salasnich1999a ; Williams1999a .

### II.2 Gaussian variational solution

Since we are exploring the parameter space of $`\{\omega _z,\omega _\rho ,U_0,\sigma ,N_c\}`$ and must solve (6) many times, it is useful to construct a simplified model solution of the left and right states $`\mathrm{\Phi }_i(𝐫)`$. Several earlier studies have shown that a Gaussian variational solution gives an accurate description of the condensate Baym1996a ; Perez-Garcia1996a ; Menotti2001a . We therefore employ this procedure by taking the following ansatz for $`\mathrm{\Phi }_\alpha (𝐫)`$

$$\mathrm{\Phi }_\alpha (𝐫)=C_\alpha e^{A_\alpha ^2\rho ^2}\left[e^{B_\alpha ^2(z+z_\alpha )^2}\pm e^{B_\alpha ^2(zz_\alpha )^2}\right].$$
(8)
Here, $`C_\alpha `$ is a normalization constant given by

$$C_\alpha =\sqrt{\frac{2^{1/2}A_\alpha ^2B_\alpha }{\pi ^{3/2}\left[1\pm \mathrm{exp}(2B_\alpha ^2z_\alpha ^2)\right]}}.$$
(9)
The three variational coefficients $`\{A_\alpha ,B_\alpha ,z_\alpha \}`$ must be determined by minimizing the Gross-Pitaevskii energy functional

$$E(\mathrm{\Phi })=𝑑𝐫\left[\frac{\mathrm{}^2}{2m}|\mathrm{\Phi }|^2+U_{\mathrm{ext}}|\mathrm{\Phi }|^2+\frac{g}{2}|\mathrm{\Phi }|^4\right]$$
(10)
with respect to variations in $`\{A_\alpha ,B_\alpha ,z_\alpha \}`$. We note that these coefficients will be different for the symmetric and anti-symmetric states. The Gaussian form of our variational ansatz permits us to evaluate the integral in (10) analytically. However, due to the nonlinear terms in the resulting expression, we carry out the minimization of $`E(\mathrm{\Phi })`$ numerically. In Section III, we compare the Gaussian variational solution with the full numerical solution of (6) and find qualitative agreement over a wide range of parameters.

### II.3 Equations for the relative phase and population

We substitute the two-mode ansatz (7) into the GP equation (1), multiply by $`\mathrm{\Phi }_i(𝐫)`$ and integrate over position to obtain equations of motion for the complex coefficients $`\psi _i(t)`$ Raghavan1999b

$`i\mathrm{}\dot{\psi }_i(t)`$ $`=`$ $`\left[E_i+U_iN_i(t)\right]\psi _i(t)`$ (11)
$``$ $`\left[K+U_{1,12}N_1(t)+U_{2,12}N_2(t)\right]\psi _j(t),`$
where we have used the orthonormality property. Here, $`E_i`$ is the bare zero-point energy in well $`i`$,

$$E_i=𝑑𝐫\mathrm{\Phi }_i(𝐫)H_0\mathrm{\Phi }_i(𝐫)$$
(12)
and $`U_iN_i(t)`$ is the mean-field interaction energy, with

$$U_i=g𝑑𝐫\mathrm{\Phi }_i^4(𝐫).$$
(13)
The term $`K`$ in (11) is the coupling energy between condensates due to the finite probability of tunneling through the barrier,

$$K=𝑑𝐫\mathrm{\Phi }_1(𝐫)H_0\mathrm{\Phi }_2(𝐫).$$
(14)
In (11), we also include the mean-field contribution to the coupling, given by $`[U_{1,12}N_1(t)+U_{2,12}N_2(t)]`$, where

$$U_{i,12}=g𝑑𝐫\mathrm{\Phi }_i^2(𝐫)\mathrm{\Phi }_1(𝐫)\mathrm{\Phi }_2(𝐫).$$
(15)
As shown recently in Ref. Zhang2001a , these factors (15) give rise to a time-dependent coupling term, but were neglected in the calculation of Ref. Raghavan1999b . The integrals in (12)-(15) can be carried out analytically using the Gaussian variational solution (8).

It is straightforward to obtain the equations of motion for the phases and populations of each condensate by substituting $`\psi _i(t)=\sqrt{N_i(t)}e^{i\varphi _i(t)}`$ into (11),

$`\dot{N}_i(t)`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}}}[U_{1,12}N_1(t)+U_{2,12}N_2(t)`$ (16)
$`+`$ $`K][N_i(t)N_jt)]^{1/2}\mathrm{sin}[\varphi _i(t)\varphi _j(t)]`$
$`\mathrm{}\dot{\varphi }_i(t)`$ $`=`$ $`E_iU_iN_i(t)+[U_{1,12}N_1(t)+U_{2,12}N_2(t)`$ (17)
$`+`$ $`K\left]\right({\displaystyle \frac{N_j(t)}{N_i(t)}})^{1/2}\mathrm{cos}[\varphi _i(t)\varphi _j(t)].`$
Using (16) and (17), we can derive two equations for the relative phase $`\varphi (t)=\varphi _1(t)\varphi _2(t)`$ and the relative population $`\eta (t)=[N_1(t)N_2(t)]/N_c`$ Raghavan1999b

$`\dot{\eta }(t)`$ $`=`$ $`\omega _J(t)\sqrt{1\eta ^2(t)}\mathrm{sin}\varphi (t),`$
$`\dot{\varphi }(t)`$ $`=`$ $`\mathrm{\Delta }\omega \omega _C\eta (t)\omega _J(t){\displaystyle \frac{\eta (t)\mathrm{cos}\varphi (t)}{\sqrt{1\eta ^2(t)}}},`$ (19)
Here we have simplified the notation by defining the Josephson frequency

$$\omega _J(t)\left[2K+N_cU_{1,12}(1+\eta (t))+N_cU_{2,12}(1\eta (t))\right]/\mathrm{},$$
(20)
the “capacitive” frequency proportional to the mean field interaction

$$\omega _CN_c(U_1+U_2)/2\mathrm{},$$
(21)
and the frequency due to the difference in zero-point energies

$$\mathrm{\Delta }\omega \left[E_1E_2+(U_1U_2)N_c/2\right]/\mathrm{}.$$
(22)
For simplicity, we restrict our calculations to a symmetric double-well, since the overall behavior of $`\eta _{\mathrm{max}}`$ does not depend much on a relative offset in the well depths. In this case, the spatially averaged quantities (12), (13), and (15) are equal for each well ($`E_1=E_2`$, $`U_1=U_2`$, $`U_{1,12}=U_{2,12}U_{12}`$), allowing us to make the simplifications $`\mathrm{\Delta }\omega =0`$ and

$$\omega _J=2(K+N_cU_{12})/\mathrm{},$$
(23)
which is now time-independent. Using these results, the equations of motion (II.3) and (19) reduce to

$`\dot{\eta }(t)`$ $`=`$ $`\omega _J\sqrt{1\eta ^2(t)}\mathrm{sin}\varphi (t),`$ (24)
$`\dot{\varphi }(t)`$ $`=`$ $`\omega _C\eta (t)\omega _J{\displaystyle \frac{\eta (t)}{\sqrt{1\eta ^2(t)}}}\mathrm{cos}\varphi (t).`$ (25)
In Ref. Raghavan1999b , Raghavan *et al.* give an analytic solution of (24) and (25) in terms of Jacobian elliptic functions. The solution depends on the initial conditions for the relative number and phase, $`\eta (0)`$ and $`\varphi (0)`$. From the solution for $`\eta (t)`$ we can define two maximum amplitudes of oscillation, one for the case $`\eta (0)=0`$, the other in the case of a finite $`\eta (0)0`$. For equal populations initially, i.e. $`\eta (0)=0`$, the maximum amplitude of oscillation $`\eta _{\mathrm{max}}`$ can be obtained from equation (B5a) of Ref. Raghavan1999b  by taking $`\varphi (0)=\pi /2`$, which gives

$$\eta _{\mathrm{max}}=\frac{1}{\mathrm{\Lambda }}\left[2\left(\sqrt{\mathrm{\Lambda }^2+1}1\right)\right]^{1/2},$$
(26)
where $`\mathrm{\Lambda }\omega _C/\omega _J`$. In the non-interacting limit ($`g=0`$), $`\omega _C=0`$ and hence $`\mathrm{\Lambda }=0`$, so that (26) gives $`\eta _{\mathrm{max}}=1`$, as it must in this limit. In the strongly interacting limit $`\mathrm{\Lambda }1`$, (26) takes the approximate form $`\eta _{\mathrm{max}}\sqrt{2(\mathrm{\Lambda }1)}/\mathrm{\Lambda }1`$. If instead the initial population difference is nonzero, then a slightly higher amplitude $`\eta _c`$ can be obtained compared to that given in (26). Taking $`\varphi (0)=0`$ in equation (4.9) of Ref. Raghavan1999b , we obtain

$$\eta _c=2\sqrt{\mathrm{\Lambda }1}/\mathrm{\Lambda }.$$
(27)
If $`\eta (0)>\eta _c`$, then the solution is self-trapped. We note that in the strongly interacting limit $`\mathrm{\Lambda }1`$, these two quantities are approximately equal and satisfy the following relation $`\eta _c\sqrt{2}\eta _{\mathrm{max}}`$. We work with $`\eta _{\mathrm{max}}`$ in our calculations rather than $`\eta _c`$ because (26) is valid over the whole range of values for $`\mathrm{\Lambda }`$.

For small-amplitude oscillations, (24) and (25) can be linearized and combined to give

$$\ddot{\eta }(t)+\omega _J(\omega _C+\omega _J)\eta (t)=0.$$
(28)
This equation tells us that the system oscillates at the “Josephson plasma” frequency

$$\omega _{\mathrm{JP}}=\sqrt{\omega _J(\omega _C+\omega _J)}.$$
(29)
In the case where the mean field interaction is large compared to the coupling energy $`\omega _C/\omega _J1`$, this reduces to $`\omega _{\mathrm{JP}}=\sqrt{\omega _J\omega _C}`$. Due to the mean field interaction, the actual frequency of oscillation between the two wells can be large $`\omega _{\mathrm{JP}}>\omega _z`$, even though the coupling is weak $`\omega _J<\omega _z`$.

### II.4 Illustration of nonlinear oscillations

It is instructive to show the time evolution of $`\eta (t)`$ in order to visualize the behavior of the nonlinear Josephson oscillations and to clarify the precise meaning of $`\eta _{\mathrm{max}}`$, $`\eta _c`$, and $`\omega _{\mathrm{JP}}`$. In Fig. 1 we graph $`\eta (t)`$ for three different initial conditions $`\eta (0)`$ and $`\varphi (0)`$ for the solution of the coupled nonlinear equations (24) and (25). We first consider the case of equal populations $`\eta (0)=0`$ and an initial relative phase of $`\varphi (0)=\pi /2`$, shown by the solid line. The population difference $`\eta (t)`$ oscillates with an amplitude of $`\eta _{\mathrm{max}}`$. In the non-interacting limit (eg. $`\omega _C0`$), the problem reduces to a linear coupled two-level system and a complete population transfer between wells would occur in this case, $`\eta _{\mathrm{max}}=1`$, at a frequency of $`\omega _J`$. The amplitude suppression introduced by the mean field interaction arises from the term $`\omega _C\eta (t)`$ in (25), which acts as a time-dependent offset in the zero-point energies of the two wells. We also indicate the period of small-amplitude oscillation $`T_{\mathrm{JP}}2\pi /\omega _{\mathrm{JP}}`$ on the graph, which is valid when the amplitude is much less than $`\eta _{\mathrm{max}}`$ (we have verified this in calculations not shown here). $`T_{\mathrm{JP}}`$ gives a useful estimate of the period, which, in general, depends on the initial conditions.

We also consider the case of unequal initial populations $`\eta (0)0`$. For the dashed line in Fig. 1, we have taken $`\eta (0)=0.9\eta _c`$ and $`\varphi (0)=0`$, where the critical amplitude is $`\eta _c=0.11`$ for this case. If the initial population difference is larger than $`\eta _c`$, then the population becomes “trapped”, as shown for the case of $`\eta (0)=0.15`$ by the dot-dashed line. It is important to realize that in the non-interacting limit ($`\omega _C0`$), this self-trapping does not occur. In the absence of interactions, a full oscillation of the dot-dashed line would occur between $`\eta =\pm 0.15`$, at a frequency of $`\omega _J`$. We remark that for $`\eta (0)\eta _c`$, higher harmonics appear in the solution due to the nonlinearity Raghavan1999b . We refer the reader to Ref. Raghavan1999b  for a more detailed discussion of this system, where a full range of solutions of (24) and (25) is explored. We also note that this same behavior is also found in Ref. Williams1999a  for a condensate with two internal states coupled by a weak external electromagnetic field.

## III Results

In this section we calculate $`\eta _{\mathrm{max}}`$ as given by (26) as a function of the trap frequencies $`\omega _z/2\pi `$ and $`\omega _\rho /2\pi `$, the barrier height $`U_0`$ and width $`\sigma `$, and the condensate population $`N_c`$. To do this, we solve for the left and right localized states $`\mathrm{\Phi }_i(𝐫)`$ and then calculate the quantities (12)-(15), so that the frequencies $`\omega _J`$ and $`\omega _C`$ can then be computed. We compare two different solutions of $`\mathrm{\Phi }_i(𝐫)`$, a direct numerical solution Holland1997a  of the stationary GP equation (6) and the Gaussian variational ansatz (8). The variational coefficients $`\{A_\alpha ,B_\alpha ,z_\alpha \}`$ are computed numerically by minimizing the energy functional $`E(\mathrm{\Phi })`$ (10) using the Nelder-Mead simplex search algorithm in the software package MATLAB.

While varying the physical parameters $`\omega _z`$, $`\omega _\rho `$, $`\sigma `$, and $`N_c`$ in our calculations, we keep $`U_0`$ at a fixed height relative to the chemical potential $`\mu _S`$, which also varies. The bottom of the double well $`U_{\mathrm{min}}`$ is obtained by minimizing $`U_{\mathrm{ext}}`$ along the $`z`$-axis to give

$$U_{\mathrm{min}}=m\omega _z^2\sigma ^2\left[1\mathrm{ln}(m\omega _z^2\sigma ^2/U_0)\right].$$
(30)
In our calculations, we define the ratio $`\chi `$ of the barrier height $`U_0`$ to the chemical potential $`\mu _S`$ of the symmetric solution, each measured relative to the bottom of the well

$$\chi =\frac{(U_0U_{\mathrm{min}})}{(\mu _SU_{\mathrm{min}})}.$$
(31)
We emphasize that in *all* of our calculations, the coupling due to tunneling is *weak*, with $`0.05\omega _{\mathrm{JP}}/\omega _z0.5`$.

We first calculate $`\eta _{\mathrm{max}}`$ for a few typical magnetic traps shown in Table 1. For the “cigar” trap of Ref. Andrews1997b , taking $`N_c=5\times 10^6`$ and a barrier width of $`\sigma =6\mu `$m, the maximum amplitude of oscillation is miniscule $`\eta _{\mathrm{max}}=5.2\times 10^6`$. Even if we greatly reduce the condensate density by lowering the population $`N_c=10^4`$ and increase the coupling by reducing the barrier width $`\sigma =1\mu `$m, the maximum amplitude is still very small, with $`\eta _{\mathrm{max}}=0.036`$. The Josephson oscillations would not be visible in this case if the experimental error in measuring $`\eta _{\mathrm{max}}`$ were a few percent. For the “pancake” trap of Ref. Matthews1999a , taking a modest condensate population $`N_c=10^4`$ and a narrow barrier $`\sigma =2\mu `$m, we find $`\eta _{\mathrm{max}}=0.054`$. A better result can be obtained by decreasing the barrier width $`\sigma =1\mu `$m, which gives $`\eta _{\mathrm{max}}=0.14`$. Finally, we consider the spherical trap of Ref. Matthews1999b , which is a much shallower trap than the other two, so that the density of the condensate is much less. Taking $`N_c=10^4`$ and $`\sigma =2\mu `$m, we find $`\eta _{\mathrm{max}}=0.12`$, and for $`\sigma =1\mu `$m, $`\eta _{\mathrm{max}}=0.26`$. For the optimistic parameters of this shallow trap, Josephson oscillations should be observable.

For comparison we also consider the experimental setup of Ref. Anderson1998a , where a condensate was loaded into a one dimensional optical lattice. Each well in the array can be modeled as being harmonic as a first approximation, for which one obtains $`\omega _z/2\pi 7200`$ Hz and $`\omega _\rho /2\pi 17`$ Hz Anderson1998a ; AndersonPrivate . For a rough comparison, we consider a double well with this geometry. Taking $`\chi =1.5`$, $`\sigma =1\mu `$m, and $`N_c=10^4`$ 87Rb atoms, we find $`\eta _{\mathrm{max}}=0.85`$. This value increases to $`\eta _{\mathrm{max}}=0.99`$ if we reduce the population to $`N_c=10^3`$.

Motivated by these results, we obtain a very crude estimate of the scaling of $`\eta _{\mathrm{max}}`$ with the parameters. In the strongly interacting limit $`\mathrm{\Lambda }1`$, (26) reduces to

$$\eta _{\mathrm{max}}\frac{1}{\sqrt{\mathrm{\Lambda }}}\sqrt{\frac{\omega _J}{\omega _C}}.$$
(32)
Taking a simple WKB approach to calculate the tunneling probability Zapata1998a ; Sols1999a , an approximate scaling of the Josephson frequency can be obtained $`\omega _J/\omega _z\mathrm{exp}[\sqrt{m\sigma ^2(U_0\mu )/\mathrm{}^2}]`$. One can also get the approximate scaling of $`\omega _C`$ on $`m`$ and $`a`$ by taking $`\omega _C\mu _{\mathrm{TF}}`$, where $`\mu _{\mathrm{TF}}/\mathrm{}\omega _z=0.5[15\lambda ^2N_ca/z_{\mathrm{ho}}]^{2/5}`$ is the chemical potential of a condensate in a harmonic trap in the Thomas-Fermi limit Dalfovo1999a . Inserting the definition of the oscillator length, $`z_{\mathrm{ho}}=\sqrt{\mathrm{}/m\omega _z}`$, the expression can be simplified to $`\omega _C/\omega _z(ma^2N_c^2\lambda ^4\omega _z/\mathrm{})^{1/5}`$. Inserting these crude estimates of $`\omega _J`$ and $`\omega _C`$ into (32), we find that $`\eta _{\mathrm{max}}`$ scales like

$$\eta _{\mathrm{max}}\frac{\mathrm{exp}[\sqrt{m\sigma ^2(U_0\mu )/4\mathrm{}^2}]}{(ma^2\lambda ^4N_c^2\omega _z/\mathrm{})^{1/10}}.$$
(33)
This indicates that the amplitude of oscillation can be increased by decreasing $`m`$, $`a`$, $`N_c`$, $`\sigma `$, or $`U_0`$. For a fixed $`\omega _z`$, $`\eta _{\mathrm{max}}`$ increases as $`\lambda `$ is decreased, while for a fixed aspect ratio $`\lambda `$, $`\eta _{\mathrm{max}}`$ decreases when the trap is tightened by increasing $`\omega _z`$. Expression (33) only gives a very crude estimate of the scaling behavior. In the following discussion, we explore the dependence of $`\eta _{\mathrm{max}}`$ on the physical parameters more quantitatively by calculating $`\omega _J`$ and $`\omega _C`$ numerically.

In Figs. 2 - 5 we plot $`\eta _{\mathrm{max}}`$ (given by (26)) as a function of the trap geometry (by varying $`\omega _z`$ and $`\omega _\rho `$), $`N_c`$, and the barrier height $`\chi `$ and width $`\sigma `$. In each figure, we show results for the Gaussian variational solution (lines) and the direct numerical solution of the stationary GP equation (6) (filled points). The Gaussian variational solution agrees with the full solution within a few percent over a broad range of parameters. The greatest discrepancy occurs in Fig. 2 as $`\omega _z`$ decreases at fixed $`\omega _\rho `$ and in Fig. 4 as $`\sigma `$ increases, with the two solutions differing by a factor of 2 at $`\sigma =4\mu `$m for 87Rb. In all cases, the Gaussian solution slightly underestimates $`\eta _{\mathrm{max}}`$.

Our most important result is shown in Fig. 2, which shows the strong dependence of $`\eta _{\mathrm{max}}`$ on the trap geometry. For a fixed value of $`\omega _\rho `$, $`\eta _{\mathrm{max}}`$ can be optimized by increasing the axial trap frequency $`\omega _z`$, that is, by making the trap more “pancake” like. In contrast, for a fixed value of $`\omega _z`$, $`\eta _{\mathrm{max}}`$ decreases as $`\omega _\rho `$ is increased. It is important to realize that both quantities $`\omega _C`$ and $`\omega _J`$ are effected by changing the aspect ratio of the trap. By making the trap more “pancake” like, one is increasing the area of contact between the two wells, which increases $`\omega _J`$. On the other hand, compressing the condensate in the $`z`$ dimension also leads to a higher density, causing $`\omega _C`$ to increase as well. However, the ratio of these two quantities in (32) increases for increasing $`\omega _z`$. In contrast, by making the trap more “cigar” like, one is decreasing the area of contact between wells, resulting in a decrease in $`\omega _J`$ as $`\omega _\rho `$ increases. Furthermore, in this case one is compressing the condensate in two dimensions rather than one, which results in a more dramatic increase in the density and a corresponding increase in $`\omega _C`$, so that the ratio $`\omega _J/\omega _C`$ decreases with increasing $`\omega _\rho `$.

In Fig. 3 we vary the condensate population $`N_c`$ from $`10^2`$ to $`10^5`$ atoms, and show $`\eta _{\mathrm{max}}`$ for both 23Na and 87Rb. We see in Figs. 2 - 4 that a much larger amplitude of oscillation $`\eta _{\mathrm{max}}`$ can be achieved using 23Na compared to 87Rb. This behavior is expected based on the rough scaling of $`\eta _{\mathrm{max}}`$ given in (33). We can understand why 23Na is favorable over 87Rb, since the mass $`m`$ of 23Na is a about a factor of four smaller than that of 87Rb, and the scattering length $`a`$ of 23Na is almost half that of 87Rb. As $`N_c`$ increases in Fig. 3, the effect of the mean field interaction becomes more pronounced causing $`\omega _C`$ to increase, while this has only a small effect on $`\omega _J`$. We also show the dependence of $`\omega _{\mathrm{JP}}`$ on the condensate number $`N_c`$ in the inset.

We also show the dependence of $`\eta _{\mathrm{max}}`$ on the Gaussian barrier parameters by varying the relative height $`\chi `$ (31) in Fig. 4 and the width $`\sigma `$ in Fig. 5. We note that varying these two parameters has very little effect on $`\omega _C`$, which remains roughly constant in these plots. The behavior of $`\eta _{\mathrm{max}}`$ is not surprising: as the barrier height and width are each increased, the coupling between wells is diminished so that $`\omega _J`$ decreases, resulting in a decrease in $`\eta _{\mathrm{max}}`$ as indicated by (32).

## IV Conclusions

In this paper we have focussed on an important property of the double-well Josephson system for a dilute atomic Bose-Einstein condensate: the maximum amplitude of oscillation $`\eta _{\mathrm{max}}`$ of the population between wells is suppressed due to the mean field interaction Milburn1997a ; Zapata1998a ; Raghavan1999b ; Salasnich1999a . For typical magnetic trap parameters, we have shown that $`\eta _{\mathrm{max}}`$ is very small, indicating that Josephson oscillations may be visible only for the very shallow trap used in the experiments described in Matthews1999b . We have shown that much larger values of $`\eta _{\mathrm{max}}`$ can be achieved by decreasing the aspect ratio $`\lambda =\omega _\rho /\omega _z`$, that is, by making the trap more “pancake” like, in particular, for small values of the radial frequency $`\omega _\rho `$. Our study was motivated by the fact that Josephson oscillations were clearly observed in the experiment reported in Ref. Anderson1998a , where a relatively small condensate was loaded into an array of wells with an extremely small aspect ratio $`\lambda 0.002`$ Anderson1998a ; AndersonPrivate . We have also shown that larger values of $`\eta _{\mathrm{max}}`$ can be attained by decreasing the condensate population $`N_c`$, and decreasing the barrier height $`\chi `$ and width $`\sigma `$. Furthermore, conditions are much more favorable for 23Na compared to 87Rb. It should be noted, however, that present traps for 23Na atoms possess the less optimal “cigar” geometry.

We have restricted our calculations to zero temperature. Our study of the optimal conditions for observing Josephson oscillations should be extended to finite temperatures, where dissipative effects associated with the thermal cloud of non-condensate atoms must be included Zapata1998a ; Ruostekoski1998a . Very few explicit calculations of damping rates of Josephson oscillations in a dilute Bose gas have been carried out, and further study of damping in this system is needed.

I would like to thank Prof. A. Griffin for suggesting this problem and for his useful comments. I thank B. Anderson for many insightful comments and E. Zaremba, B. Jackson, and R. Spekkens for useful discussions. This work was supported by a grant from NSERC.