text:
# Radio Profile Instability in the Millisecond Pulsar PSR J0437-4715 due to Spiky Emission

## 1 Introduction

The coherent radio emission mechanism of rotation powered pulsars is as yet an unsolved problem. The broad similarity in the coherent radio emissions from normal and millisecond pulsars (McConnell et al. , Gil & Krawczyk , Jenet et al. , Vivekanad et al. , Kramer et al. \[1999b\], Vivekanand ) is intriguing; that a similar mechanism operates over three to four orders of magnitude difference in their periods and period derivatives is probably an important clue to the coherent radio emission mechanism of rotation powered pulsars.

Of equal importance are the differences in the radio emission from these two classes of pulsars. One of them is the slow variation in the IP of millisecond pulsars (Camilo , Backer & Sallmen , Vivekanad et al. , Kramer et al. \[1999a\]), over time scales ranging from several hundred periods to hours and even days, which is quite unlike the fast mode changing phenomenon observed in some normal pulsars (Backer , Helfand et al. , Manchester & Taylor , Rathnasree & Rankin ). Neither phenomenon has a satisfactory explanation so far.

Here it is shown that the slow IP variations seen in one of the components of PSR J0437$``$4715 (Vivekanand et al. ) are correlated with variations of the spiky emission in that component (Vivekanand ).

## 2 Profile Instability in PSR J0437$``$4715

Figure 1 shows IP of PSR J0437$``$4715 at five different epochs, observed at 326.5 MHz using Ooty Radio Telescope (ORT); only a single linear polarization was available. From top to bottom the five panels correspond to (1) 9 files observed from UT 15:27:59 to UT 15:42:08 on 7 Mar 1995, (2) 2 files observed at UT 15:43:57 and UT 15:45:49 on 7 Mar 1995, (3) 12 files observed from UT 16:04:17 to UT 16:23:22 on 7 Mar 1995, (4) 13 files observed from UT 15:39:25 to UT 16:04:21 on 9 Mar 1995, and (5) 3 files observed at UT 15:51:21, UT 15:53:02 and UT 15:54:40 on 14 Mar 1995. Each file consists of 12 125 useful periods of data. The instrument, the method of observation and data analysis are described in Ables et al. , Vivekanand et al. , Vivekanand , and Vivekanand .

In fig. 1 the first panel from top corresponds to the profile labelled type B in fig. 8 of section 5 of Vivekanand et al. . The second panel contains the two files that show a transition from type B to type A, and the third and fourth panels together correspond to the profile labelled type A. The main difference between the profiles of fig. 1 here and fig. 8 of Vivekanand et al.  is the more accurate period, differing by a few nanoseconds, that is used for folding the data (Vivekanand ).

On 7 Mar 1995, the peak flux of the third component of the IP of PSR J0437$``$4715 decreased from 3.3 to 2.6 over a period of $``$ 36 minutes (panels 1 to 3 of fig. 1). During this time the peak flux of its central component increased from 5.2 to 5.6, and that of its first component increased marginally from 2.1 to 2.3. These correspond to fractional flux variations of 21.2%, 7.7% and 9.5% respectively, which are much larger than the $``$ 1% to 2% maximum allowed by purely random noise. Two days later PSR J0437$``$4715 had an IP similar to that observed earlier (panel 4). Five days later its first component became stronger than its third component (panel 5), exactly contrary to its behavior a week earlier (panel 1).

Vivekanand et al.  argued that these profile variations are intrinsic to PSR J0437$``$4715, and are not an artifact due to the single polarization observations of ORT coupled with Faraday rotation in the variable ionosphere. In the next section it will be shown that the IP variations of component 3 are highly correlated with the occurrence of the spiky emission in that component.

## 3 Probability of Occurrence of Spikes

On account of limited signal-to-noise ratio of the data (Vivekanand ), only the 500 strongest spikes in each data file were analyzed; ideally one would have liked to completely separate the spiky from the non-spiky emission for this exercise.

First, the 500 strongest spikes in each data file were identified and their positions (phases) were obtained as described in Vivekanand . Next, these positions were plotted against the numbers of the periods in which they occurred, for those spikes which satisfied following two selection criterion: (1) their peak signal to noise ratio should be greater than 12.5, and (2) the error on their position should be less than 20 micro seconds ($`\mu `$s). However the final result is not very sensitive to these two selection criterion.

Figure 2 plots the positions of 375 of the 500 strongest spikes in the data file labelled 50661527, as a function of the period number. This belongs to the profile category labelled B by Vivekanand et al. ; the rms errors on the positions are less than the width of the dots. There are 52, 314 and 9 spikes in the three main components of the IP, respectively. If the two selection criterion are dropped, one obtains 66, 421 and 13 spikes in the corresponding components of the IP, which are entirely consistent with earlier numbers. In file labelled 50661538 belonging to the same category, that was observed at UT 15:38:41 on 7 Mar 1995, PSR J0437$``$4715 was at one of its lowest average fluxes. The number of spikes in the three main components of the IP, with and without the selection criterion, are 39, 301 and 1, and 56, 442 and 2, respectively, which are again fractionally similar. Therefore the two selection criterion do not significantly alter the results of this article.

Figure 3 plots the positions of the 500 strongest spikes in the data file labelled 50661612, which belongs to the category labelled A by Vivekanand et al. . Here the number of spikes in the three main components of the IP are 46, 453 and 1, respectively; no spikes have been selected out due to the high average pulsar flux in this data file. Clearly the fraction of spikes in the third component of the IP has decreased significantly compared to category B. This is verified in a visual check of as many periods as is reasonably possible. Even in file labelled 50661538, several periods had spikes in the third component of the IP, that were not part of the 500 strongest spikes.

To test this more rigorously, the above numbers were totalled for all files in each category of the IP; the results are shown in Table 1. The 9 files of panel 1 in fig. 1 yielded totally N0 = 3598 spikes; the rest of the 4500 spikes were ignored due to the two criterion mentioned above. Out of these, N1 = 408 spikes occurred in the first component of the IP, N2 = 3155 spikes occurred in the second component, and the N3 = 35 spikes occurred in the third component. This was for the type B profile. N0 is 12 405 for the 25 files of type A profile (panels 3 and 4 in fig. 1). From this one would have expected N1 = 408 / 3598 $`\times `$ 12 405 $``$ 1407 $`\pm `$ 38, whereas one obtained the smaller number 1255. Similarly one expected N2 = 3155 / 3598 $`\times `$ 12 405 $``$ 10 878 $`\pm `$ 104 in the second component, but obtained the larger number 11 144. But the most interesting number is the expected N3 = 35 / 3598 $`\times `$ 12 405 $``$ 121 $`\pm `$ 11, which is significantly larger than the observed value of 6; and this occurs in spite of the average pulsar flux increasing by about 14% (Table 1). Clearly the decrease in rate of occurrence of spikes in the third component of the IP of PSR J0437$``$4715 is correlated with the decrease in average flux of that component.

This trend continues in the last row of table 1 (panel 5 of fig. 1); but caution is required here. Since one selects the 500 strongest spikes in each data file, the results should not depend upon the average flux of PSR J0437$``$4715 to the zeroth order. But for very low average flux of PSR J0437$``$4715 one expects to see no spike in the third component of the IP, since the probability of occurrence there is very low anyway (Vivekanand ). This is probably what is happening in the last row of table 1; the average number of spikes per file is significantly lower here.

The expected values of N1, N2 and N3 in the last panel of fig. 1, based on the type B profile, are 39 $`\pm `$ 6, 302 $`\pm `$ 17 and 3 $`\pm `$ 2; these are consistent with the observed values. One would have liked to see a larger value for N1 here, to be consistent with the increased flux of the first component of the IP; but one is probably dealing with small numbers due to limited data, as well as a significantly weaker pulsar flux in those files.

Table 1 has also been derived after removing the two selection criterion, and the results remain similar.

The N3 values in the 9 files of row 1 of table 1 range from 9 to 1. The first three files observed chronologically, labelled 50661527, 50661529 and 50661531, have N3 = 9, 5 and 7, respectively. Now PSR J0437$``$4715 was strong in these files and weakened systematically with time. Not surprisingly, N3 vales in the next six files are 3, 2, 4, 1, 2 and 2; these files are labelled 50661533 through 50661542. It is clear that there was no hope of seeing any spikes in the third component of the IP in the next two files, viz. 50661543 and 50661545 in panel 2 of fig. 1, even if they existed, due to the low average flux of 2.54 $`\pm `$ 0.04 in these two data files; which is why the single pulse mode of observation was discontinued after that time. However, after about 18 minutes, the average flux of PSR J0437$``$4715 became high again; which is why the single pulse mode of observation was re-started once again. But N3 for the next 12 data files, belonging to panel 3 of fig. 1, was a mere 5; only five data files had N3 = 1, and the rest had N3 = 0. More convincingly, N3 was a mere 1 for the 13 data files belonging to panel 4 of fig. 1; and these contain some of the highest average flux data on PSR J0437$``$4715.

It is therefore concluded that the decrease in N3 between the first two rows of table 1 is a feature intrinsic to PSR J0437$``$4715. That this decrease correlates with the decrease in the flux of the third component of the IP proves the basic hypothesis of this article. The clinching evidence would have been if one had actually seen the decrease in the number of spikes as time progresses (in figures similar to fig. 2 and fig. 3), in the two data files labelled 50661543 and 50661545 of panel 2 of fig. 1, which are believed to be in transition between categories B and A. Unfortunately, this is not possible due to the low average flux of PSR J0437$``$4715 in these two data files.

## 4 Summary and Discussion

Vivekanand et al.  showed that the IP of the milli second pulsar PSR J0437$``$4715 changes significantly in shape over time scales of minutes and hours and even days. This is unlike the much faster mode changing known in normal pulsars. They argued that this behavior was intrinsic to PSR J0437$``$4715, and not an instrumental effect. Here it is shown that profile variations in the third component of the IP of PSR J0437$``$4715 are correlated with the rate of occurrence of the spiky emission in that component. Each period that contained any of the 35 + 6 = 41 spikes, in the third component of the IP in Table 1, was visually checked for confirmation.

Before proceeding further, two comments are in order:

1. The first concerns the availability of only a single polarization at ORT. Navarro et al.  show that the percentage of linear polarization of PSR J0437$``$4715 at 438 MHz is highest at the peaks of components 2 and 3 in fig. 1 (also their fig. 1); component 1 has very low linear polarization. From their fig. 1 the position angle of linear polarization differs by less than $`40^{}`$ at these two peaks. Therefore, it is possible that, in principle, differential Faraday rotation in a variable ionosphere at ORT can reduce the effective flux density of spikes in component 3 by $`1\mathrm{cos}^2(40^{})`$ 41% or less, in spite of the arguments of Vivekanand et al. . This issue can be resolved satisfactorily only when full polarization observations of PSR J0437$``$4715 are available at 326.5 MHz. However, one would have expected an as frequent reduction of number of spikes in component 2 also, with respect to the number of spikes in component 3; this is not observed, although the statistical basis of such a conclusion is not strong. Furthermore, it is quite possible that the mean position of spikes, in components 2 and 3 of fig. 1, is not coincident exactly with the corresponding position of the peak of linear polarization. In such an eventuality, the effect discussed above would be less serious, since the linear polarization profiles in components 2 and 3 are fairly sharply peaked. Finally, similar profile variations have been noted by Bell et al.  in PSR J0437$``$4715 at 430 MHz in full polarization observations. This author feels they might have wrongly attributed the observed profile variations to “relative gain variation in the two polarization channels and the high degree of polarization of the emission observed across the profile”.
2. Vivekanand  has shown that the position of occurrence of spikes is uncorrelated with their height; no simple relation exists between these positions and the integrated profile (figures 5 and 6 of that paper). Therefore this author can not think of any obvious selection effect involving these parameters that can produce the result of this paper.

One would have liked to find a quantitative correlation between the rate of occurrence of spikes (Table 1) and the average flux in the other two components of the IP (fig. 1) also. Unfortunately this is complicated by two factors. Firstly, the average pulsar flux in a data file could be so low that the observed rate of occurrence of spikes in any component may be biased to a lower value. Secondly, the average flux of the spiky emission can also change in principle, from file to file and from component to component. To illustrate the problem, the fraction of N3 spikes in Table 1 decreases down the rows, consistent with decrease in flux of the third component of the IP in fig. 1. However, the fraction of N1 spikes in Table 1 also decreases, although marginally; this is inconsistent with the increase in flux of the first component of the IP; this can easily arise due to the reasons mentioned above, as also due to a small change in the occurrence rate of spikes in the second component of the IP, which has about 90% of the spikes anyway.

PSR J0437$``$4715 has also been studied at $``$ 1400 MHz by several workers, but this kind of behavior has not been reported so far. One possible reason is the low signal to noise ratio of their data. The other possible reason could be that PSR J0437$``$4715 does not show such behavior at other frequencies. There is a precedence for this – PSR B1821-24 shows radio profile instability at 1400 MHz but not at 800 MHz (Backer & Sallmen ).

It is likely that the IP variations of the first component of PSR J0437$``$4715, as well as any other component, are also due to the rate of spiky emission in that component, although that has not been proved rigorously here.

It is further speculated here that the exact shape of the three principal components of the IP of PSR J0437$``$4715 is a consequence of the spiky emission, and not its cause. This is supported by fig. 4, where about 16% of the periods containing the 2000 strongest spikes in the entire data have an IP that has intense peaks in the three components; small variations in the spiky emission in any one of these three components is quite likely to change the shape of the overall IP of PSR J0437$``$4715.

This further enhances the doubt, expressed in the last section of Vivekanand , about the interpretation of the individual components of the IP of PSR J0437$``$4715 as individual beams or regions of emission. It is worth exploring the suggestion in that paper regarding the independent method of resolving the IP into independent beams or regions of emission on the pulsar.


format: text

text:
# Correlated electrons and generalized statistics

## I Introduction

In this work, we attempt to relate generalized statistical theories to correlated electron or quasiparticle systems which, since 1980’s, have been subjected to intensive theoretical and experimental studies due to their important roles in optical, mechanical, electric and magnetic  properties of metal compounds with partially filled d- or f-electrons. In general, with complex correlations, the conventional Boltzmann-Gibbs statistics (BGS) for free fermions is no longer valid. Most of the theoretical study up to now was numerical simulations based on approximate models including essential physics of the systems of interest. Although these theoretical works are very important and fruitful, the need of general understanding and theoretical formulation of the distribution laws begins to be felt. Due to the complexity of the correlations, we can conjecture that generalization of Fermi-Dirac (FD) theory may be useful. In this sense, some efforts have been made to understand localized electron behavior on the basis of fractal geometry (see  and references there-in), a generalization of normal space-time.

In what follows, I present a work in this direction based, not on the fractal or chaotic assumption, but on generalized statistical theories proposed for complex systems showing correlated phenomena. The theoretical results will be compared to numerical and experimental ones for correlated electrons. It is expected that this attempt can help to better understand different generalized theories and identify valid theories for this special circumstance.

## II Some important generalizations of BGS and correlated electrons

Indeed, in the last decade, BGS theory has experienced a turbulent period with the rapid development of some anomalous theories to treat complex systems. These theories are anomalous in the sense that they may perturb our old conceptions concerning, e.g., information, energy, quantum states, additivity etc. The main characters of these theories can be resumed as follows :

1) Consideration of nonadditivity in entropy, energy, quantum occupation number etc. depending on empirical parameters, e.g. $`0<q<\mathrm{}`$ in Tsallis nonextensive statistical mechanics (NSM) and in its incomplete statistics (IS) version as well as in quantum group theory (QGT), $`\kappa 0`$ in $`\kappa `$-statistics (KS).

2) Bosonization of fermions or fermionization of bosons allowing intermediate occupation number related to empirical parameter, e.g. $`0\alpha 1`$ in fractional exclusion statistics(FES) and $`1\eta 1`$ for the FES in the frame of KS.

3) The empirical parameters take some particular values when the theories recover BGS ($`q=1`$, $`\alpha =0,1`$, $`\kappa =0`$ and $`\eta =1,0,1`$).

The philosophy of these generalizations, explicit or not, is to introduce empirical parameters to “absorb” the effects of complex correlations. So interacting systems can be mathematically treated as noninteracting or conventional ones. For example, the total energy of an interacting “free particle” can be written as $`p^2/2m`$ where $`p`$ is the momentum and $`m`$ the mass of the particle. The extra interaction energy is in this way “absorbed” in the empirical parameter when they are different from the values corresponding to conventional cases.

About the generalization of FD statistics, Figure 1 shows the different fermion distributions at $`T`$=100 K of the generalized statistical mechanics mentioned above. We see that the discontinuity in occupation number $`n`$ at Fermi energy $`e_f`$ depends almost only on temperature for all these distributions. The sharp $`n`$-drop hardly changes for whatever value of the empirical parameters. In addition, we note that the approximate distribution (AD) $`n=\frac{1}{[1+(q1)\beta (ee_f)]^{1/(q1)}+1}`$ ($`0<q1`$ in order that the inverse temperature $`\beta `$ can go to $`\mathrm{}`$) of NSM obtained with factorization approximation by neglecting correlation energy between the particles is only slightly different from FD one. So its Fermi energy is almost identical to that of FD. The situation of KS is similar. At low temperature, the $`\kappa `$-distribution of standard fermion ($`\eta =1`$) $`n=\frac{1}{[\sqrt{1+\kappa ^2\beta ^2(ee_f)^2}+\kappa \beta (ee_f)]^{1/\kappa }+1}`$ is hardly different from FD for $`\kappa `$ values which give finite internal energy (e.g. $`0\kappa <1/2`$ for two dimensional fermion gas). For FES fermion distribution $`(1/n\alpha )^\alpha (1/n\alpha +1)^{1\alpha }=e^{(ee_f)/kT}`$ (or equivalently $`n=\frac{1}{e^{q\beta (ee_f)}1}\frac{(1+\alpha )/\alpha }{e^{q\beta (ee_f)(1+\alpha )/\alpha }1}`$), the $`n`$-discontinuity depends only on temperature just as for FD. The influence of $`\alpha `$ smaller than unity (for FD) is to decrease $`e_f`$ ($`\alpha `$ times) and to increase the maximal occupation number ($`1/\alpha `$). As for the nonextensive fermion distribution of IS $`n=\frac{1}{[1+(q1)\beta (ee_f)]^{q/(q1)}+1}`$, though $`n`$ drop at $`e_f`$ is always very sharp, $`e_f`$ increase considerably as $`q`$ decreases from unity. And consequently, there is an important wide $`n`$ decrease with decreasing $`q`$. It was shown that, for $`q0`$, $`e_f`$ can increase up to two times the $`e_f`$ of FD with $`n=1/2`$ ($`n=0`$) for all energy below (above) $`e_f`$. Another common character of NSM and IS distributions is that, due to the energy cutoff with $`q<1`$, there are few electrons above $`e_f`$ at low temperature (large $`\beta `$).

Now let see the distributions of correlated electrons. From some results of experiments and numerical simulations based on low dimension Kondo lattice models (KLM), we notice two important effects of correlations. First, the wide decrease of $`n`$ and the sharp $`n`$ drop at $`e_f`$ much larger than that without correlation have indeed been observed for very strongly correlated electrons (with the coupling parameter $`J4`$) (see Figure 3 and 4). So it is possible to describe correlated electron systems in strong coupling regime by IS fermion distribution with $`q`$ very smaller than unity (close to zero). On the other hand, another effect of correlation, in the weak coupling regime, is the flattening of $`n`$ drop at $`e_f`$. That is that the correlation, even at low temperature, drives electrons above $`e_f`$ so that the $`n`$ discontinuity becomes less and less sharp as the correlation increases. Curiously, this flattening of $`n`$ discontinuity at $`e_f`$, though confirmed by experimental results, is completely absent in the fermion distributions given by all the statistics, extensive (FES) or nonextensive (NSM, IS and KS), mentioned above.

Fermi surface plays an essential role in the physical properties of metals and metal compounds. So it is of great interest for generalized statistical theories to take into account the essential physics relevant to Fermi energy. In what follows, we will show that it is possible to cover the observed weak correlation effects upon $`e_f`$ within an additive generalization of BGS based on the IS principles. The reader can find a complete review of IS theory in references .

## III Additive incomplete statistics

The concept of IS is inspired by incomplete probability distribution (i.e. $`\text{Tr}\rho =Q1`$) as well as the theoretical difficulties one encounters with $`\text{Tr}\rho ^q1`$ in NSM. The basic assumption of IS is that our knowledge about physical systems is in general incomplete due to unknown interactions or their effects which can not be studied explicitly. In this case, there are always missing informations about the physical systems. So probability distribution is in general inexact and incomplete and can not sum to one. One should write $`\text{Tr}(\rho /Q)=\text{Tr}F(\rho )=1`$ where $`F`$ is certain function to be determined. In the case of complete or approximately complete distribution (such as BGS), $`F`$ is identity function. In my previous paper, in order to keep Tsallis framework of NSM, I proposed $`F(\rho )=\rho ^q`$ so that

$$\text{Tr}\rho ^q=1,$$
(1)
where $`q`$ is Tsallis entropy index. Since $`\rho <1`$, we have to set $`q[0,\mathrm{}]`$. When $`q=1`$, $`\rho `$ becomes complete, which implies a complete or quasi-complete knowledge about the system of interest. The nonextensive IS formalism on the basis of Tsallis entropy yields Tsallis $`q`$-exponential distributions and the fermion distribution discussed above. An successful application of Eq.(1) to the study of some power laws on the basis of Rényi entropy can be found in reference . One can even find a plausible justification of Eq.(1) in a work of Tsallis on the probability distribution on some simple fractal supports like Koch or Cantor sets. Indeed, the phase space of the systems including complex interactions is in general chaotic or fractal structure. In these structures, exact probability distribution is impossible due to the fact that, first, the structures are not differentiable and integrable, second, as mentioned above, some state points in the structures are not accessible. So the treatment of the incomplete distributions, not in the fractal space of $`d_f`$ dimension, but in a differentiable and integrable $`d`$-dimension space, is inevitable. This consideration yields in a natural way the incomplete normalization Eq.(1).

Now we suppose that the system of interest has additive information and physical quantities (e.g. $`H=_iH_i`$ where $`H`$ is the hamiltonian of the compound system and $`H_i`$ that of the $`i^{th}`$ subsystem). This case may happen if there are only weak and short range interactions, so that the Hartley formula $`\mathrm{ln}(1/\rho )`$ for information measure should apply. The entropy of the system can be defined à la Shannon as follows :

$$S=k\text{Tr}\rho ^q\mathrm{ln}(1/\rho ).$$
(2)
$`S`$ obviously becomes Gibbs-Shannon entropy ($`S_{GS}`$) when $`q=1`$, which identifies $`k`$ to Boltzmann constant. It is straightforward to see that $`S`$ verifies all the properties of $`S_{GS}`$. For microcanonical ensemble ($`\rho ^q=\frac{1}{w}`$), we have $`S=\frac{k}{q}lnw`$ which decreases with increasing $`q`$ value. In general, $`\mathrm{\Delta }S=SS_{GS}<0`$ (or $`>0`$) if $`q>1`$ (or $`q<0`$).

For canonical ensemble, we postulate

$$U=\text{Tr}\rho ^qH$$
(3)
where $`H`$ is the Hamiltonian and $`U`$ the internal energy of the system. The maximum entropy of $`S`$ with Eq.(1) and (3) as constraints leads to :

$$\rho =\frac{1}{Z}e^{\beta H}$$
(4)
with $`Z=(\text{Tr}e^{q\beta H})^{1/q}`$. The Lagrange parameter $`\beta `$ can be determined by $`\frac{S_q}{U_q}=k\beta =\frac{1}{T}`$ where $`T`$ is the absolute temperature. It is easy to see that all thermodynamic relations in this theory are identical to those in BGS.

We indicate in passing that this generalization of BGS is obviously the consequence of the “power-normalization” Eq.(1). Other useful normalizations, if any, of incomplete distributions may lead to different generalization of BGS. For example, if we use Hartley formula as information measure and define expectation value for additive entropy and energy with $`F(p_i)`$ satisfying $`\frac{\mathrm{ln}F(x)}{x}=\frac{1/x}{\alpha \mathrm{ln}x(\mathrm{ln}Zx)^{1/\gamma }}`$ where $`\gamma `$ is an empirical parameter, maximum entropy will lead to the famous stretched exponential distribution $`p_i=\frac{1}{Z}e^{(\beta E_i)^\gamma }`$ where $`E_i`$ is positive energy of the system at state $`i`$, $`\alpha `$ and $`\beta `$ the Lagrange multipliers related respectively to normalization and energy constraint $`U=_{i=1}^wF(p_i)E_i`$. When $`\gamma =1`$, we can recover $`\frac{\mathrm{ln}F(x)}{x}=\frac{1}{x}`$ as in BGS.

## IV Additive IS Fermion distribution

For grand canonical ensemble, the same machinery as above leads to

$$\rho =\frac{1}{Z}e^{\beta (H\mu N)}$$
(5)
with $`Z=[\text{Tr}e^{q\beta (H\mu N)}]^{1/q}`$. Supposing $`U=_jn_je_j`$ and $`N=_jn_j`$ in accordance with the additivity assumption, the average occupation number $`n`$ of one-particle state of energy $`e`$ can be straightforwardly calculated :

$$n=\frac{1}{e^{q(ee_f)/kT}+1}$$
(6)
which recovers FD distribution for $`q=1`$. The Fermi energy can be calculated in the standard way. For 2-D systems, e.g., we have $`e_f=e_{f_0}+k_qT\mathrm{ln}(1e^{e_f/k_qT})`$. In Figure 2 is plotted the $`T`$-dependence of $`e_f`$ for different $`q`$ values for 1-D fermion systems. We see that $`e_f`$ considerably changes with decreasing $`q`$. From Eq.(6), we see that $`q`$-dependence of $`e_f`$ at given $`T`$ will be similar to the $`T`$-dependence. The effect of $`q`$ on other physical quantities (specific heat, electrical conductivity, susceptibility, effective mass etc.) is discussed in .

The $`n`$ distribution around the Fermi momentum $`k_f`$ can be estimated from Eq.(6) using $`e=\mathrm{}^2k^2/2m`$. We get :

$$n=\frac{1}{2}[1+\frac{q}{kT}\frac{\mathrm{}^2k_f^2}{m}(kk_f)]$$
(7)
which is similar to the analog given by a Monte-Carlo calculation based on a tight-binding Kondo lattice ($`KL`$). A comparison with  leads to $`Z=\frac{q\mathrm{}k_f}{2}`$, where $`Z`$ is the quasiparticle weight of photoemission spectrum. Note that $`Z`$ is $`k_f`$-dependent here. Eq.(7) tells us that decreasing $`q`$ yields a flattening of $`n`$ drop at $`e_f`$.

In Figure 3 and 4, we compare the momentum distribution given by Eq.(6) for $`T`$=50 K to some numerical results of the momentum distribution of 1-D correlated electrons. We see that, for about $`J<1`$ and $`q>0.003`$, KLM calculations can be well reproduced by Eq.(6). But for smaller $`q`$ or stronger coupling with $`J>1`$, a long tail in the KLM distributions begins to develop at high energy and can not be reproduced by Eq.(6). At the same time, a new Fermi surface at $`k_f+\pi /2`$ starts to appear and a sharp $`n`$ drop (energy cutoff) takes place at the new Fermi momentum. All these strong correlation effects turn out to be completely absent within the present IS Fermion distribution. It is clear that this discrepancy marks the limit of this formalism and perhaps implies that nonextensivity is no more negligible in very strong-coupling case. It is worth noticing that the strong correlation effects, absent in additive IS fermion distribution, are just what we noted in nonextensive IS fermion distribution plotted in Figure 1. This suggests that a combination of these two IS generalizations of BGS, representing respectively two different aspect of correlation, may be an interesting approach leading to more generally valid statistical mechanics suitable for correlated particle systems.

## V conclusion

We have compared several important generalizations of FD distribution to numerical and experimental results for correlated electron systems. It is found that the generalized fermion distributions based on incomplete information can be useful for describing these systems. It is shown that the extensive IS fermion distribution gives a good description of correlated electrons in weak coupling regime. On the other hand, it fails to describe strongly correlated conduction electrons and localized $`f`$-electrons which are possible to be described by the nonextensive IS fermion distribution. The latter, though powerless for the description of weakly correlated electrons, shows similar behavior as the heavy electrons in strong coupling regime : strong increase of $`e_f`$ and $`n`$-cutoff accompanied by a general decrease in $`n`$ at all energy up to $`e_f`$. In a current work, we are trying to combine these two partially valid generalizations. Further results will be presented in a future paper.

Figure caption :


text:
# The nucleation, growth and stabilization of vortex lattices

## Abstract

We give a simple unified theory of vortex nucleation and vortex lattice formation which is valid from the initiation process up to the final stabilization of the lattice. We treat the growth of vortex lattices from a rotating thermal cloud, and their production using a rotating trap. We find results consistent with previous work on the critical velocity or critical angular velocity for vortex formation, and predict the initial number of vortices expected before their self assembly into a lattice. We show that the thermal cloud plays a crucial role in the process of vortex lattice nucleation.

The opportunity that Bose Einstein condensates present for precisely controlled experimental and theoretical study of superfluids has revived an intense interest in the properties and behavior of quantized vortices. Many different theoretical calculations have been made of stationary single vortex states and their excitations (Fetter and Svidinsky (2001) and references therein), and vortex lattices (Feder et al. (1999a) and references therein) in Bose condensates, but the recent experimental observations of vortex lattices Madison et al. (2000, 2001); Haljan et al. (2001); Raman et al. (2001); Abo-Shaeer et al. (2001, 2002); Hodby et al. (2001) have focussed attention on the mechanisms of vortex formation.

One experimental method for creating vortices is to stir a condensate with an anisotropic potential, and a number of theoretical analyses of this scenario have been made in terms of the Gross-Pitaevskii equation, e.g. Fetter (2001); Dalfovo and Stringari (2000); Anglin (2001); Muryshev and Fedichev . The view of these treatments is that the perturber causes a mixing of the condensate ground state and excited condensate states (with angular momentum values determined by the stirrer geometry). Typically, a perturbative calculation is used to obtain a critical rotational speed (or a critical linear speed at the Thomas-Fermi radius) of the stirrer at which the mixing becomes effective. The argument is then made that an instability will lead to growth of the vortex state. A more complete non-perturbative calculation for a localized rotating stirrer Caradoc-Davies et al. (1999) shows that in this case the mixing is predominantly between the ground state and an $`l=1`$ vortex state, and that as the relative amplitudes evolve by coherent “Rabi cycling”, the full condensate exhibits a vortex cycling from infinity to the central regions of the condensate. No coherent mixing mechanisms, however, can explain the formation of a vortex lattice, since an energy barrier exists between the superposition state and the vortex lattice state, and some additional mechanism is required to remove the energy liberated when the vortex lattice is formed.

An alternative method of producing vortices, in which a vapor of cold atoms is evaporatively cooled so as to preserve its angular momentum has been demonstrated by Haljan et al. Haljan et al. (2001). In this experiment, which involves no stirring potential, the mechanism of formation of the vortices must involve a modification of the theory of condensate growth (which now exists in a reasonably good quantitative form; Davis and Gardiner (2002) and references therein) to take account of the non-zero angular momentum of the vapor from which the condensate is formed. In this paper we will show that the mechanism thus demonstrated is the fundamental process in those experiments in which the stirring of the condensate apparently generates a vortex lattice. Weak stirring, in the absence of any thermal cloud, only produces “Rabi cycling” of angular momentum in and out of the condensate. However, if the stirring also produces a rotating thermal cloud, then, by the same process responsible for the growth of a condensate Gardiner et al. (1997), angular momentum is transferred irreversibly into the condensate, leading in equilibrium to a rotating vortex lattice. In general terms, this has already been noted in Zhuravlev et al. (2000), and in Gardiner et al. (2002) we derived a very simplified vortex growth equation—given below as Eq. (3)—embodying the essentials of this concept.

Other theoretical treatments of condensate dynamics at finite temperature have been applied to the case of vortex dynamics. Zhuravlev et al. Zhuravlev et al. (2000) studied the interaction of vortex lattices with a thermal cloud in the presence of a trap asymmetry, and from broad conservation principles they obtained equations for the time evolution of the angular velocity of the condensate and thermal cloud, and the transfer of energy between them. Williams et al. Williams et al. (2002) showed how the inclusion of interactions with a thermal vapor cloud would give rise to an instability of surface modes, as required to produce vortices. Our own treatment Gardiner et al. (2002) had a basis broadly equivalent to that of Williams et al., but showed how the interaction with the vapor cloud could produce an equation which should describe the full process of vortex nucleation, growth and stabilization. Using this equation, we will show in this paper how the growth of vortex lattices can be understood as a condensate growth process, and that the condition that the gain in this growth process be positive is the same as the criteria noted above for vortex nucleation. We also demonstrate explicitly how the vortices produced in the “Rabi cycling” caused by a rotating trap nucleate the growth process, which then proceeds to dominate the vortex lattice formation process.

The Gross-Pitaevskii equation  for $`\psi _R=\mathrm{exp}(i𝛀\mathbf{}𝐋t/\mathrm{})\psi `$, the condensate wavefunction in a frame of reference rotating with angular velocity $`𝛀`$, takes the form

$`i\mathrm{}{\displaystyle \frac{\psi _R}{t}}=\left({\displaystyle \frac{\mathrm{}^2}{2m}}^2+V_T^R+u\left|\psi _R\right|^2𝛀\mathbf{}𝐋\right)\psi _R,`$ (1)
where $`V_T^R(𝐱)`$ is the trap potential in the rotating frame. A vortex lattice results from the stationary solutions of Eq. (1) satisfying $`i\mathrm{}\psi _R/t=\mu _C\psi _R`$ Feder et al. (1999b, a), which are normally obtained by integrating the Gross-Pitaevskii equation along an imaginary time direction Feder et al. (1999b). Indeed, noting that the achievement of the vortex lattice stationary state requires dissipation in order to allow the system to settle into the energy minimum associated with the lattice, Tsubota et al. Tsubota et al. (2002) proposed an ad hoc modification of this method as a phenomenological model of the process of vortex lattice formation. They modified the left-hand side of Eq. (1) to $`(1+i\gamma )i\mathrm{}\psi _R/t`$, and compensated for the consequent loss of probability by continuously adjusting the norm of the wavefunction.

Origin of dissipation:  When the condensate is grown from a rotating vapor cloud the growth process itself is dissipative, and in the experiments of Madison et al. (2000, 2001); Haljan et al. (2001); Raman et al. (2001) dissipation can arise by transfer of atoms between the thermal cloud and the condensate. The vortex growth equation derived in Gardiner et al. (2002) uses quantum kinetic theory for the growth of the wavefunction for a condensate trapped in a potential $`V_T(𝐱,t)`$ rotating with angular velocity $`𝛀`$, from a thermal cloud held at temperature $`T`$ and chemical potential $`\mu `$ and rotating with angular velocity $`𝜶`$. A dissipative term arises from collisions between atoms in a thermal cloud of atoms trapped by the same potential as the condensate, in which one of the colliding atoms enters the condensate after the collision. The net rate at which atoms enter and leave the condensate as the result of such collisions is characterized by a coefficient

$`W^+(𝐱)`$ $`=`$ $`{\displaystyle \frac{u^2}{(2\pi )^5\mathrm{}^2}}{\displaystyle d^3𝐊_1d^3𝐊_2d^3𝐊_3\delta \left(\omega _{123}\right)}`$ (2a)
$`\times \delta (𝐊_1+𝐊_2𝐊_3)F_1F_2(1+F_3)`$
$``$ $`g\times 4m(akT)^2/\pi \mathrm{}^3,`$ (2b)
and in practice $`g3`$ fits most growth experiments. Here $`F(𝐱,𝐊)`$ is the Bose-Einstein distribution function for a thermal cloud in an effective potential in the frame rotating at the cloud’s angular velocity $`𝜶`$, i.e. $`V_\alpha (𝐱)\overline{V}_T(𝐱)m\left(𝜶\mathbf{\times }𝐱\right)^2/2`$, where $`\overline{V}_T(𝐱)`$ is the time averaged trap potential over the differential rotation between the trap and the cloud. The vortex growth equation takes the form in the frame rotating with the trap potential

$`i\mathrm{}{\displaystyle \frac{\psi _R}{t}}`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{}^2}{2m}}^2+V_T^R(𝐱)+u\left|\psi _R\right|^2𝛀\mathbf{}𝐋\right)\psi _R`$ (3)
$`+i\gamma \left(\mu +(𝜶𝛀)\mathbf{}𝐋i\mathrm{}{\displaystyle \frac{}{t}}\right)\psi _R,`$
where $`\gamma \mathrm{}W^+/kT`$. It is important to distinguish between the cloud’s angular velocity $`𝜶`$, and the trap’s angular velocity $`𝛀`$. If these are different the distribution of the thermal cloud is not truly thermodynamic equilibrium, and the form chosen represents an approximation expressing the experimental observation that it is very difficult to spin up a thermal cloud by means of a rotating potential. One therefore expects that $`\alpha <\mathrm{\Omega }`$. In this paper, however, we shall consider mainly the cases in which either the trap is cylindrically symmetric, so that we may set $`𝛀=0`$, or in which the cloud and the trap rotate at the same velocity, so $`𝜶=𝛀`$. Eq. (3) is essentially equivalent in terms of the physical assumptions required for its derivation to that of Williams et al. (2002), although its appearance is very different, and it is in practice easier to solve. In contrast, although Eq. (3) looks superficially very similar to that of Tsubota et al. (2002), it is in fact very different, as explained in Gardiner et al. (2002).

Simulation results for a rotating vapor cloud:  We have simulated Eq. (3) in two dimensions for many different scenarios, and in Fig. 1 we present a time sequence of results from a representative case in which the condensate is initially in the ground state of a rotationally symmetric harmonic trap of frequency $`\omega `$, the thermal cloud has $`\mu =12\mathrm{}\omega `$, and $`\alpha =0.65\omega `$. The initial condensate has unit norm and chemical potential $`\mu _C=12.7\mathrm{}\omega `$, although they can be any non-zero number. We simulate the non-stimulated collisions which start the process by adding to the intial wavefunction a uniform superposition of angular momentum states with $`l=1`$ to $`30`$, on a Gaussian radial profile centered at the Thomas-Fermi radius and with maximum amplitude $`2\times 10^7`$.

Fig. 1(a) is a density plot of the condensate near the end of the first stage of the evolution, in which an imperfect ring of $`19`$ vortices arrives from infinity to just outside the Thomas-Fermi radius. As this ring shrinks further (Fig. 1(b)), several vortices are shed, and a dominant ring of 16 vortices passes through the Thomas-Fermi radius into the interior of the condensate. These vortices then distribute themselves irregularly but quasi-uniformly over the condensate which expands and picks up angular velocity (Fig. 1(c)). Subsequently, over a long period, further vortices leave (and for larger values of $`\alpha `$ may enter) the dense region, until finally a regular lattice of 12 vortices rotating at angular velocity $`\alpha `$ in the lab frame remains (Fig. 1(d)).

The overall time scale of the process is illustrated in Fig. 2(a) which shows the total angular momentum increasing sharply around $`t40`$, when the 16 vortices enter the dense part of the condensate. For a given seeding, the timescale for this onset is determined by $`\gamma `$ and to good approximation scales as $`\gamma ^1.`$ The local chemical potential $`\mu _{\mathrm{loc}}(𝐱)`$ (defined as the absolute value of the right-hand side of Eq. (1) divided by $`\psi _R`$) provides a useful characterisation of the temporal development. In Fig. 2(b) , we plot the evolution of the spatial mean of $`\mu _{\mathrm{loc}}(𝐱)`$, which exhibits a rapid initial adjustment (to $`\mu `$), as the condensate roughly equilibrates with the thermal cloud, and a later dip at $`t48`$ during which the number of atoms in the condensate increases (see Fig. 2(a)). The spatial variance $`\sigma _\mu `$ of $`\mu _{\mathrm{loc}}(𝐱)`$ quantifies the deviation of the solution from equilibrium (where it is zero). It becomes large as the vortices enter the dense region of the condensate, and its subsequent slow decay indicates the long time needed for stabilization of the vortex lattice. It is not easy to characterize the latter time scale, since the local chemical potential can be almost constant everywhere, with possibly only a single vortex not in its final place. From a large number of simulations we have found that the number of vortices in the stable final lattice increases with $`\alpha `$ and $`\mu `$, and is generally different from the number seen in the initial ring. It also exhibits a very weak dependence on the exact form of the seeding. For the case of $`\mu =12\mathrm{}\omega `$ the critical value for the appearance of any vortices is found from the simulations to be $`\alpha =0.444\omega `$, at which speed an irregular ring of $`10`$ vortices appears initially, settling down eventually to three vortices in a lattice inside the Thomas-Fermi radius. Three dimensional simulations of Eq. (3) also confirm the general behavior seen in Fig. 1.

Analytic treatment:  An analytic understanding of the nucleation and growth process of the vortices can be obtained from Eq. (3) by considering the linearized Bogolyubov excitations above the initial state, which have angular momentum $`l`$ and energy eigenvalues $`ϵ_{n,l}`$ measured relative to $`\mu _C`$. From the final term of Eq. (3) we see that excitations for which $`ϵ_{n,l}<\mu \mu _C+\mathrm{}\alpha l`$ will experience a positive growth rate $`G\gamma (\mu \mu _C+\mathrm{}\alpha lϵ_{n,l})/\mathrm{}.`$ Initially this causes rapid growth of the condensate (the $`n=l=0`$ component, which has the lowest eigenvalue) until the condensate and thermal cloud approximately equilibrate $`\left(\mu \mu _C\right)`$. Subsequently, the gain of the other components is $`\gamma (\alpha lϵ_{n,l}/\mathrm{})`$, and thus the critical value of $`\alpha `$ for positive gain is given by $`ϵ_{0,l}=\mathrm{}\alpha l,`$ the same condition for criticality found in Williams et al. (2002). Above criticality, the dominant value $`l_v`$ is that which maximizes the gain, and is given by $`(\mathrm{}\alpha lϵ_{0,l})/l=0`$, and a crude estimate for $`l_v`$ can be obtained by minimizing $`\mathrm{}^2l^2/2mr^2\mathrm{}\alpha l`$ at the Thomas-Fermi surface, to give $`l_v\pi R_{\mathrm{TF}}^22m\alpha /h`$. This is the number of vortices which would result from filling a disk of radius $`R_{\mathrm{TF}}`$ with a vortex lattice Donnelly (1991). In the initial stages of vortex formation, when the process is still linear, the wavefunction of the condensate takes the form ( for simplicity including only the excitation with maximum gain, which has angular momentum $`\mathrm{}l_v`$ and energy $`\mathrm{}\omega _v\mathrm{}\alpha l_vϵ_{0,l_v}`$)

$`\psi _Re^{i\mu _Ct/\mathrm{}}\left\{\xi _0+e^{il_v\varphi +Gt}\left[ue^{i\omega _vt}+ve^{i\omega _vt}\right]\right\}.`$ (4)
Here $`\xi _0(r)`$ is the initial rotationally symmetric condensate wavefunction, $`\varphi `$ is the azimuthal angle, and $`u(r),v(r)`$ are obtained by solving the Bogolyubov-de Gennes equations. The essential behavior can be seen by neglecting $`v(r)`$, so that $`\psi `$ has $`l_v`$ zeroes—that is vortices—given by the interference of the two terms. These vortices all occur at the same radius, initially at infinity. Since the long distance behavior of $`u(r)`$ must be less rapid than that of $`\xi _0(r)`$, the ring will steadily shrink as the excitation grows. The time dependence of the coefficient means that the ring of vortices rotates at the angular frequency $`\omega _v`$ in the rotating frame. In a more detailed picture, more values of $`l`$ should be included, and the corresponding superposition will be an imperfectly circular ring of vortices, perhaps containing more than $`l_v`$ vortices.

This analytic treatment can be verified by decomposing the condensate into angular momentum components, and plotting the occupation values $`P_l`$ against time, as in Fig. 3(a).

We see that the $`l`$ components above threshold grow exponentially up to the time $`t40`$ (when the vortex ring passes through the Thomas-Fermi radius), and the gain is maximum for $`l_v=16`$. In Fig. 3(b) we compare the gain prediction $`\gamma (\alpha lϵ_{0,l}/\mathrm{})`$ (obtained by calculating the surface modes $`ϵ_{0,l}`$ for a condensate with $`\mu _C=12\mathrm{}\omega `$ ) against the gain rates measured from the simulation. Not only is our gain prediction very accurate, but we also find that the spatial particle density for each $`l`$ component very accurately matches the prediction from the corresponding Bogolyubov wavefunction.

Rotating trap, no thermal cloud:  Finally, we consider the connection between stirring and growth. As an example of a case of pure stirring we have simulated a rotating elliptical trap in the absence of a thermal cloud. The general behavior is well represented by the case $`\mu _C=12.7\mathrm{}\omega `$ and $`\mathrm{\Omega }=0.65\omega `$, for which we find that the condensate deforms into an elliptical shape rotating at the trap frequency, and a similarly deformed ring of about 28 vortices shrinks until it meets the dense part of the condensate at its narrowest radius. The resulting state is a distorted version of Fig. 1(b) Pen , and this subsequently undergoes an oscillatory behavior, but with no further penetration of the vortices into the condensate. Only the even angular momentum components of the condensate become significantly occupied and undergo periodic cycling (Fig. 4(a)). This can be well understood in terms of the Rabi cycling model Caradoc-Davies et al. (1999), noting that an elliptical potential can connect the initial ground state only to states of even $`l`$. Resonant mixing to the surface modes will occur when $`\mathrm{}\mathrm{\Omega }lϵ_{0,l}=0`$, and the lowest value of $`\mathrm{\Omega }`$ for which this occurs (when also $`(\mathrm{}\mathrm{\Omega }lϵ_{0,l})/l=0`$) gives the same condition for criticality as found in Refs. Dalfovo and Stringari (2000); Anglin (2001); Muryshev and Fedichev . This is also exactly the same condition as for the rotating cloud.

Rotating trap with a thermal cloud:  Now repeating this simulation, but with a thermal cloud (Eq. (3) with $`\alpha =\mathrm{\Omega }`$ and $`\gamma 0`$), we get dramatically different results. The same elliptical distortion of the condensate occurs, but the effect tends to smooth out with time. A distorted vortex ring appears, with an initial number of vortices as expected from the growth mechanism (here $`14`$), and then shrinks into the condensate and nucleates the vortex lattice. The angular momentum projections in Fig. 4(b) show the difference very clearly. The vortex lattice is the result of the growth process, which is essentially superposed on the effects of stirring, eventually masking them. This scenario leads to the conclusion that the principal effect of the rotating trap is to produce a rotating thermal cloud, which then transfers angular momentum into the condensate, eventually forming the vortex lattice.

Comparison with other work:  The main point of this paper is the analysis of vortex formation in terms of the growth in occupation of modes with non-zero angular momentum. This requires the presence of a rotating thermal cloud, but stirring the condensate can seed angular momentum components which may then grow from the thermal cloud by stimulated collisions. The critical angular velocity of the thermal cloud to allow growth is the same as the critical angular velocity required for stirring to resonantly mix surface modes with the condensate, and agrees with the latter values calculated by previous workers Dalfovo and Stringari (2000); Anglin (2001); Muryshev and Fedichev ; Williams et al. (2002). However, our treatment goes beyond the critical point, after which the behavior of the two processes becomes quite different, and includes the full dynamics of lattice stabilization. Our vortex growth equation (3) can be solved without linearization or using hydrodynamic methods, and its solutions can be interpreted in a straightforward quantum mechanical way. The dynamics of the thermal cloud are suppressed in this paper, but their inclusion is necessary for the full understanding of vortex lattice dynamics, including lattice decay, which can also be treated by our formalism.

###### Acknowledgements.

This work was supported by the Marsden Fund of New Zealand under contract PVT-902.


text:
# Machine Learning with Lexical Features: The Duluth Approach to Senseval-2

## 1 Introduction

The Duluth systems in Senseval-2 take a supervised learning approach to the Spanish and English lexical sample tasks. They learn decision trees and Naive Bayesian classifiers from sense-tagged training examples where the context in which an ambiguous word occurs is represented by lexical features. These include unigrams and bigrams that occur anywhere in the context, and co-occurrences within just a few words of the target word. These are the only types of features used. There are no syntactic features, nor is the structure or content of WordNet employed. As a result these systems are highly portable, and can serve as a foundation for systems that are tailored to particular languages and sense inventories.

The word sense disambiguation literature provides ample evidence that many different kinds of features contribute to the resolution of word meaning. These include part–of–speech, morphology, verb–object relationships, selectional restrictions, lexical features, etc. When used in combination it is often unclear to what degree each type of feature contributes to overall performance. It is also unclear to what extent adding new features allows for the disambiguation of previously unresolvable test instances. One of the long term objectives of our research is to determine which types of features are complementary and cover increasing numbers of test instances as they are added to a representation of context.

## 2 Experimental Methodology

The training and test data for the English and Spanish lexical sample tasks is split into separate training and test files per word. A supervised learning algorithm induces a classifier from the training examples for a word, which is then used to assign sense tags to the test instances for that word.

The context in which an ambiguous word occurs is represented by lexical features that are identified using the Bigram Statistics Package (BSP) version 0.4. This is free software that extracts unigrams and bigrams from text using a variety of statistical methods. Each unigram or bigram that is identified in the training data is treated as a binary feature that indicates whether or not it occurs in the context of the word being disambiguated. The free software package SenseTools (version 0.1) converts training and test data into a feature vector representation, based on the output from BSP. This becomes the input to the Weka suite of supervised learning algorithms. Weka induces classifiers from the training examples and applies the sense tags to the test instances.

The same software is used for the English and Spanish text. BSP and SenseTools are written in Perl and are freely available from www.d.umn.edu/~tpederse/code.html. Weka is written in Java and is freely available from www.cs.waikato.ac.nz/~ml.

## 3 System Descriptions

There were eight pairs of Duluth systems in the English and Spanish lexical sample tasks. The only language dependent components are the tokenizers and stop–lists. For both English and Spanish a stop–list is made up of all words that occur ten or more times in five randomly selected word training files of comparable size. All Duluth systems exclude the words in the stop–list from being features.

Each pair of systems is summarized below. All performance results are based on accuracy (correct/total) using fine-grained scoring. The name of the English system appears first, followed by the Spanish system.

Duluth1/Duluth6 create an ensemble of three Naive Bayesian classifiers, where each is based on a different set of features. The hope is that these different views of the training examples will result in classifiers that make complementary errors, and that their combined performance will be better than any of the individual classifiers.

Separate Naive Bayesian classifiers are learned from each representation of the training examples. Each classifier assigns probabilities to each of the possible senses of a test instance. These are summed and the sense with the largest value is used. This technique is used in many of our ensembles and will be referred to as a weighted vote.

The first feature set is made up of bigrams, i.e., consecutive two word sequences, that can occur anywhere in the context with the ambiguous word. To be selected as a feature, a bigram must occur two or more times in the training examples and have a log-likelihood ratio ($`G^2`$) value $``$ 6.635, which is associated with a p-value of .01.

The second feature set is based on unigrams, i.e., one word sequences, that occur five or more times in the training data.

The third feature set is made up of co-occurrence features that represent words that occur on the immediate left or right of the target word. In effect, these are bigrams that include the target word. They must also occur two or more times and have a log-likelihood ratio $``$ 2.706, which is associated with a p-value of .10.

These systems are inspired by \[Pedersen, 2000\], which presents an ensemble of eighty-one Naive Bayesian classifiers based on varying sized windows of context to the left and right of the target word that define co-occurrence features. However, the current systems only use a three member ensemble to capture the spirit of simplicity and portability that underlies the Duluth approach to Senseval-2.

English accuracy was 53%, Spanish was 58%.

Duluth2/Duluth7 learn an ensemble of decision trees via bagging. Ten samples are drawn, with replacement, from the training examples for a word. A decision tree is learned from each of these permutations of the training examples, and each of these trees becomes a member of the ensemble. A test instance is assigned a sense based on a weighted vote among the members of the ensemble. In general decision tree learning can be overly influenced by a small percentage of the training examples, so the goal of bagging is to smooth out this instability.

There is only one kind of feature used in these systems, bigrams that occur two or more times and have a log-likelihood ratio $``$ 6.635. This is one of the three feature sets used in the Duluth1/Duluth6 systems.

The set of bigrams that meet these criteria become candidate features for the J48 decision tree learning algorithm, which is the Weka implementation of the C4.5 algorithm. The decision tree learner first constructs a tree of features that characterizes the training data exactly, and then prunes features away to avoid over–fitting and allow it to generalize to the previously unseen test instances. Thus, a decision tree learner performs a second cycle of feature selection and is not likely to use all of the features that we identify prior to learning with BSP. The default C4.5 parameter settings are used for pruning.

These systems are an extension of \[Pedersen, 2001\], which learns a single decision tree where the representation of context is based on bigrams. This earlier work does not use bagging, and the top 100 bigrams according to the log-likelihood ratio are the candidate features.

English accuracy was 54%, Spanish was 60%.

Duluth3/Duluth8 rely on the same features as Duluth1/Duluth6, but learn an ensemble of three bagged decision trees instead of an ensemble of Naive Bayesian classifiers. There is a strong contrast between these techniques, since decision tree learners attempt to characterize the training examples and find relationships among the features, while a Naive Bayesian classifier is based on an assumption of conditional independence among the features.

The feature set used in these systems is from Duluth1/Duluth6 and consists of bigrams, unigrams and co-occurrences. A bagged decision tree is learned for each of the three kinds of features. The test instances are classified by each of the bagged decision trees, and a majority vote is taken among the members to assign senses to the test instances.

These are the most accurate of the Duluth systems for both English (57%) and Spanish (61%). These are within 7% of the most accurate overall approaches for English (64%) and Spanish (68%).

Duluth4/Duluth9 uses a Naive Bayesian classifier based on a bag of words representation of context, where each unigram that occurs in the training data is taken as a feature. This is a common benchmark in word sense disambiguation studies and text classification problems.

In the English training examples any word that occurs five or more times is used as a feature, and in the Spanish data any word that occurs two or more times is used. These features are used to estimate the parameters of a Naive Bayesian classifier. This will assign the most probable sense to a test instance, given the surrounding context.

Accuracy for English was 54%, and for Spanish 56%. This Naive Bayesian classifier was one of the three member classifiers in the ensemble approach of Duluth1/Duluth7, which was 1% less accurate for English and and 2% more accurate for Spanish.

Duluth5/Duluth10 add a co–occurrence feature to the Duluth2/Duluth7 systems. In every other respect they are identical. The co-occurrence feature was also used in Duluth1/Duluth6, and is essentially a bigram where one of the words is the ambiguous word. These must occur two or more times in the training examples and have a log-likelihood ratio $``$ 2.706 to be included as a feature. In addition to the co-occurrence feature the bigram feature from Duluth2/Duluth7 is used, where a bigram must occur two or more times and have a log-likelihood ratio $``$ 6.635.

Accuracy for English was 55%, and for Spanish 61%. This was a slight improvement over Duluth2 (54%) and Duluth7 (60%).

DuluthA/DuluthX build an ensemble of three different classifiers that are induced from the same representation of the training examples. A weighted vote is taken to assign senses to test instances. The three classifiers are a bagged J48 decision tree, a Naive Bayesian classifier, and the nearest neighbor classifier IB$`k`$, where the number of neighbors parameter $`k`$ is set to 1.

The context in which the ambiguous word occurs is represented by bigrams that may include zero, one, or two intervening words that are ignored. To be considered as features these bigrams must occur two or more times and have a log-likelihood ratio $``$ 10.827, i.e., a p-value of .001. The log-likelihood ratio threshold is set to 0 for the Spanish data due to the smaller volume of data.

English accuracy was 52%, Spanish was 58%.

DuluthB/DuluthY are identical to Duluth5/Duluth10, except that rather than learning an entire decision tree they stop the learning process once the root of the decision tree is selected. The resulting one node decision tree is called a decision stump. At worst a decision stump will reproduce the most common sense baseline, and may do better if the selected feature is particularly informative. In previous work we have observed that decision stumps can serve as a very aggressive lower bound on performance \[Pedersen, 2001\].

Decision stumps are the least accurate method for both English (DuluthB, 51%) and Spanish (DuluthY, 52%), but are more accurate than the most common sense baseline for English (48%) and Spanish (47%).

DuluthC/DuluthZ take a kitchen sink approach to ensemble creation, and combine the seven systems for English and Spanish into ensembles that assign senses to test instances by taking a weighted vote among the members.

Accuracy for English was 55%, and for Spanish 59%. This is less than the accuracy of some of the members systems, suggesting that the members of the ensemble are making redundant errors.

## 4 Discussion

There are several hypotheses that underly and motivate these systems.

### 4.1 Features Matter Most

This hypothesis is at the core of much of our recent work. It holds that variations in learning algorithms matter far less to disambiguation performance than do variations in the features used to represent the context in which an ambiguous word occurs. In other words, an informative feature set will result in accurate disambiguation when used with a wide range of learning algorithms, but there is no learning algorithm that can overcome the limitations of an uninformative or misleading set of features.

There are a number of demonstrations that can be made from the Duluth systems in support of this hypothesis, but perhaps the clearest is found in comparing the systems Duluth1/Duluth6 and Duluth3/Duluth8. The first pair learns three Naive Bayesian classifiers and the second learns three bagged decision trees. Both use the same feature set to represent the context in which ambiguous words occur. There is a 3% improvement in accuracy when using the decision trees. We believe this modest improvement when moving from a simple learning algorithm to a more complex one supports the hypothesis that the true dividends are to be found in improving the feature set.

### 4.2 50/25/25 Rule

We hypothesize a 50/25/25 rule for supervised approaches to word sense disambiguation. This loosely holds that given a classifier learned from a sample of sense–tagged training examples, about half of the test instances are easily disambiguated, a quarter are harder but still possible, and the remaining quarter are extremely difficult. This is a minor variant of the 80/20 rule of time management, which holds that 20% of effort accounts for 80% of results.

When the two highest ranking systems in the official English lexical sample results are compared there are 2180 test instances (50%) that both disambiguate correctly using fine-grained scoring. There are an additional 1183 instances (28%) where one of the two systems are correct, and 965 instances (22%) that neither system can resolve. If these two systems were optimally combined, their accuracy would be 78%. If the third-place system is also considered, there are 1939 instances (44.8%) that all three systems can disambiguate, and 816 (19%) that none could resolve.

For all the Duluth systems for English, there are 1705 instances (39%) that all eight systems got correct. There are 1299 instances (30%) that none can resolve. The accuracy of an optimally combined system would be 70%. The most accurate individual system is Duluth3 with 57% accuracy.

For the Spanish Duluth systems, there are 856 instances (38%) that all eight systems got correct. There are 478 instances (21%) that none of the systems got correct. This results in an optimally combined result of 79%. The most accurate Duluth system was Duluth8, with 1369 correct instances (62%). If the top ranked Spanish system (68%) and Duluth8 are compared, there are 1086 instances (49%) where both are correct, 737 instances (33%) where one or the other is correct, and 402 instances (18%) where neither system is correct.

This is intended as a rule of thumb, and suggests that a fairly substantial percentage of test instances can be resolved by almost any means, and that a hard core of test instances will be very difficult for any method to resolve.

## 5 Acknowledgments

This work has been partially supported by a National Science Foundation Faculty Early CAREER Development award (#0092784).

The Bigram Statistics Package and SenseTools have been implemented by Satanjeev Banerjee.


text:
# Evidence for Cloud Disruption in the L/T Dwarf Transition

## 1 Introduction

The recent discovery of a vast population of free-floating brown dwarfs in the vicinity of the Sun (Kirkpatrick et al., 1999, 2000; Martín et al., 1999; Burgasser et al., 2002b; Geballe et al., 2002) has led to the definition of two new spectroscopic classes, L dwarfs and T dwarfs. These are the first additions to the standard stellar sequence in over 60 years, encompassing objects cooler than the M spectral class. L dwarfs (Kirkpatrick et al., 1999; Martín et al., 1999) are characterized by absorption bands of metal hydrides (e.g., FeH, CrH, and MgH), strong alkali lines (Li, Na, K, Rb, and Cs), and H2O and CO absorption in the near-infrared. Below effective temperatures $`T_{\mathrm{eff}}1500\mathrm{K}`$, CO reduces to $`\mathrm{CH}_4`$ (Fegley & Lodders, 1996), a molecule which has distinct absorption bands at 1.6 and $`2.2\mu \mathrm{m}`$. Objects exhibiting these $`\mathrm{CH}_4`$ features have grossly different spectroscopic and photometric properties as compared to L dwarfs, and are designated T dwarfs (Kirkpatrick et al., 1999, 2000; Burgasser et al., 2002b). The spectroscopic sequence M to L to T is the evolutionary cooling sequence for a typical brown dwarf.

L dwarfs are further characterized by the presence of condensates near or just below their photospheres. These condensates include liquid iron; solid VO; and aluminum, calcium, magnesium, and titanium-bearing minerals such as enstatite ($`\mathrm{MgSiO}_3`$), grossite ($`\mathrm{CaAl}_4\mathrm{O}_7`$), and perovskite ($`\mathrm{CaTiO}_3`$) (Burrows & Sharp, 1999; Lodders, 1999). The opacity of the condensate grains and droplets gives L dwarfs fairly red near-infrared colors, with $`JK_s2`$ (Kirkpatrick et al., 2000). Conversely, condensates in mid- and late-type T dwarfs appear to lie deep below the photosphere. These objects have neutral near-infrared colors, $`JK_s0`$, governed by strong absorption features of $`\mathrm{H}_2\mathrm{O}`$, $`\mathrm{CH}_4`$, and collision-induced $`\mathrm{H}_2`$, particularly at K-band (Burgasser et al., 2002b).

Three classes of models have been developed to account for condensates in brown dwarf atmospheres. The dusty model (Lunine et al., 1989; Chabrier et al., 2000; Allard et al., 2001) assumes that condensates remain dispersed in the gas and in chemical equilibrium throughout the atmosphere. The clear model (Burrows et al., 1997) removes all condensates from the photosphere as they form, presumably through gravitational settling. The cloudy model (Ackerman & Marley, 2001; Marley et al., 2002) incorporates condensate cloud formation based on an assumed sedimentation efficiency, $`f_{\mathrm{rain}}`$. The dusty and clear models are extreme cases of the cloudy model with very low and very high values of $`f_{\mathrm{rain}}`$, respectively.

In this paper, we present an exploratory model of partial cloud clearing in cool dwarf atmospheres derived from the cloudy and clear models. In $`\mathrm{\S }`$2 we describe our model and compare its derived photometry to the near-infrared color-magnitude diagram of late-type dwarfs. This analysis indicates that a rapid clearing of clouds takes place across the L/T transition. In $`\mathrm{\S }`$3 we present further evidence of cloud clearing in the persistance and strengthening of FeH absorption at 9896 Å in the early- and mid-type T dwarfs. We discuss our results in $`\mathrm{\S }`$4.

## 2 Photometric Evidence for Cloud Clearing

### 2.1 The Color-Magnitude Diagram for Late-type Dwarfs

Figure 1 shows the near-infrared color-magnitude diagram for a sample of M, L, and T dwarfs with known parallaxes. Photometry were obtained from the Two Micron All-Sky Survey (Skrutskie et al., 1997, hereafter 2MASS) and parallaxes from the USNO Parallax Program (Dahn et al., 2002) or from Hipparcos (Perryman et al., 1997) in the case of companions to nearby stars. As described in Dahn et al. (2002), late-type M and L dwarfs are increasingly redder and fainter toward later spectral types, but the transition to the bluer T dwarfs is marked by a distinct brightening at J-band, with both the T2 dwarf SDSS J1254-0122 (Leggett et al., 2000) and the T5 dwarf 2MASS 0559-1404111Burgasser (2001) has suggested that 2MASS J0559-1404 may be an unresolved binary brown dwarf. However, even if it were an equal-magnitude double, it would still remain brighter at J-band than the latest L dwarfs. (Burgasser et al., 2000) being brighter at J-band than the latest-type L dwarfs. Note that neither of these objects is expected to be young and hence less compact (R $``$ 1 RJup; Burrows et al. 1997) based on extremely weak or absent H$`\alpha `$ emission (Burgasser et al. 2002a; Kirkpatrick et al., in preparation). Late-type T dwarfs are increasingly bluer in $`JK_s`$ color and fainter at J-band.

Also shown in this figure are color-magnitude tracks for the three atmosphere models described above. The dusty model (Chabrier et al., 2000, dashed line) reproduces the colors of early-type L dwarfs, but continues to redden to colors far greater than those observed for the mid- and late-type L dwarfs. This is due to the continued influence of dust opacity, which mutes the otherwise strong absorption bands of H2O, CH4, and H2. The clear model (solid line, left) reproduces the colors of mid- and late-type T dwarfs, but not those of L dwarfs. Finally, the cloudy model (solid line, right) successfully reproduces the near-infrared color limit for the latest-type L dwarfs, $`JK_s2`$-$`2.5`$ (Kirkpatrick et al., 2000; Marley et al., 2002), and subsequently turns to bluer near-infrared colors. However, the predicted transition occurs over a greater decrease in brightness, $`\mathrm{\Delta }M_J=23`$ mag from $`JK_s=2`$ to $`JK_s=0`$, than is observed. Increasing $`f_{\mathrm{rain}}`$ in this model (i.e., moving toward a more cleared photosphere) allows the color transition to occur at brighter magnitudes, but fails to match the colors and magnitudes of the latest L dwarfs. Therefore, the condensate cloud model, while adequately reproducing the brightness and color evolution of L dwarfs, cannot reproduce the observed rapid transition to the T dwarfs.

Is it possible to resolve this apparent discrepancy between the cloudy model and the observations? All three atmosphere models assume horizontal homogeneity and therefore uniform condensate distribution over the visible disk of the brown dwarf. However, clouds in the planets of the solar system, for instance Jupiter’s upper cloud decks of solid NH3 and NH4SH, display discrete structures of spots and bands caused by tropospheric weather patterns. Breaks between the clouds on Jupiter coincide with bright “hot-spots” at 5 $`\mathrm{\mu m}`$ (Westphal, Matthews, & Terrile, 1974), as light emerges from warmer regions well below the upper cloud decks. Since horizontally non-homogeneous clouds are ubiquitous in planetary atmospheres, it is reasonable to expect that similar features could be present in the condensate cloud layers of brown dwarfs as well (Ackerman & Marley, 2001).

### 2.2 Cloud Clearing Model

Connecting the cloudy and clear tracks in Figure 1 are the predicted fluxes for an alternative model we have developed that allows for partial clearing of condensate clouds. These were constructed by first linearly interpolating between the fluxes of the clear model and a modified cloudy model in which the cloud opacity has been removed, weighing by the fraction of cloud coverage. This combination approximates the effect of partial cloud opacity on the underlying atmospheric pressure-temperature profile. We then linearly interpolated the resulting fluxes with those of the standard cloudy model, again weighing by the fractional cloud coverage, to derive the relative brightness contributions of cloudy and cloud-free regions. Our simple approach assumes that many clouds are distributed randomly across the observed disk.

When compared to the data, we see that the evolution from the red and dusty L dwarfs to the blue and clear T dwarfs can be traced (heavy line) down the cloudy track to an effective temperature $`T_{\mathrm{eff}}1200\mathrm{K}`$, across the 1200 K partly cloudy model, and subsequently down the cloud-free track. SDSS J1254-0122 is a true transition object, lying close to the 1200 K partly cloudy track at roughly 40% cloud coverage. 2MASS J0559-1404, however, appears to have a cloud-free photosphere. The Teff at which the color transition appears to take place is only slightly cooler than that estimated for one of the latest-type L dwarfs, Gliese 584C ($`T_{\mathrm{eff}}=1350\mathrm{K}`$; Kirkpatrick et al. 2000; Dahn et al. 2002). The clearing of clouds at nearly constant Teff explains the apparent brightening of the early-type T dwarfs at 1 $`\mathrm{\mu m}`$, as the emitting layers below the clouds in these objects would be warmer than the layers above the cloud decks in late-type L dwarfs. Note that there is some intrinsic scatter among the empirical data, likely due to the influence of gravity, rotation, and metallicity on cloud coverage222Magnetic activity may also be a source of photometric scatter. Magnetically active low-mass dwarfs have extended cool spots (Alekseev & Gershberg, 1997) which could exhibit stronger H2O and/or CH4 absorption features, resulting in later-type spectral morphologies and photometric colors for a given Teff. This phenomenon would not explain the rapid L/T transition, however, again due to the lack of strong H$`\alpha `$ emission in either SDSS J1254-0122 or 2MASS J0559-1404, as well as the very low photospheric magnetic Reynolds numbers in these cool, neutral atmospheres (Gelino et al., 2002).. Nonetheless, the apparently sudden color transition from L to T clearly cannot be explained by the progressive sinking of a uniform cloud deck with decreasing $`T_{\mathrm{eff}}`$, but is possible through the disruption of condensate clouds.

## 3 Spectroscopic Evidence for Cloud Clearing: FeH Absorption

There is additional evidence that light emitted at $`1\mu \mathrm{m}`$ by early and mid-type T dwarfs is seen through opening holes in condensate cloud decks. Figure 2 plots the optical spectra of a series of late-type L and T dwarfs (Kirkpatrick et al. 1999, 2000; Burgasser et al., in preparation), obtained using the Low Resolution Imaging Spectrograph (Oke et al., 1995) on Keck. Spectral types listed are from Kirkpatrick et al. (1999, 2000) and Burgasser et al. (2002b), and are consistent with a monotonic Teff sequence (Burgasser et al., 2002b). Here we focus on the region centered at the $`9896\mathrm{\AA }`$ FeH band. This feature weakens in the later-type L dwarfs, in concert with the observed disappearance of the $`8962\mathrm{\AA }`$ FeH band (Kirkpatrick et al., 1999). This is consistent with the segregation of liquid iron into condensate clouds, which depletes the atmosphere of iron-bearing gases above the cloud layer. As the clouds settle below the photosphere with decreasing temperature, the amount of detectable FeH drops as well. However, the $`9896\mathrm{\AA }`$ FeH band unexpectedly reappears in SDSS J1254-0122 and is quite pronounced in 2MASS J0559-1404. The retention of FeH absorption from late-type L dwarfs to early- and mid-type T dwarfs is a general trend, as shown in Figure 3. The FeH band clearly weakens from L4 to L8, but persists and even strengthens slightly toward type T5.5, subsequently weakening in the latest T dwarfs.

The plateau in band strength across the L/T transition cannot be explained with a cloud deck that gradually sinks below the photosphere with decreasing Teff. In this case, we would only see the atmosphere above the cloud deck in which FeH is depleted. However, if holes are present in the clouds, deeper and warmer regions in which FeH is not depleted become observable. These deep layers can be detected at $`1\mu \mathrm{m}`$, where $`\mathrm{H}_2\mathrm{O}`$, $`\mathrm{CH}_4`$, and $`\mathrm{H}_2`$ absorption is relatively weak. Since the reappearance of FeH coincides with the turnover in $`JK_s`$ color, the hypothesis of cloud clearing simultaneously explains both phenomena.

## 4 Discussion

The disruption of condensate clouds in brown dwarf atmospheres likely begins when the cloud deck forms within the atmospheric convection zone and becomes subject to coherent vertical motions. At higher effective temperatures, the clouds lie primarily above the radiative-convective boundary and would likely remain spatially uniform, perhaps like the stratospheric hazes on Titan and the giant planets. Evidence for vertical transport has already been seen in cool brown dwarfs, with the detection of CO absorption at $`4.7\mu \mathrm{m}`$ in the T6.5 dwarf Gliese 229B (Nakajima et al., 1995) in excess of its equilibrium abundance (Noll, Geballe, & Marley, 1997; Oppenheimer et al., 1998). Photometric studies of late-type L and T dwarfs also indicate far greater absorption in the $`5\mu \mathrm{m}`$ band than predicted by chemical equilibrium models (Leggett et al., 2002; Reid & Cruz, 2002), likely caused by enhanced CO abundances brought up by vertical flows. Since the principal cloud decks lie within 0.5 pressure scale heights of one another, the convective motion need not be particularly vigorous to disrupt all layers. Note that upwelling brings substantial CO into the photosphere because of its stable bond and hence relatively slow conversion to $`\mathrm{CH}_4`$ (Lodders & Fegley, 2002). The fragility of the FeH bond (1.63 eV versus 11.09 eV for CO) precludes any similar non-equilibrium enhancement of this species above the cloud layer.

The picture of a brown dwarf’s atmospheric evolution presented here has a number of observable consequences. First, the Teff scale between the late-type L and early-type T dwarfs should be very narrow; this has already been surmised by a number of studies (Kirkpatrick et al., 2000; Burgasser et al., 2002b). A narrow temperature scale also implies fewer early-type T dwarfs than mid- or late-type T dwarfs, given the shorter time a brown dwarf spends in this phase of its thermal evolution. The statistics are currently insufficient to adequately test this prediction. We also expect substantial photometric and/or spectroscopic variability in late-type L and early-type T dwarfs, caused by the formation, evolution, and motion (rotational or wind-driven) of cloud holes in their upper atmospheres. Such weather phenomena have been previously cited to explain non-periodic variability in late-type M and early-type L dwarfs (Bailer-Jones & Mundt, 2001; Martín, Zapatero Osorio, & Lehto, 2001). Similar variability could be significant at $`1\mu \mathrm{m}`$ in the L/T transition.

Filling in the color-magnitude sequence of late-type L and early-type T dwarfs will enable us to further characterize the behavior of condensate clouds in brown dwarf atmospheres. Regardless, the observed spectral and photometric properties of brown dwarfs indicate that the planetary-like phenomena of cloud formation and disruption are important in the evolution of free-floating brown dwarfs.

We thank A. Burrows and I. McLean for valuable discussions, M. Putman and K. Zahnle for assistance in the preparation of the manuscript, and G. Chabrier for an electronic version of his dusty models in the 2MASS J and Ks bands. We also thank our anonymous referee for bringing up the issue of magnetic activity in the L/T transition. AB acknowledges support by NASA through Hubble Fellowship grant HST-HF-01137.01 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. MM acknowledges support from NSF grant AST 00-86288 and NASA grant NAG5-9273. Work by KL and DS is supported by NSF grant AST 00-86487 and NASA grant NAG5-4988 respectively. Portions of the data presented herein were obtained at the W. M. Keck Observatory which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. This publication makes use of data from the Two Micron All Sky Survey, which is a joint project of the Unversity of Massachusetts and the Infrared Processing and Analysis Center, funded by the National Aeronautics and Space Administration and the National Science Foundation.