text:
# Magnetism and Phase Separation in the Ground State of the Hubbard Model

## 1 Introduction

Originally proposed for the description of ferromagnetism in transition metals, the Hubbard model hgk

$$H=\underset{i,j}{}\underset{\sigma }{}t_{ij}c_{i\sigma }^{}c_{j\sigma }^{}+U\underset{i}{}n_in_i$$
(1)
is the simplest model to describe the interplay between delocalization or band formation in solids and the local Coulomb correlations. Despite its simplicity, the phase diagram of the Hubbard model (1) reveals a surprising richness. One finds Mott-Hubbard type metal-insulator transitions mit , antiferromagnetism afm , ferromagnetism fm  and incommensurate magnetic phases incomm . More recently, the Hubbard model has also become one of the most promising candidates to describe the low-energy properties and possibly the superconductivity in the high-$`T_c`$ cuprates hightc .

The Hubbard model at half filling $`n=1`$ has been investigated thoroughly using various approximate and exact techniques hub\_pre ; hub\_rigo  and its properties are understood to a large extent. Off half filling, the model is well understood for $`d=1`$ but the situation is less clear in dimensions $`d>1`$. Basically the only rigorous result is due to Nagaoka nagaoka , who proved that a ferromagnetic ground state is possible under certain conditions.

The introduction of the limit $`D\mathrm{}`$ mv89  in principle allows to solve models like the Hubbard model exactly without loosing the competition between kinetic energy and local Coulomb repulsion review\_dmft . This surprising insight has subsequently triggered a large amount of investigations of the infinite dimensional Hubbard model review\_dmft . In addition, $`D\mathrm{}`$ turned out to be a reasonable starting point for weak-coupling expansions PvD\_wk ; PvD\_ps . Within this approach, in addition to the expected magnetic order for a bipartite lattice, phase separation was found for the whole region of the magnetically ordered phase PvD\_ps . Since this result is based on a weak-coupling expansion, it is far from clear whether it holds for finite values of the interaction $`U`$ as well. Results from a numerically exact solution of the Hubbard model in $`D=\mathrm{}`$ based on Quantum Monte-Carlo simulations for example showed no evidence for phase separation jim\_mark , but these calculations were done in the paramagnetic phase and at finite, comparatively high temperatures.

The question whether phase separation in the Hubbard model occurs in a certain parameter regime is of some importance for two reasons. First, from a model theoretical point of view, it is of course interesting to explore the stability of the different possible ordered phases which might be unstable with respect to phase separation. Second, a vicinity to phase separation has been discussed as one of the possible ingredients to the superconductivity in the high-$`T_c`$ cuprates kivelson ; castellani . Moreover, a tendency towards phase separation together with the long-range part of the Coulomb interaction may in principle lead to charge ordered states such as stripe-phases.

Phase separation has long been predicted kivelson  and indeed been observed for the $`t`$-$`J`$ model in $`D=1,\mathrm{\hspace{0.17em}2}`$ dagotto . Since the $`t`$-$`J`$ model for vanishing $`J`$ is connected to the Hubbard model in the limit $`U/t\mathrm{}`$ fulde , additional information about phase separation in the strong coupling limit could thus be obtained. The critical $`J`$, below which phase separation vanishes in the $`t`$-$`J`$ model is typically of the order of $`t`$ dagotto ; this seems to rule out phase separation in the Hubbard model at least for $`U/t\mathrm{}`$ in $`D=2`$. However, the situation has not been clarified yet, because direct inspection of the Hubbard model in the limit $`U/t\mathrm{}`$ leads to contradictory results su ; Tandon99 .

The results for the $`2D`$ Hubbard model for finite $`U`$ available so far have not revealed signs for phase separation dagotto ; becca . However, these results are typically based on Quantum Monte-Carlo or related techniques, which have severe problems in the interesting parameter regime close to half filling and at very low temperatures. Consequently, one either has to restrict oneself to rather small system sizes dagotto  or use further approximations becca . Thus, to our present knowledge a detailed study of the ground state phase diagram of the Hubbard model in the thermodynamic limit and in the vicinity of half filling, comprising weak, intermediate and strong coupling within a non-perturbative approach is not available.

Such an approach is provided by the limit $`D\mathrm{}`$, which allows for in principle exact calculations in the thermodynamic limit for all model parameters, even at $`T=0`$. The price one pays is the neglect of non-local dynamics, which of course is most severe for $`D2`$. Nevertheless, the theory can give valuable information about whether phase separation is possible at all.

The paper is organized as follows. In section 2, we give a brief description of the solution of the Hubbard model in the limit $`D\mathrm{}`$. Results for the phase diagram and the dynamics in the different phases are presented in section 3. The main conclusions of the paper are summarized in section 4.

## 2 Theoretical background

### 2.1 General remarks

The dynamical mean-field theory (DMFT) to exactly solve the Hubbard model in the limit $`D\mathrm{}`$ is based on the work by Metzner and Vollhardt mv89  and is by now well-established review\_dmft . The basic ingredient is that for $`D\mathrm{}`$ the proper single-particle self energy $`\mathrm{\Sigma }(\stackrel{}{k},z)`$ becomes purely local or momentum independent, i.e. $`\mathrm{\Sigma }(\stackrel{}{k},z)\stackrel{D\mathrm{}}{}\mathrm{\Sigma }(z)`$ mv89 ; muha89 . This can be used to map the Hubbard model (1) onto an equivalent quantum impurity problem supplemented by a self-consistency condition review\_dmft . The remaining problem (the solution of a quantum impurity model) is, however, highly nontrivial. Several approximate and numerically exact techniques are currently available review\_dmft ; nrg .

Most of these methods cannot access $`T0`$ or are restricted to the weak-coupling regime of the Hubbard model. The most reliable technique to solve the quantum impurity problem for all interaction strenghts $`U`$ and fillings $`n`$ at $`T=0`$ and low $`T`$ is the numerical renormalization group (NRG) nrg ; krish . Originally, this method was set up to treat the paramagnetic problem only krish , but recent extensions have shown that calculations with a symmetry breaking field are possible with a similar level of accuracy, too nrg\_magn\_costi ; nrg\_magn\_walter . Hence we are able to study magnetically ordered phases directly at $`T=0`$.

In contrast to the standard NRG, a more refined approach has to be used to calculate dynamical quantities in the presence of a magnetic field. This has first been noted by Hofstetter, who observed discrepancies in the magnetization calculated from the spectral functions and the ground state occupation numbers nrg\_magn\_walter . To resolve this problem, he proposed a modification of the standard method sakai  to calculate the spectral function. A more detailed discussion of this technical point and its physical background will be presented elsewhere.

There are in principle two ways to determine the phase boundary between the paramagnetic and a magnetically ordered state. First, one can calculate the susceptibility corresponding to the anticipated order and look for a divergence. Second, one can allow for a proper symmetry breaking in the one-particle Green function and search for the region in parameter space where a solution with broken symmetry becomes stable. Especially for $`T=0`$ the first method is rather cumbersome in general and, by construction, also makes no statement about the thermodynamic stability of phases beyond the critical point.

We thus use the second approach as our method of choice. However, this prohibits the search for incommensurate phases, because only broken symmetries with a commensurate wave vector can be implemented that way. Since we are interested mainly in standard Néel type antiferromagnetic order, the proper way is to introduce an AB-lattice structure and allow for different sublattice magnetizations. The resulting Green function then becomes a $`2\times 2`$ matrix, which within the DMFT has the form review\_dmft

$$𝖦_{\stackrel{}{k}\sigma }(z)=\left(\begin{array}{cc}z+\mu \mathrm{\Sigma }_\sigma ^\mathrm{A}(z)& ϵ_\stackrel{}{k}\\ ϵ_\stackrel{}{k}& z+\mu \mathrm{\Sigma }_\sigma ^\mathrm{B}(z)\end{array}\right)^1,$$
(2)
where $`\stackrel{}{k}`$ is a vector in the magnetic Brillouin zone (MBZ). For the Néel state on an AB-lattice a further simplification arises from the symmetry $`\mathrm{\Sigma }_\sigma ^\mathrm{A}(z)=\mathrm{\Sigma }_{\overline{\sigma }}^\mathrm{B}(z)`$. For the calculation this means that we do not have to solve independent quantum impurity models for the two sublattices, but only one for say sublattice A.

But how to control the filling if the homogeneous solution (homogeneous concerning the charge distribution) turns out to be unstable towards phase separation? Fixing the chemical potential $`\mu `$ is not sufficient as the system will be driven to a filling corresponding to a stable solution, such as $`n=1`$. To enforce a metastable state with finite doping we adopt a procedure which has been already used in calculations for the half-filled Hubbard model in a homogenous magnetic field laloux . The schematic flow diagram of the resulting DMFT self-consistency cycle is shown in Fig. 1. Starting from a paramagnetic solution for the desired doping, a homogenous or staggered

symmetry breaking is introduced and the corresponding effective medium, $`\mathrm{\Gamma }_\sigma (\omega )`$, for the DMFT cycle is determined review\_dmft . Keeping the medium $`\mathrm{\Gamma }_\sigma (\omega )`$ fixed, one now varies the on-site energy of the effective SIAM until the desired doping has been reached. This result is used to obtain a new effective medium, and the procedure is repeated until convergency is reached. It should be noted that, for a metastable state, no true convergency can be reached in the sense that the solution, when iterated further without adjusting the filling properly, will flow again into the phase separated one. Typically, for a converged calculation, the on-site energy between successive DMFT iterations shows a weakly damped oscillatory behavior but does not change by more than 2-3%. In order to minimize errors introduced by these oscillations we average all quantities over several iterations.

To find the correct ground state, we need to calculate the ground state energy

$$\frac{E}{N}=\frac{1}{N}H=\frac{1}{N}H_t+\frac{U}{N}\underset{i}{}n_in_i,$$
(3)
where $`H_t`$ is the kinetic part of the Hamiltonian (1). The expectation value $`n_in_i`$ can be determined within the NRG directly. The quantity $`H_t`$, on the other hand, depends on the phase we are looking at. For the para- and ferromagnetic phases it is simply given by review\_dmft

$$\frac{1}{N}H_t=\underset{\sigma }{}\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑ϵϵ\rho ^{(0)}(ϵ)\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑\omega f(\omega )A_\sigma (ϵ,\omega ),$$
(4)
with $`\rho ^{(0)}(ϵ)`$ the density of states (DOS) for the non-interacting system, $`f(\omega )`$ the Fermi function and

$$A_\sigma (ϵ,\omega )=\frac{1}{\pi }\mathrm{}m\frac{1}{\omega +\mu ϵ\mathrm{\Sigma }_\sigma (\omega +i0^+)}$$

the spectral function of the Hubbard model in the DMFT, i.e. with $`\stackrel{}{k}`$-independent one-particle self energy.

For an antiferromagnetic state with Néel order one has to take into account the AB-lattice structure and the formula becomes review\_dmft

$$\frac{1}{N}H_t=2\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑ϵϵ\rho ^{(0)}(ϵ)\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑\omega f(\omega )B(ϵ,\omega )$$
(5)
instead, with

$$B(ϵ,\omega )=\frac{1}{\pi }\mathrm{}m\frac{1}{\sqrt{\zeta _\sigma (\omega )\zeta _{\overline{\sigma }}(\omega )}ϵ}$$

and $`\zeta _\sigma (\omega )=\omega +\mu \mathrm{\Sigma }_\sigma (\omega +i0^+)`$. Obviously, expression (5) reduces to (4) without magnetic order, i.e. $`\zeta _\sigma (\omega )=\zeta _{\overline{\sigma }}(\omega )`$.

Throughout the rest of the paper we concentrate on results for a simple hypercubic lattice and

$$t_{ij}=\{\begin{array}{cc}t\hfill & \text{if }i\text{}j\text{ are nearest neighbors}\hfill \\ 0\hfill & \text{else}\hfill \end{array}.$$

The resulting DOS $`\rho ^{(0)}(ϵ)`$ then becomes mv89 ; muha89

$$\rho ^{(0)}(ϵ)=\frac{1}{t^{}\sqrt{\pi }}e^{(ϵ/t^{})^2}.$$
(6)
In the following we use $`t^{}=2\sqrt{D}=1`$ as our unit of energy.

### 2.2 Weak-coupling results

Let us briefly review some weak-coupling results as these will be frequently referred to in section 3. Since the hypercubic lattice is a bipartite lattice, one obtains in lowest order, i.e. in Hartree approximation, a transition into a Néel state for any $`U>0`$ at $`T=0`$ below a critical doping $`\delta _c^\mathrm{H}(U)`$. For small $`U0`$ the magnetization $`m`$ as well as the critical doping depend non-analytically on U, i.e. $`m,\delta _c^\mathrm{H}\mathrm{exp}\left(1/(U\rho ^{(0)}(0))\right)/U`$ independent of the dimension.

A quantity of particular interest in the DMFT is the single-particle Green function. The general structure of the Green function in the Néel phase for both Hartree theory and DMFT is given by expression (2), where in the Hartree approximation $`\mathrm{\Sigma }_\sigma (z)`$ reduces to $`\mathrm{\Sigma }_\sigma ^\mathrm{H}(z)=Un_{\overline{\sigma }}=\frac{1}{2}U(n\sigma m)`$ with $`n`$ the filling and $`m`$ the magnetization. The local Green function is obtained from (2) by summing over $`\stackrel{}{k}\mathrm{MBZ}`$, which yields for example for spin up

$$G_{}(\omega )=\frac{\zeta _{}(\omega )}{\sqrt{\zeta _{}(\omega )\zeta _{}(\omega )}}G^{(0)}\left(\sqrt{\zeta _{}(\omega )\zeta _{}(\omega )}\right)$$
(7)
with $`\zeta _\sigma (\omega )=\omega +i0^++\mu \frac{U}{2}n+\sigma \frac{U}{2}m`$ and

$$G^{(0)}(z)=\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑ϵ\frac{\rho ^{(0)}(ϵ)}{zϵ}.$$
(8)
For the further discussion let us define

$$\begin{array}{cc}\omega _{}=\hfill & \frac{U}{2}n\mu \frac{U}{2}m\hfill \\ \omega _+=\hfill & \frac{U}{2}n\mu +\frac{U}{2}m\hfill \end{array}$$

Then, as long as $`\omega \omega _{}`$ or $`\omega \omega _+`$, the radicant in (7) is positive and the resulting DOS can be expressed as

$$\rho _{}(\omega )=\frac{\zeta _{}(\omega )}{\sqrt{\zeta _{}(\omega )\zeta _{}(\omega )}}\rho ^{(0)}\left(\sqrt{\zeta _{}(\omega )\zeta _{}(\omega )}\right).$$

For $`\omega _i<\omega <\omega _+`$, on the other hand, the radicant in (7) is negative, i.e. $`\sqrt{\zeta _{}(\omega )\zeta _{}(\omega )}=i\sqrt{|\zeta _{}(\omega )\zeta _{}(\omega )|}`$. Since for the particle-hole symmetric DOS (6) the Green function $`G^{(0)}(z)`$ defined in (8) for purely imaginary arguments is purely imaginary, too, one finds

$$\rho _{}(\omega )=0,$$

i.e. the DOS has a gap between $`\omega _{}`$ and $`\omega _+`$.

As one approaches $`\omega _{}`$ from below or $`\omega _+`$ from above, it is easy to confirm that

$$\rho _{}(\omega )\{\begin{array}{cc}\sqrt{\frac{Um}{|\omega \omega _{}|}}\rho ^{(0)}(0)\hfill & \omega \omega _{}\hfill \\ \sqrt{\frac{|\omega \omega _+|}{Um}}\rho ^{(0)}(0)\hfill & \omega \omega +\hfill \end{array}.$$
(9)
The corresponding DOS for $`\sigma =`$ has a similar behavior. Here, however, the DOS diverges like $`1/\sqrt{|\omega \omega _+|}`$ at the upper gap edge, and vanishes like $`\sqrt{|\omega \omega _{}|}`$ at the lower one.

In order to determine the thermodynamically stable phase one has to calculate the ground state energy as function of the doping $`\delta =1n`$. The result up to second order in $`U`$ is PvD\_ps

$$E(\delta )E(0)=\frac{U}{2}\delta +\alpha ^\mathrm{H}\mathrm{\Phi }^\mathrm{H}(\delta /\delta _1),$$
(10)
where

$$\mathrm{\Phi }^\mathrm{H}(x)=\{\begin{array}{cc}\frac{1}{2}x\left(1\frac{1}{4}x\right)\hfill & x<1\hfill \\ \frac{1}{4}\left(1+\frac{1}{2}x^2\right)\hfill & x>1\hfill \end{array}$$
(11)
and $`\delta _1`$ is the critical doping for antiferromagnetism in Hartree approximation. The coefficient $`\alpha ^\mathrm{H}`$ is given by $`\alpha ^\mathrm{H}=2\delta _1^2/\rho ^{(0)}(0)`$.

The function $`\mathrm{\Phi }^\mathrm{H}(\delta /\delta _1)`$ appearing in expression (10) leads to the full line in Fig. 3. Apparently, this function is not convex for small $`\delta `$, i.e. the resulting phase is thermodynamically unstable towards phase separation for dopings less than $`\delta _c=\sqrt{2}\delta _1`$. The resulting ground state energy is then obtained from a Maxwell construction, given by the straight dashed line in Fig. 3.

## 3 Results

Let us start with a short

overview of the behaviour at half filling, $`n=1`$. Here, the Néel phase is energetically stable. The variation of the DOS for increasing $`U`$ from $`U=1`$ (full curve) to $`U=6`$ (dashed curve) is shown in Fig. 4. As expected, the DOS for small $`U`$ resembles the form (9) predicted for weak-coupling, i.e. one sees the remnants of the characteristic square-root divergency in the spin up DOS at the lower gap edge and a corresponding power law at the upper gap edge. These characteristic features however vanish rapidly with increasing $`U`$, and already for $`U=3`$ the DOS mainly consists of the Hubbard peaks at $`\omega =+U/2`$ and $`\omega =U/2`$ for $`\sigma =`$ and $`\sigma =`$, respectively; reminiscent of the behavior expected for the Mott-Hubbard insulator, where only the incoherent charge excitation peaks at high energies are present mit ; review\_dmft ; jtp92 ; bul99 . Note that neither from the spectra in Fig. 4 nor from the behaviour of the magnetic moment in the inset of Fig. 4 one can infer that at $`U4.1`$ the Mott-Hubbard metal-insulator transition occurs in the paramagnetic state jtp92 ; bul99 .

Keeping $`U`$ fixed at $`U=3`$ and increasing $`\delta `$ leads to the spectra shown in Fig. 5. Quite interestingly, the typical weak-coupling characteristics reappear in the spectra for small doping and are still recognizable for $`\delta =13\%`$. Note also that upon variation of doping and hence of the magnetization the spectra are not shifted in the same way as in Hartree theory. Instead, the dominant effect is a strong redistribution of spectral weight from the Hubbard bands to the Fermi level. Eventually, in the paramagnetic phase one recovers the well-known three peak structure of the doped Hubbard model in the DMFT review\_dmft .

The evolution of the spectra both at and off half filling can be understood within a simple picture. In Fig. 6 we show a sketch of

the Hartree bandstructure of the Hubbard model in the Néel state, which has two branches in the MBZ and a gap of width $`U|m|`$ between them. If, on the other hand, we inspect the paramagnetic solution, one for example finds at half filling and for small values of $`U`$ a Fermi liquid with quasiparticles defined on an energy scale larger than $`U|m|`$. This situation is indicated by the arrow labeled (a) on the left side of Fig. 6. Here we expect, and indeed find for the antiferromagnetic solution (see full curve in Fig. 4), a DOS that shows the characteristic van-Hove singularities of Fig. 2. Increasing $`U`$ eventually leads to a situation, where the energy scale for the quasiparticles in the paramagnetic state is finite but much smaller than $`U|m|`$ (arrow (b) in Fig 6). The self-energy in the energy region of the van-Hove singularities then has a large imaginary part and will completely smear out the characteristic structures. Further increasing $`U`$ into the Mott-Hubbard insulator will then not change the picture qualitatively, explaining the similarity between the curves for $`U=3<U_{\mathrm{MIT}}`$ and $`U=6>U_{\mathrm{MIT}}`$ in Fig. 4. With finite doping, we move the chemical potential into, e.g., the lower band; this means that even for a relatively small quasiparticle energy scale one again sees the van-Hove singularities at the band edge, which results in the well defined singularities in the spectra for small doping in Fig. 5.

From the occupation numbers $`n_\sigma `$ obtained after convergence of the DMFT calculation one can

calculate the magnetization per electron, $`m=(n_{}n_{})/(n_{}+n_{})`$, as function of the doping $`\delta `$. The results for $`U=1`$ and $`U=3`$ are shown in Fig. 7a and b together with a fit to a power law

$$m(\delta )=m_0\left|1\frac{\delta }{\delta _c^{\mathrm{AF}}}\right|^\nu .$$
(12)
The resulting fit parameters are summarized in Table 1. As expected for a mean-field theory, the value for the critical exponent is $`\nu =1/2`$.

Finally, with the converged DMFT self-energy $`\mathrm{\Sigma }_\sigma (z)`$ we can calculate the expectation value $`H/N`$ according to Eqs. (3) and (4) respectively (5) for the paramagnetic and antiferromagnetic phase. The results for the characteristic function

$$\mathrm{\Phi }(\delta )=E(\delta )+\frac{U}{2}\delta E_{\mathrm{mag}}(0)$$

are summarized in Fig. 8a and b. In Fig. 8 the energies of the antiferromagnetic phase are represented by the circles, those of the paramagnetic phase by squares. The full lines interpolating the antiferromagnetic data are fits to the function

$$\mathrm{\Phi }(\delta )=\alpha \mathrm{\Phi }^\mathrm{H}(\delta /\delta _1)+\gamma \left(\frac{\delta }{\delta _1}\right)^3$$
(13)
with $`\mathrm{\Phi }^\mathrm{H}(x)`$ according to (11). The fit parameters are summarized in Table 1. The use of the function $`\mathrm{\Phi }^\mathrm{H}(x)`$ in

(13) is motivated by the results of van Dongen PvD\_ps . The lines interpolating the paramagnetic data are meant as guides to the eye only. The dotted vertical lines denote the value $`\delta _c^{\mathrm{AF}}`$ as obtained from Fig. 7.

The antiferromagnet obviously has the lower energy as compared to the paramagnet in the region $`0\delta \delta _c^{\mathrm{AF}}`$. However, in both cases $`U=1`$ and $`U=3`$ we find a clear non-convex behavior in $`E(\delta )`$ in that region, i.e. the aforementioned signature of an instability towards phase separation. The true ground state energy as function of $`\delta `$ is obtained again via a Maxwell construction, leading to the dashed lines in Fig. 8 and the values $`\delta _c^{\mathrm{PS}}`$ given in Table 1. Note that in both cases $`\delta _c^{\mathrm{AF}}\delta _c^{\mathrm{PS}}`$ within the accuracy of the fitting procedure.

While for $`U=1`$ the function $`\mathrm{\Phi }(\delta )`$ nicely follows the weak-coupling prediction (10) with renormalized constant $`\alpha `$ one finds a sizeable contribution $`\delta ^3`$ for $`U=3`$. This additional term results in a much weaker non-convex behavior of $`E(\delta )`$ for $`U=3`$.

For values $`U>4`$ we were not able to find a stable solution with Néel order and well-defined doping $`\delta >0`$, although for $`\delta <\delta _c(U)`$ the paramagnetic phase becomes unstable. However, the numerical calculations rather produce a cycle encompassing a range of fillings instead of one solution with definite filling here. It might be interesting to note that at least each of the fillings in this cycle has a unique magnetization associated with it and that all spectra in this cycle correspond to an insulator. Currently it is neither clear what type of magnetic solution we find here, nor whether the breakdown of the Néel solution is a true physical effect or due to numerical problems. Since at half filling the Néel state is present at these values of $`U`$, incommensurate structures or again a phase separated state seem to be possible.

For values of $`U`$ beyond $`U_c25`$ yet another magnetic phase appears, namely the ferromagnet. The existence of this phase has been observed in the case of a hypercubic lattice and $`U=\mathrm{}`$ obermeier  and for a generalized fcc lattice ulmke  before. Since these calculations had to be done at finite and comparatively high temperatures, questions regarding the ground state magnetization and, especially in the case of a hypercubic lattice, the actual extent of the ferromagnetic phase in $`(\delta ,U)`$ space could not be discussed satisfactorily.

As an example for the ferromagnetic phase at $`T=0`$ Fig. 9 shows the magnetization per electron, $`m(\delta )=(n_{}n_{})/(n_{}+n_{})`$ as function of doping (Fig. 9a) and the local DOS for two dopings (Fig. 9b)

($`U=50`$). The data for $`m(\delta )`$ in Fig. 9a are fitted to the function

$$m(\delta )=m_0\sqrt{1\left(\frac{\delta }{\delta _c}\right)^\nu },$$
(14)
and the result is given by the dotted line. The parameters for the fit are $`m_0=1`$, $`\delta _c=14.6`$% and $`\nu =2.75`$. While for $`\delta \delta _c`$ the typical mean-field behavior, i.e. $`m(\delta )\sqrt{1\frac{\delta }{\delta _c}}`$, is obtained, the result for $`\delta 0`$ is rather unconventional, viz $`m(\delta )1\frac{1}{2}\left(\frac{\delta }{\delta _c}\right)^{2.75}`$. This fit assumes that a fully polarized state is only reached as $`\delta 0`$ naga\_stab . Note, however, that the numerical results for the magnetization $`m(\delta )`$ for small $`\delta `$ are also consistent with a fully polarized ferromagnet at finite $`\delta `$.

It is also quite apparent from the DOS in Fig. 9b, that the ferromagnetism found here cannot be understood on the basis of the typical Stoner theory. In contrast to the shifts of the spectrum expected in the latter, we observe a strong redistribution of spectral weight instead, but retain otherwise the typical structures due to the strong correlations. Only in the case $`\delta 0`$ the spectrum again resembles that of a free system for the (almost completely polarized) majority spins. The minority spins become strongly depleted below the Fermi energy, the spectral weight can be found almost completely in the upper Hubbard band situated around $`\omega U/2`$ (not shown in the figure). Nevertheless, we observe a tiny resonance just above the Fermi energy even as $`\delta 0`$.

## 4 Summary and conclusion

In this paper, we used the dynamical mean-field theory together with Wilson’s numerical renormalization group to investigate the ground-state properties of the Hubbard model on a hypercubic lattice with nearest-neighbor hopping both at and off half-filling. While at half-filling the ground-state is antiferromagnetic for all $`U>0`$, at least for the weak and intermediate coupling regime this magnetic order can only be realized in a phase-separated state for any finite doping, thus supporting and extending earlier weak-coupling predictions.

The mapping of the Hubbard model for large $`U`$ to an antiferromagnetic $`t`$-$`J`$ model strongly suggests the dominance of antiferromagnetism in the ground state. The results for the Hubbard model in this paper show, however, that the type of magnetic order for intermediate values of the Coulomb repulsion $`U`$ off half filling is still an open issue; furthermore, the role of phase-separation (which is observed for $`U3`$) has still to be clarified for larger values of $`U`$.

The results are summarized in the schematic $`(\delta ,U)`$ ground state phase diagram of Fig. 10.

To allow the inclusion of all values $`0U<\mathrm{}`$, the ratio $`U/(1+U)`$ is used on the abcissa. Close to half filling we find a phase separated Néel antiferromagnet up to a certain value of $`U<4`$. The magnetization as function of doping follows a typical mean-field behavior in all cases studied and the spectra show the characteristic van-Hove singularities of the band structure in the Néel state in cases where the characteristic energy scale of the paramagnet is large enough. Most important is the observation that, as typical for correlated ordered systems, the spectra are not strongly shifted, as e.g. predicted by Hartree theory, but rather show a strong redistribution of spectral weight.

For values of $`U>4`$ the system shows the tendency towards a magnetic instability, which could not be further identified due to technical problems in the solution of the DMFT self-consistency. However, we at least can exclude ferromagnetism here and a speculative possibility would be the occurence of incommensurate phases or magnetic phases with additional charge order. While the former cannot be addressed easily within the present method, the latter possibility will be investigated further.

At very large values of $`U>25`$, there is a region of ferromagnetism, extending between $`0<\delta <30`$% as $`U\mathrm{}`$ obermeier . For a fixed value of $`U`$, the magnetization per electron in the ferromagnetic state shows the tendency to saturate; from the numerical data it is of course impossible to reliably conclude whether the system is fully polarized at a finite $`\delta `$ already or only as $`\delta 0`$. The data are consistent with both scenarios, but the latter is supported by analytical treatments of the case $`\delta 0`$. As in the case of the antiferromagnet, the spectrum shows a rather strong redistribution of spectral weight, not simply a shift of the features, as would be expected from Stoner theory.

The phase diagram shows a peculiarity, which has alread been pointed out by Obermeier et al. obermeier . In the region of very large $`U`$ and $`\delta 0`$ there exists the possibility of a direct transition between the “antiferromagnetic” phase and the ferromagnet. As at the point $`(\delta ,U)=(0,\mathrm{})`$ all possible spin configurations are degenerate, one can speculate how the phase diagram looks like as $`(\delta ,U)(0,\mathrm{})`$ priv\_com . Generic

possibilities are sketched in Fig. 11. There can either be a direct transition between the two phases (Fig. 11a), which quite likely would then be of first order, a gap filled by a paramagnetic phase (Fig. 11b) or a new phase, e.g. a ferrimagnet interpolating between the two extremes. The a priori exclusion or verification of any of these structures is, without a detailed knowledge of the analytic behaviour of the relevant quantities as function of $`(\delta ,U)`$ in the vicinity of $`(\delta ,U)=(0,\mathrm{})`$, not possible.

While there is a consensus about the magnetic properties of the Hubbard model in a qualitative sense, the direct inspection of details still reveals unexpected surprises. Even within the DMFT, where one can safely state that the paramagnetic phase diagram including the Mott-Hubbard metal insulator transition is now understood, the investigation of the magnetic properties is far from complete. Obvious open questions are the magnetic properties at intermediate values of $`U`$ and the behavior when the antiferromagnetic and ferromagnetic phases meet. Furthermore, the behaviour of the magnetic phases, especially the spectral properties in the ordered phases, in the presence of a frustration due to longer range hopping, has not been addressed yet. This might be of some interest regarding the question how the first order Mott-Hubbard transition manifests itself in the magnetically ordered state. Work along these lines is in progress.

###### Acknowledgements.

We acknowledge useful discussions with P.G.J. van Dongen, M. Jarrell, J. Keller, D. Logan, M. Vojta and D. Vollhardt. This work was in part supported by the SFB 484 Kooperative Phänomene im Festkörper and the Leibniz Computer center.


format: text

text:
# Stellar Velocity Dispersion and Black Hole Mass in the Blazar Markarian 501

## 1 Introduction

The tight correlation recently discovered between stellar velocity dispersion and black hole mass in nearby galaxy bulges (Ferrarese & Merritt, 2000; Gebhardt et al., 2000a) has become the key to our understanding of black hole demographics, as well as a new tool for probing the evolution of galaxies and quasars. An equally important aspect of the $`M_{}\sigma _{}`$ correlation is its predictive power. While dynamical mass measurements are observationally difficult and only possible for a limited number of galaxies, the $`M_{}\sigma _{}`$ relation makes it possible to obtain a black hole mass estimate good to $`40\%`$ accuracy or better from a single measurement of bulge velocity dispersion.

Having such a straightforward method to estimate black hole masses is a tremendous boon to studies of active galactic nuclei (AGNs), because $`M_{}`$ is a fundamental parameter affecting the energetics and emission properties of AGNs. Gebhardt et al. (2000b) and Ferrarese et al. (2001) have recently shown that black hole masses derived for Seyfert galaxies via the $`M_{}\sigma _{}`$ relation are consistent with masses determined by reverberation mapping, providing added confidence that AGNs do follow the same correlation as inactive galaxies. For some classes of AGNs, using the correlations between $`M_{}`$ and host galaxy properties may be the *only* reliable way to estimate the central masses. BL Lac objects fall in this category, since more direct methods for determining $`M_{}`$ (stellar dynamics or reverberation mapping) cannot be applied.

Markarian 501 is one of the nearest known BL Lac objects. Its redshift is $`z=0.0337`$ (Ulrich et al., 1975), corresponding to a distance of 144 Mpc for $`H_0=70`$ km s-1 Mpc-1. The host of Mrk 501, also known as UGC 10599, is a giant elliptical that appears morphologically normal in ground-based images (Hickson et al., 1982; Abraham, McHardy, & Crawford, 1991; Stickel, Fried, & Kühr, 1993; Wurtz, Stocke, & Yee, 1996; Nilsson et al., 1999). Stellar absorption lines are visible in the nuclear spectrum, along with weak emission lines (Ulrich et al. 1975; Stickel et al.1993). Mrk 501 is also a $`\gamma `$-ray source, and is one of the few extragalactic objects from which TeV emission has been detected (Quinn et al., 1996; Bradbury et al., 1997).

In this *Letter*, we present new optical spectra of Mrk 501. The Ca II triplet lines are clearly detected, albeit substantially diluted by the nonthermal continuum. The stellar velocity dispersion measured from these lines is $`372\pm 18`$ km s-1. Using this result in conjunction with the $`M_{}\sigma _{}`$ correlation, we derive an estimate of the mass of the black hole in Mrk 501.

## 2 Observations and Reductions

Observations were obtained with the Palomar Hale 200-inch telescope and Double Spectrograph (Oke & Gunn, 1982) on 2001 June 24 UT. On the red side, we used a 1200 grooves mm-1 grating blazed at 9400 Å, and a wavelength setting giving coverage of 8435–9075 Å at 0.635 Å pixel-1. For the blue side, we used a 600 grooves mm-1 grating blazed at 4000 Å, covering 4200–5950 Å at 1.72 Å pixel-1. A 2″-wide slit was used for all observations. The total exposure time for Mrk 501 was 3.5 hours, broken into individual 30-minute exposures.

The airmass ranged between 1.01 and 1.20 during the observations, and the slit was held at a fixed position angle of 49°. There were intermittent thin clouds during the observations, and the seeing was typically $`1\stackrel{}{\mathrm{.}}5`$. We observed Mrk 501 again on June 26, but due to thick clouds the data were of very poor quality and we do not include these exposures in our analysis. We also observed 24 stars as velocity templates during the observing run. The stars included giants, subgiants, and dwarfs with spectral types between F3 and K5.

Each frame was bias-subtracted and flat-fielded using standard techniques. Spectral extractions were performed using a width of 3$`\stackrel{}{\mathrm{.}}`$7, corresponding to 8 pixels on the red side and 6 pixels on the blue side. The extractions were flux-calibrated using observations of the standard star BD +174708, and wavelength-calibrated using exposures of a He-Ne-Ar lamp, with a final linear shift based on the wavelengths of night sky emission lines. Telluric absorption bands longward of 8925 Å were removed by dividing by the normalized spectrum of the standard star, but the corrected spectra are extremely noisy longward of 8950 Å. A few narrow, weak telluric features with equivalent widths of $`<0.04`$ Å remain in the object spectra in the region 8800-8925 Å. These features could not be removed by standard star division since the standard star has intrinsic absorption features in this wavelength range that are comparable to, or stronger than, the telluric features.

The spectrum of Mrk 501 is displayed in Figure 1, along with the nearby elliptical galaxy NGC 4278 as a comparison. Several stellar absorption features are clearly visible in the Mrk 501 spectrum, despite the strong dilution by nonthermal emission. These include the Ca II triplet lines at $`\lambda \lambda 8498,8542,8662`$ Å, and on the blue side, the G band and Mg I*b*. Weak emission lines of \[O III\] $`\lambda \lambda 4959,5007`$ are also detected in Mrk 501, but H$`\beta `$ emission is absent. The Ca II $`\lambda 8542`$ absorption line has an equivalent width of only $`1`$ Å (in comparison with 3–4 Å in our K-giant stellar templates), indicating that the spectrum is dominated by nonthermal emission.

## 3 Measurement of Velocity Dispersion

To measure the velocity dispersion, we used a template fitting routine based on the method described by van der Marel (1994), which performs a direct fit of broadened stellar templates to the galaxy spectrum. The object and template spectra are first rebinned to a wavelength scale that is linear in $`\mathrm{log}\lambda `$. The template is broadened by convolution with a Gaussian, and its spectral shape is adjusted with an additive dilution, representing the nonstellar continuum, and multiplication by a quadratic polynomial, to allow the continuum shape of the template to match the object. The galaxy’s radial velocity is an additional free parameter in the fit. This procedure is repeated using a downhill simplex search algorithm (Press et al., 1988), allowing all parameters to vary freely, to find the parameter values that minimize $`\chi ^2`$, yielding the best-fitting Gaussian width $`\sigma _{}`$. The fitting routine was tested with spectra of velocity dispersion standard galaxies from the catalog of McElroy (1995) that were observed during the same run, and in all cases our results agreed well with the previously measured values.

The assumption of a Gaussian for the line-of-sight velocity distribution is a reasonable choice given the relatively low signal-to-noise ratio of the Mrk 501 spectrum; our data do not have the sensitivity needed to detect higher-order moments of the velocity distribution. The high-velocity wings on the velocity profiles due to the central black hole should not be significant given the large spectroscopic aperture, which corresponds to 1.3 kpc $`\times `$ 2.0 kpc at the distance of Mrk 501.

For the red side spectrum, we chose a fitting region just large enough to contain the Ca lines, 8475–8680 Å. The spectrum becomes extremely noisy beyond 8680 Å and we truncated the fit at that wavelength; the red wing of the $`\lambda 8662`$ line is essentially lost in this portion of the spectrum due to the telluric water absorption band, even after division by the standard star. We excluded a few remaining weak telluric features from the fit, including one located in the blue wing of the $`\lambda 8542`$ absorption feature. The best-fitting template is HD 188056, a K3III star (Figure 2). With this template, the derived velocity dispersion is $`\sigma _{}=372\pm 7`$ km s-1, and the standard deviation among all templates observed during the run is 17 km s-1. Since template mismatch is the dominant source of uncertainty in the measurement, we combine the template-matching and fitting uncertainties in quadrature for a final result of $`\sigma _{}=372\pm 18`$ km s-1.

The Ca triplet region is expected to be less sensitive to stellar population variations than the blue spectral region (e.g., Dressler, 1984), so the fit to the Ca lines should give the best estimate of the velocity dispersion. We also performed fits to several portions of the blue spectrum, and found that the best results were obtained for a region just redward of the Mg I*b* line. The Mg I*b* line itself did not match the template spectra well, however, and it was excluded from the final fits. Over the range 5200-5600 Å, the best-fitting template was the K3III star HD 125560, giving $`\sigma _{}=386\pm 9`$ km s-1, in reasonable agreement with the red side results. However, several of the template stars gave very poor fits, and the standard deviation among all templates for this region was 81 km s-1.

## 4 The Mass of the Black Hole in Mrk 501

If we assume that Mrk 501 follows the $`M_{}\sigma _{}`$ relation defined for nearby galaxies, then it is straightforward to convert our measurement of $`\sigma _{}`$ into an estimate of $`M_{}`$. This assumption is supported by the fact that the Mrk 501 host galaxy appears morphologically undisturbed, without any indications of recent major merger or interaction events that might cause the galaxy to deviate from the correlation.

The most recent fits to the $`M_{}\sigma _{}`$ relation are given by Kormendy & Gebhart (2001), who find $`M_{}\sigma _{}^{3.65}`$, and by Merritt & Ferrarese (2001), who find $`M_{}\sigma _{}^{4.72}`$. Given the disagreement between these two groups, the most conservative approach is to apply both relations and derive a range of possible black hole masses for Mrk 501. Another issue is the fact that Gebhardt et al. (2000a) use the luminosity-weighted velocity dispersion within the half-light radius, while Ferrarese & Merritt (2000) base their fit on the “central” velocity dispersion measured within an aperture of radius $`r_e/8`$. For our purposes, the distinction between these two choices is not large; Gebhardt et al. (2000a) show that variations in luminosity-weighted velocity dispersion for apertures of different size rarely exceed 10% for $`r<4r_e`$. One possibility would be to use the prescription of Jørgensen, Franx, & Kjærgaard (1995) to normalize our measured aperture dispersion to $`r_e`$ or $`r_e/8`$. However, this would introduce additional uncertainty as the value of $`r_e`$ for Mrk 501 is not well constrained by existing ground-based imaging data. Values for $`r_e`$ compiled from the literature by Nilsson et al. (1999) range from 9″ to 20″. Thus, we use our measured velocity dispersion without applying any aperture corrections.

Propagating the uncertainty in $`\sigma _{}`$ as well as the uncertainties in the coefficients of the $`M_{}\sigma _{}`$ relation, the black hole in Mrk 501 has $`M_{}=(1.3\pm 0.4)\times 10^9`$ $`M_{}`$ with the Kormendy & Gebhart (2001) fit, or $`(2.4\pm 1.0)\times 10^9`$ $`M_{}`$ with the Merritt & Ferrarese (2001) fit, for a $`1\sigma `$ allowed range of $`(0.93.4)\times 10^9`$ $`M_{}`$. Given this result, the projected radius of the gravitational sphere of influence ($`r_g=GM/\sigma ^2`$) of the black hole in Mrk 501 is expected to be $`0\stackrel{}{\mathrm{.}}040\stackrel{}{\mathrm{.}}16`$. This raises the possibility that a genuine stellar-dynamical mass measurement could be done with the *Hubble Space Telescope* (*HST*), if $`M_{}`$ is in the upper half of this uncertainty range.

The main caveat to this result is that the upper end of the $`M_{}\sigma _{}`$ correlation is rather uncertain at present. In the range $`M_{}>10^9`$ $`M_{}`$, there are as yet only three galaxies with dynamical $`M_{}`$ measurements derived from *HST* data. *HST* programs currently in progress will rectify this situation during the next few years with many additional measurements, giving better constraints on the amount of intrinsic scatter in the $`M_{}\sigma _{}`$ relation in this mass range.

## 5 Discussion

The correlation between black hole mass and galaxy bulge luminosity can also be used to obtain black hole mass estimates. The Mrk 501 host galaxy has a total $`B`$ magnitude of 14.4 (Stickel et al. 1993). Correcting for Galactic extinction of $`A_B=0.084`$ mag (Schlegel, Finkbeiner, & Davis, 1998) and applying a $`K`$-correction of 0.16 mag (Pence, 1976), the galaxy has $`M_B=21.6`$ mag for $`H_0=70`$ km s-1 Mpc-1. From the fit to the $`M_{}L_{\mathrm{bulge}}`$ correlation given by Kormendy & Gebhart (2001), the expected mass is $`M_{}=6.1\times 10^8`$ $`M_{}`$. While this value is outside our $`1\sigma `$ uncertainty range on $`M_{}`$ from the $`M_{}\sigma _{}`$ relation, the scatter in the $`M_{}L_{\mathrm{bulge}}`$ correlation is more than an order of magnitude in $`M_{}`$ at fixed $`L_{\mathrm{bulge}}`$ (Kormendy & Gebhart, 2001). Also, the host galaxy luminosity is rather uncertain; the values compiled from the literature by Nilsson et al. (1999) show a range of $`0.5`$ mag among different authors. (We note that for the adopted luminosity, the galaxy lies very close to the mean Faber-Jackson relation for nearby ellipticals.) Thus, the level of disagreement between these two $`M_{}`$ estimates is not a cause for concern; M87 is an outlier by a similar amount in the $`M_{}L_{\mathrm{bulge}}`$ correlation. The $`M_{}\sigma _{}`$ relation is a much more reliable predictor of $`M_{}`$ because its intrinsic scatter is much smaller, less than $`40\%`$ (Gebhardt et al., 2000a) and possibly near zero (Ferrarese & Merritt, 2000).

Some previous estimates of the black hole mass in Mrk 501 have been obtained by examination of the spectrum and variability of its nonthermal emission. For example, Fan, Xie, & Bacon (1999) derive $`M_{}`$ by combining the variability timescale (from which they estimate the physical size of the emitting region) with the assumption that the $`\gamma `$-rays in blazars are produced at $`200`$ Schwarzschild radii from the black hole. For Mrk 501 and seven other $`\gamma `$-ray blazars, they find central masses in the range $`(17)\times 10^7`$ $`M_{}`$. Rieger & Mannheim (2000) interpret the periodic behavior observed in X-ray and $`\gamma `$-ray light curves as evidence for a binary black hole in Mrk 501, and propose a geometric model for a jet originating from the less massive black hole to estimate $`M_{}10^8`$ $`M_{}`$ and $`(442)\times 10^6`$ $`M_{}`$ for the two components of the binary. Recently, Wang, Xue, & Wang (2001) have devised a method to estimate $`M_{}`$ in blazars from the peak luminosity and peak frequency. Using data for a large sample (but not including Mrk 501), they show that their method implies masses of $`10^{911}`$ $`M_{}`$ for low-frequency peaked blazars, and $`10^{58}`$ $`M_{}`$ for high-frequency peaked blazars.

However, these methods have not been calibrated for galaxies whose black hole masses have been measured dynamically, so it is difficult to assess their accuracy. Black hole masses below $`10^8`$ $`M_{}`$ in BL Lac objects would be difficult to reconcile with the $`M_{}`$$`L_{\mathrm{bulge}}`$ correlation, since *HST* imaging has conclusively demonstrated that most BL Lac hosts are luminous ellipticals (Urry et al., 2000). In addition, masses of $`10^8`$ $`M_{}`$ would pose a problem for unified models of radio-loud AGNs, in which BL Lacs are interpreted as radio galaxies oriented with the radio jet along our line of sight (e.g., Urry & Padovani, 1995). Nearby FR I radio galaxies, which would presumably appear as BL Lac objects if viewed pole-on, have black holes with masses in the range $`(0.33)\times 10^9`$ $`M_{}`$ (e.g., Harms et al., 1994; Ferrarese, Ford, & Jaffe, 1996). Thus, for Mrk 501, the central mass derived from either $`\sigma _{}`$ or $`L_{\mathrm{bulge}}`$ is consistent with the range of masses expected from unified schemes.

The Eddington ratio of the active nucleus in Mrk 501 can be determined if its bolometric luminosity is known. Only indirect estimates of $`L_{\mathrm{bol}}`$ are possible, since the observed nonthermal emission is highly beamed. One way to estimate $`L_{\mathrm{bol}}`$ is by comparison with FR I radio galaxies, the unbeamed counterparts of BL Lac objects. Bolometric luminosities for four nearby FR I radio galaxies (M84, M87, NGC 6251, and NGC 4261) are available from Ho (1999), and nuclear emission-line measurements for these objects are listed by Shuder & Osterbrock (1981) and Ho, Filippenko, & Sargent (1997). For these four galaxies, the ratio $`L`$(H$`\alpha `$+\[N II\])/$`L_{\mathrm{bol}}`$ has a mean value of $`4.4\times 10^3`$ with a scatter of only $`25\%`$ over an order of magnitude range in $`L_{\mathrm{bol}}`$. Mrk 501 has $`L`$(H$`\alpha `$+\[N II\]) $`=1.1\times 10^{41}`$ erg s-1 (Stickel et al. 1993). If its $`L`$(H$`\alpha `$+\[N II\])/$`L_{\mathrm{bol}}`$ ratio is similar to these local FR I galaxies, as might be expected from unified models, then it has $`L_{\mathrm{bol}}2.4\times 10^{43}`$ erg s-1 and $`L/L_{\mathrm{Edd}}10^4`$. Mrk 501 is evidently an extremely sub-Eddington accretor; this strengthens previous indications that nearby TeV-emitting blazars are intrinsically weak AGNs (e.g., Cavaliere & Malquori, 1999).

The correlations between $`M_{}`$ and host galaxy properties make it possible to consider new tests of AGN unification scenarios, by examining the distributions of black hole masses and Eddington ratios between different classes of AGNs. Such tests have been performed using the host galaxy luminosity to estimate $`M_{}`$ (Treves et al., 2001), but a comparison based on stellar velocity dispersion would provide improved constraints on the central masses, and new probes of the demographics of AGN classes such as FR I and FR II radio galaxies and blazars. Also, it would be worthwhile to test whether the black hole masses differ systematically between high- and low-frequency peaked blazars, as advocated by Wang et al. (2001). It should be possible to measure stellar velocity dispersions in several other low-redshift BL Lac objects. At recession velocities beyond 12,000 km s-1, the Ca $`\lambda 8542`$ line will be redshifted into telluric H2O absorption bands and difficult to measure accurately, but the Mg I*b* spectral region will still be observable. We are beginning a spectroscopic survey of additional nearby BL Lac objects in order to address these questions.

Research by A.J.B. is supported by NASA through Hubble Fellowship grant #HST-HF-01134.01-A 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, and by a postdoctoral fellowship (through September 2001) from the Harvard-Smithsonian Center for Astrophysics. Research by W.L.W.S. is supported by NSF grant AST-9900733. We thank Tom Matheson for writing some of the software used for spectroscopic data reduction, Todd Small for assistance during the observing run, and Roger Blandford for illuminating conversations.


text:
# A Search for Variability in the Active T dwarf 2MASS 1237+6526

## 1 Introduction

A previously undetected population of T dwarfs, brown dwarfs exhibiting CH4 absorption bands at 1.6 and 2.2 $`\mathrm{\mu m}`$ (Kirkpatrick et al., 1999; Burgasser et al., 2002a), have recently been identified by various deep optical and near-infrared surveys (Strauss et al., 1999; Burgasser et al., 1999, 2000b, 2000c, 2002a; Cuby et al., 1999; Tsvetanov et al., 2000; Leggett et al., 2000; Geballe et al., 2002). One of these objects, 2MASSI J1237392+652615 (Burgasser et al., 1999, hereafter, 2MASS 1237+6526), identified in the Two Micron All Sky Survey (Skrutskie et al., 1997, hereafter 2MASS), was found to have H$`\alpha `$ in emission, with a relative luminosity $`\mathrm{log}_{10}`$(L/Lbol) = $``$4.3 (Burgasser et al., 2000a). While many low-mass stars and young brown dwarfs are seen to exhibit this feature, likely generated by chromospheric magnetic activity, the frequency and relative luminosity of H$`\alpha `$ emission drops rapidly beyond spectral types M7 V (Kirkpatrick et al., 2000; Gizis et al., 2000). Indeed, at the time of its detection, the emission of 2MASS 1237+6526 was unique amongst objects later than type L5 V. Kirkpatrick et al. (2002) have recently reported weak H$`\alpha `$ emission in the bright T2 V SDSSp J125453.90-012247.4 (Leggett et al., 2000, hereafter SDSS 1254-0122).

In order to explain the activity of 2MASS 1237+6526, Burgasser et al. (2000a) have hypothesized that this object may actually be an interacting brown dwarf binary system. In this scenario, the binary mass ratio and orbital separation are such that the lower-mass secondary fills its Roche lobe and steadily loses mass to the primary. For a 70 MJup111We adopt the definitions MJup = 1.90$`\times `$1030 g = 9.55$`\times `$10-4 M, RJup = 7.15$`\times `$109 cm = 0.10 R (Cox, 2000; Tholen, Teiffel, & Cox, 2000). primary and mass ratio $`q`$ $``$ $`M_2/M_1`$ $`>`$ 0.07 (i.e., $`M_2`$ $`>`$ 5 MJup), sustained mass loss can occur for physical separations $`a`$ $``$ 6 RJup and orbital periods $`p`$ $``$ 5.5 hours. The close orbit required for this scenario currently rules out direct imaging of the binary pair ($`a`$ $`<`$ 0$`\stackrel{}{\mathrm{.}}`$0002 assuming a distance of 14 pc; Burgasser et al. 1999); however, partial eclipsing could be observed in a system with an inclination $`i`$ $``$ 60$`\mathrm{°}`$.

We have obtained followup spectroscopic and photometric observations of 2MASS 1237+6526 in order to test the binary hypothesis, as well as examine the possibility that this object was originally observed during a prolonged flare. In $`\mathrm{\S }`$2, we present red optical (6300–10100 Å) spectral data obtained over 1.6 yr using the Keck 10m Low Resolution Imaging Spectrograph (Oke et al., 1995, hereafter LRIS), and discuss the behavior of the H$`\alpha `$ emission line over this period. J-band monitoring observations obtained over two nights using the Palomar 60” Near-Infrared Camera (Murphy et al., 1995, hereafter IRCam) are presented in $`\mathrm{\S }`$3. We discuss how the imaging observations constrain the binary scenario in $`\mathrm{\S }`$4. Finally, in $`\mathrm{\S }`$5 we discuss how these follow-up observations constrain various emission mechanisms for this unique object.

## 2 Spectroscopic Observations

2MASS 1237+6526 has been observed on seven separate occasions in the red optical regime starting in 1999 July. Early data, obtained on 1999 July 16-18 (UT) using LRIS, are described in Burgasser et al. (2000a). Here, we present our latest LRIS observations obtained on 2000 March 5 and 2001 February 20 (UT). Data acquisition and reduction procedures for the first run are described below, and are similar to those used for the 2001 February data as described in Kirkpatrick et al. (2002).

### 2.1 Data Acquisition and Reduction

Conditions during our 2000 March 5 (UT) observations of 2MASS 1237+6526 were clear with 0$`\stackrel{}{\mathrm{.}}`$8–0$`\stackrel{}{\mathrm{.}}`$9 seeing. Because of its optical faintness (R $``$ 26; Matthews et al. 1996, Burgasser et al. 1999), the target was invisible on the guider camera, so acquisition into a 1$`\mathrm{}`$ slit was done by blind offset from a nearby optical source using positions determined from 2MASS coordinates. Spectral data were obtained using the 400 lines mm-1 grating blazed at 8500 Å, covering the spectral range 6300 to 10100 Å at 9 Å resolution. The OG570 order blocking filter was used to reduce second-order light. Two exposures of 1800s each were obtained; however, the data from the first exposure were fainter than expected, so the target was repositioned prior to the second exposure. Observations of the DZA4 white dwarf WD 1225-079 (Green, Schmidt, & Liebert, 1986; Kilkenny, 1986) and the B1 V flux standard Hiltner 600 (Hiltner, 1956; Hamuy et al., 1994) were also obtained for calibration using the same instrumental configuration.

Science data were reduced using standard IRAF222Image Reduction and Analysis Facility (IRAF) is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. routines. A series of one-second dark exposures were median-combined and subtracted from the science data to remove detector bias. Median-combined quartz-lamp flat-field exposures (reflected off of the interior dome) were used to normalize the response of the detector. Spectra were then extracted using the APEXTRACT routine. Curvature of the dispersion line was determined from bright standard star exposures and applied to the extraction of the 2MASS 1237+6526 spectrum. Wavelength calibration was done using Ne-Ar arc lamp exposures taken immediately after the target observations. We then computed a telluric correction from the white dwarf spectrum by interpolating over telluric O2 (6850–6900 Å A-band, 7600–7700 Å B-band) and H2O (7150–7300, 8100–8350, and 9300–9650 Å) bands, and we applied this correction to our flux standard and 2MASS 1237+6526 spectra333Despite the presence of Ca II HK absorption in this object (Kilkenny, 1986; Liebert, Wehrse, & Green, 1987), the triplet lines at 8498, 8542, and 8662 Å were not present and hence did not affect the telluric correction.. Finally, flux calibration was done by correcting our Hiltner 600 observations to the spectrophotometric data given by Hamuy et al. (1994), being careful to interpolate over telluric features and the Balmer H$`\alpha `$ line. This correction was then applied to the spectrum of 2MASS 1237+6526.

Reduced spectral data from our 2000 March and 2001 February observations are shown in Figure 1. Full spectra with flux plotted on a logarithmic scale are shown in the right panel, while a close up of the H$`\alpha `$ line on a linear scale is shown on the left. A combined spectrum from our 1999 July 16 (UT) observations is also plotted for comparison (Burgasser et al., 2000a). The 2000 March spectrum is roughly 2.6 times fainter than our 1999 July observations, likely the result of poor targeting onto the slit due to the blind offset and proper motion of 2MASS 1237+6526 (first observed by 2MASS on 1999 March 13 UT). The 2001 February spectrum, on the other hand, is brighter by a factor of 1.3; atmospheric transparency on this night was excellent, and the target was visible on the facility guide camera, enabling proper slit placement. Features present in the spectra are noted, including the strong 3($`\nu _1`$,$`\nu _3`$) and 2($`\nu `$1,$`\nu `$3)+2$`\nu `$2 H2O bands starting at 9250 and 9540 Å (Auman, 1967), a possibly weak 4($`\nu _1`$,$`\nu _3`$) CH4 band centered at 8950 Å (Dick & Fink, 1977), and lines of Cs I at 8521 and 8943 Å ($`6s^2S_{1/2}6p^2P_{1/2,3/2}`$). The dominant feature in these spectra is the strongly pressure-broadened K I resonance doublet ($`4s^2S_{1/2}4p^2P_{1/2,3/2}`$), centered at 7665 and 7699 Å, but extending over most of the spectral regime shown. Liebert et al. (2000) have shown that the red wing of the K I doublet is the primary contributor to the steep slope seen from 8000–10000 Å in the spectrum of SDSSp J162414.37+002915.6 (Strauss et al., 1999), and this feature is likely responsible for shaping the red optical spectrum of 2MASS 1237+6526 as well.

### 2.2 H$`\alpha `$ Emission in 2MASS 1237+6526

H$`\alpha `$ emission at 6563 Å is clearly seen in all of the spectra of Figure 1, and its presence in our recent data suggests that emission is persistent in 2MASS 1237+6526 over a period of at least 1.6 years. We have investigated the temporal behavior of H$`\alpha `$ emission by measuring the line flux in all of our spectra to date. Flux data were integrated within 20 Å of the emission peak after subtracting off the neighboring pseudo-continuum, the scatter of which was used to estimate uncertainties. Corrections have also been computed to compensate for slit losses, by scaling the individual spectra to match data from 1999 July 16 (UT) in the 8000–9500 Å region. The initial and corrected line fluxes are listed in Table 1.

Figure 2 plots these line fluxes versus UT time. All of the corrected values are in excellent agreement with the mean line flux from our 1999 July observations, (7.0$`\pm `$0.5)$`\times `$10-17 erg cm-2 s-1, with the exception of data taken on 1999 July 17 (UT). As reported by Burgasser et al. (2000a), these data were of extremely low signal-to-noise (again, likely due to poor slit placement), and it is possible that the formal errors are too low. Excepting this data point, we find that the H$`\alpha `$ emission in 2MASS 1237+6526 is very stable, firmly ruling out any kind of flaring event. The weighted mean value over all observations is $`f_{H\alpha }`$ = (7.44$`\pm `$0.08)$`\times `$10-17 erg cm-2 s-1, consistent with previous results.

### 2.3 Comparison to Other Cool Stars and Brown Dwarfs

To gauge how the strength of H$`\alpha `$ emission in 2MASS 1237+6526 compares to other active, cool dwarfs, we computed the relative H$`\alpha `$ luminosity for each spectrum, $`\mathrm{log}_{10}`$(L/Lbol) = $`\mathrm{log}_{10}`$($`f_{H\alpha }`$/$`f_{bol}`$), where

$$\mathrm{log}_{10}f_{bol}=0.4m_{bol}4.61$$
(1)
in units of erg cm-2 s-1 (Drilling & Landolt, 2000), and $`m_{bol}`$ = $`m_J`$ \+ BCJ is the apparent bolometric magnitude. Using a revised J-band magnitude for 2MASS 1237+6526 from the 2MASS Second Incremental Release Catalog and BCJ = 2.09$`\pm `$0.10 (Leggett et al., 1999, determined from measurements of Gliese 229B444Gliese 229B and 2MASS 1237+6526 have the same spectral type, T6.5 V (Burgasser et al., 2002a); hence, differences in the bolometric correction are likely to be negligible.), we derive $`m_{bol}`$ = 18.12$`\pm `$0.14 and a mean $`\mathrm{log}_{10}`$(L/Lbol) = $``$4.3, consistent with the results of Burgasser et al. (2000a).

We compare this value to the emission ratios of late-M, L, and T dwarfs in Figure 3. Using spectral data from Kirkpatrick et al. (1999, 2000, 2002) and Burgasser et al. (2000a, 2002b), we computed H$`\alpha `$ line fluxes directly from the spectra as described above, with bolometric fluxes for late-M and L dwarfs computed from 2MASS J-band photometry and BCJ values from Reid et al. (2001). For T dwarfs, we used 2MASS J-band photometry and BCJ = 2.09, yielding an upper limit to the true bolometric correction for the early- and mid-T dwarfs (Burgasser et al., 2002a), and hence a lower limit to the bolometric flux; i.e., overestimating the relative H$`\alpha `$ luminosity. Detections are plotted as open symbols, while upper limits (determined from pseudocontinuum noise levels) are indicated by downward arrows. For comparison, we have included M and L dwarf samples from Hawley, Gizis, & Reid (1996, crosses) and Gizis et al. (2000, diamonds). The trend of decreasing activity beyond M5–M7 V reported by Gizis et al. (2000) is readily apparent in this plot, and appears to continue well into the T dwarf regime. Indeed, only two objects have reliable detections beyond spectral type L5 V, the T2 V SDSS 1254-0122, with $`\mathrm{log}_{10}`$(L/Lbol) = $``$5.7, and 2MASS 1237+6526. The former object is over 25 times less luminous in H$`\alpha `$ then 2MASS 1237+6526, and upper limits for the other T dwarfs do not exceed $`\mathrm{log}_{10}`$(L/Lbol) = $``$5.2.

Burgasser et al. (2000a) compares 2MASS 1237+6526 to another unusually active cool dwarf, the M9.5 Ve PC 0025+0047 (Schneider et al., 1991), which is indicated in Figure 3 by the solid square. This object, with $`\mathrm{log}_{10}`$(L/Lbol) = $``$3.4, is also more than ten times brighter in H$`\alpha `$ than other objects in its class, although it is not extremely active in comparison to earlier emission M dwarfs (Hawley, Gizis, & Reid, 1996). PC 0025+0047 also shows long-term stability in its emission (Schneider et al., 1991; Mould et al., 1994; Martín, Basri, & Zapatero Osorio, 1999). These common characteristics provide further evidence that emission processes in these two objects may be related.

### 2.4 A Possible Velocity Shift in the 2000 March Spectrum?

Examination of the spectra of Figure 1 also shows an apparent blueshift in the H$`\alpha `$ line in our 2000 March data as compared to the other epochs. We measured a line center deviation of 3.2$`\pm `$0.5 Å in this spectrum with respect to the 2001 February data, where the uncertainty is derived from the scatter of telluric OH emission line wavelengths between these two data sets. This deviation corresponds to a velocity of 150$`\pm `$25 km s-1, which cannot be due to differences in the Earth’s rotational and orbital motion between the two epochs, amounting to only 2.4 km s-1.

As with our initial detection of H$`\alpha `$, we carefully scrutinized the reality of this subtle (1–2 pixel) line shift. Figure 4 shows the H$`\alpha `$ region in the raw image data from 2000 March and 2001 February. Telluric OH emission lines in the 6563 Å region appear as vertical bands, while the bright spots are H$`\alpha `$ emission in 2MASS 1237+6526; note the lack of continuum in this region. There is a slight offset in the emission between these two epochs with respect to the telluric lines, confirming the spectral results in Figure 1. However, this shift could be due to poor centering of the object on the slit. Both individual exposures taken on 2000 March 5 (UT) show a similar offset, despite repositioning of the object. The second exposure was brighter and probably better positioned, but was still 2.6 times fainter than our 1999 July results; hence, it too may have been misaligned. On the other hand, seeing on the nights of 2000 March 5 and 2001 February 20 (UT) was roughly equivalent to the width of the slit employed, so significant shifting of the spectrum due to poor slit placement is less likely. We were unable to perform an independent check using the weak 8521 and 8943 Å Cs I lines due to the low signal-to-noise of both spectra. We therefore conclude that there is some evidence for a velocity shift in the H$`\alpha `$ line which must be confirmed with higher resolution data.

Assuming its reality, the H$`\alpha `$ line shift (which appears to be asymmetric toward the blue) could arise from a single hot spot on the surface of a rapidly rotating brown dwarf. The velocity of the line is comparable to but less than the break-up speed for a substellar object, roughly 330 (520) km s-1 for a 30 (75) MJup brown dwarf. In this case, however, the integration time of each individual observation made (typically 20–30 minutes) would be a significant fraction of the rotation period (no less than about 50 minutes), and one would expect the line to be broadened rather than shifted. A more reasonable source would be a hot spot in orbit around the brown dwarf, possibly arising from an optically thin accretion disk. Discrete regions in an edge-on Keplerian disk would have line-of-sight radial velocities as high as 150 km s-1 at separations of roughly 5 (12) RJup for a 30 (75) MJup central source. The hot spot would have to be small, due to the lack of blue continuum seen in at shorter wavelengths (Burgasser et al., 2000a). It is worth noting that the magnitude and direction of the velocity shift is similar to that seen in some solar flaring events during an impulsive evaporation phase (Antonucci et al., 1985); however, we do not see the typically brighter redshifted condensation component (Fisher, Canfield, & McClymont, 1985).

## 3 Imaging Observations

In addition to our spectroscopic data, we searched for photometric variability in 2MASS 1237+6526, for the purpose of detecting an eclipse by a closely-separated companion.

### 3.1 Data Acquisition and Reduction

J-band monitoring observations were obtained over two nights, 2000 May 19–20 (UT), using the Palomar 60” IRCam. Conditions on both nights were generally hazy but clear, with average seeing of 1$`\stackrel{}{\mathrm{.}}`$2 to 1$`\stackrel{}{\mathrm{.}}`$8; some scattered cirrus was present during the observations of May 20. Data were acquired in sets of five 45s exposures dithered by 10$`\mathrm{}`$ on May 19 and 15$`\mathrm{}`$ on May 20 in a square pattern. A total of 155 and 165 images were obtained, spanning periods of 2.50 and 2.57 hours on May 19 and 20, respectively (Table 2).

Science data were initially pairwise differenced to eliminate bias and sky background, and then divided by a normalized flat-field image, constructed by differencing median-combined sets of lamp-on and lamp-off dome observations made on May 19. Bad pixels were identified in the flat-field image and corrected in the science images by interpolation. The four separate readouts of the IRCam NICMOS-3 HgCdTe array occasionally introduced additive quadrant biases that we removed by subtracting off the median from each quadrant. Finally, individual images in each set were shifted and coadded to produce a final reduced image (Figure 5).

### 3.2 Differential Photometry and Analysis of the Time Series

For each image, we computed aperture photometry ($`f_i`$) for 2MASS 1237+6526 and other sources in the field. Flux within a 2$`\stackrel{}{\mathrm{.}}`$5 (4.04 pixel) radius aperture was integrated after subtracting the mean background measured in a 6 to 12$`\mathrm{}`$ annulus around each object. Uncertainty estimates of the photometry, $`\delta f_i`$, included contributions by shot noise in the source and variations in the background (essentially flat field uncertainties). Instrumental magnitudes, $`m_i=2.5\mathrm{log}_{10}f_i`$, were then computed with corresponding uncertainties, $`\sigma _i=1.086(\frac{\delta f_i}{f_i})`$. The scatter in instrumental magnitudes for 2MASS 1237+6526 was roughly 0.03 mag on May 19 and 0.01 mag on May 20, with formal photometric uncertainties of $`\pm `$0.025 mag.

To compute a relative magnitude light curve, we first examined sources in the field for variability, using the procedure outlined in Bailer-Jones & Mundt (2001). We first determined relative magnitudes and uncertainties for each source:

$$\widehat{m}_i^{(j)}=m_i^{(k)}m_i^{ref}$$
(2)
$$\widehat{\sigma }_i^{(j)^2}=\sigma _{i}^{(j)}{}_{}{}^{2}+\sigma _{i}^{ref}{}_{}{}^{2},$$
(3)
where the reference magnitudes and uncertainties were computed as

$$m_i^{ref}=2.5\mathrm{log}_{10}\left(\underset{kj}{\overset{K}{}}f_i^{(k)}\right)$$
(4)
$$\sigma _{i}^{ref}{}_{}{}^{2}=\underset{kj}{\overset{K}{}}\sigma _{i}^{(k)}{}_{}{}^{2},$$
(5)
summing over the $`K`$ other sources in the field. After subtracting off the mean relative magnitude, we computed the reduced $`\chi ^2`$ for each source to determine if photometric variations were significant:

$$\chi _r^2=\frac{1}{n1}\underset{i}{\overset{n}{}}(\frac{\widehat{m}_i\widehat{m}_i}{\widehat{\sigma }_i})^2;$$
(6)
here, $`n`$ is the total number of images. We constrained calibrator sources to have $`\chi _r^2`$ $``$ 0.5, corresponding to a probability of variability $`P`$ $``$ 1% for $`n`$ = 30. Two calibrators present in images from both nights passed this constraint, 2MASSI J1237367+652600 (Calibrator 1) and 2MASSI J1237372+652536 (Calibrator 2). These sources have 2MASS J magnitudes of 15.07$`\pm `$0.05 and 15.34$`\pm `$0.06, respectively, and are indicated on Figure 5.

Relative magnitudes were then computed for 2MASS 1237+6526 using these two calibrator stars. The resulting time series is shown for both nights in Figure 6. The mean relative magnitude has been subtracted off in both datasets. Time series for the differenced calibrator magnitudes ($`m_i^{(1)}m_i^{(2)}`$) are also shown (offset) for comparison. Scatter in the relative magnitudes of 2MASS 1237+6526 was 0.020 and 0.014 mag on May 19 and 20, respectively, less than the formal photometric uncertainties. Corresponding $`\chi _r^2`$ values are 0.505 and 0.314, which are greater than those of the calibrator stars (0.271 and 0.237, respectively), but nevertheless reject variability (at levels above our detection threshold) at the $`>`$ 99% confidence level.

To examine to possibility of a modulated signal (e.g., multiple eclipsing events or geometric modulation), we constructed periodograms for each night’s data using the formalism of Lomb (1976) and Scargle (1982) as summarized in Press et al. (1992). Given the duration of the observations, we are sensitive to periods ranging from roughly 0.1–1.3 hr on each night. The computed periodograms (oversampled by a factor of 10) for both 2MASS 1237+6526 (black lines) and calibrator stars (grey lines) are shown in Figure 7. May 19 and 20 data are identified by solid and dashed lines, respectively. None of the peaks in the 2MASS 1237+6526 periodograms are significant, as they do not exceed the 50% false alarm probability level and are roughly equal in magnitude to insignificant peaks seen in the calibrator periodograms.

Variability over two nights also appears to be absent. There is a 0.024 mag offset in the mean relative magnitudes between the two nights; this is greater than the difference between the mean calibrator values (0.003 mag), but comparable to the photometric uncertainties. We compute $`\chi _r^2`$ = 0.593 over both nights, again rejecting variability at the $`>`$ 99% confidence level. Because of the large temporal gap in the time series data, periodicity analysis requires more sophisticated algorithms (e.g., the CLEAN algorithm; Roberts, Léhar, & Dreher 1987); however, given the lack of statistical evidence for variability, further analysis appears to be unwarranted.

Given the null results from these statistical tests, we conclude that there is no short-term variability in 2MASS 1237+6526 at the $`\pm `$0.025 mag level over periods of 0.1–2.5 hours and no evidence of periodicity for periods of up to 1.3 hours. There also appears to be no compelling indication of variability at this level for periods of up to $``$ 14 hours.

## 4 The Binary Hypothesis Revisited

In order to determine what constraints our null variability results impose on the interacting binary hypothesis, we have updated the basic theory described in Burgasser et al. (2000a) by comparison with theoretical models from Burrows et al. (1997). Our observables are (1) a period limit from the duration of monitoring observations, (2) a magnitude limit at which we are sensitive to variability, and (3) the properties of the primary. We therefore seek to derive a relation $`\mathrm{\Delta }`$J = $`f(p,i;M_1)`$, where $`p`$ and $`i`$ are the period and inclination of the system, and $`M_1`$ the mass of the primary. For this discussion, we normalize masses to MJup and distances to RJup, and assume a circular orbit and negligible angular momentum loss (e.g., to gravitational radiation or disk interaction). The period (in hours) is determined from the standard Keplerian relation:

$$p^2=4\pi ^2\frac{a^3}{G(M_1+M_2)}=8.79\frac{a^3}{M_1(1+q)}$$
(7)
in Jupiter-based units. In the case of an interacting binary, the secondary fills its Roche lobe, and hence its radius $`R_2`$ is determined by

$$R_2=R_L=a\frac{0.49q^{\frac{2}{3}}}{0.6q^{\frac{2}{3}}+\mathrm{ln}(1+q^{\frac{1}{3}})}.$$
(8)
(Eggleton, 1983). Brown dwarf radii, however, are fairly constant from 5–70 MJup (Burrows & Liebert, 1993), and are well-described by the relation (Zapolsky & Salpeter, 1969; Stevenson, 1991):

$$R=\alpha _1\frac{\pi M^{\frac{1}{3}}}{(1+1.8M^{\frac{1}{2}})^{\frac{4}{3}}},$$
(9)
where $`\alpha _1`$ is an age-dependent scale factor (essentially a correction to partial degeneracy) chosen to match the interior models of Burrows et al. (1997). Because $`R_2`$ depends only on the primary mass and mass ratio $`q`$, both $`a`$ and $`p`$ are functions of M1 and $`q`$ exclusively. Therefore, these two parameters define the dynamics of the system.

To determine the observable change in J-band magnitude caused by a nearly edge-on eclipse, we assumed a simple geometry of overlapping spheres along the line of sight, as diagrammed in Figure 8. In general, because dlnR/dlnM $`\frac{1}{3}<`$ 0 for brown dwarfs, the lower-mass secondary has a larger radius ($`R_2>R_1`$) for $`M_2>`$ 5 MJup. The relative area of the primary, $`A`$, obscured by the secondary during maximal eclipse is:

$$\pi A=\theta _10.5\mathrm{sin}2\theta _1+(\frac{R_2}{R_1})^2(\theta _20.5\mathrm{sin}\theta _2)$$
(10)
where

$$R_1\mathrm{sin}\theta _1=R_2\mathrm{sin}\theta _2$$
(11)
$$R_1\mathrm{cos}\theta _1+R_2cos\theta _2=a\mathrm{cos}i,$$
(12)
with $`\theta _1`$ and $`\theta _2`$ defined in Figure 8. These equations allow us to uniquely determine $`A`$ in terms of $`M_1`$, $`q`$, and $`i`$. Note that no eclipse occurs when $`\mathrm{cos}i>(R_1+R_2)/a`$, while a complete eclipse occurs for $`\mathrm{cos}i<(R_2R_1)/a`$.

The observed magnitude difference at maximal eclipse is then

$$\mathrm{\Delta }J\begin{array}{c}=2.5log_{10}(\frac{F_1(1A)+F_2}{F_1+F_2})\hfill \\ =2.5log_{10}(\frac{A}{1+\beta }),\hfill \end{array}$$
(13)
where $`\beta =F_2/F_1`$ is the flux ratio of the secondary to primary at J-band. We can determine the change in J-band magnitude as a function of Teff empirically from the two known companion T dwarfs, Gliese 229B (Nakajima et al., 1995) and Gliese 570D (Burgasser et al., 2000b). The former object is a T6.5 V with MJ = 15.52$`\pm `$0.06 (Leggett et al., 1999) and Teff $``$ 950 K (Allard et al., 1996; Marley et al., 1996; Tsuji et al., 1996); the latter object is a T8 V with MJ = 16.47$`\pm `$0.07 (Burgasser et al., 2000b) and Teff $``$ 800 K (Geballe et al., 2001). These values yield the very approximate relation dJ/dT = $``$0.0063 mag/K, which may only be accurate over a limited Teff range555This rough relation also ignores the effects of gravity and metallicity, which could be substantial.. Interior models from Burrows et al. (1997) show a power law relation between Teff and mass which is roughly dependent on age ($`\tau `$). A fit to the models yields the relations:

$$\frac{dT}{dM}=\{\begin{array}{cc}88.7M^{0.43}\hfill & \tau =0.5Gyr\hfill \\ 73.4M^{0.44}\hfill & \tau =1.0Gyr\hfill \\ 45.4M^{0.44}\hfill & \tau =5.0Gyr\hfill \end{array}K/M_{Jup}.$$
(14)
Combining Equation 14 with the expression for dJ/dT yields the secondary flux ratio:

$$2.5log_{10}\beta \begin{array}{c}=\frac{dJ}{dT}\frac{dT}{dM}\mathrm{\Delta }M\hfill \\ =\frac{dJ}{dT}\frac{dT}{dM}M_1(q1)\hfill \\ =\alpha _2\times (1q)\hfill \end{array}$$
(15)
where $`\alpha _2`$ is a constant dependent on age and primary mass. We have assumed that 2MASS 1237+6526, a T6.5 V similar to Gliese 229B, has a Teff $``$ 950 K, and masses 23, 30, and 55 MJup for $`\tau `$ = 0.5, 1.0, and 5.0 Gyr, respectively (Burrows et al., 1997). The various parameters are summarized in Table 3.

Using the above equations, we have calculated the expected magnitude drop ($`\mathrm{\Delta }`$Jecl) at minimum light for an eclipsing, interacting brown dwarf binary with Tprimary = 950 K, 0.5 $`\tau `$ 5.0 Gyr, 0 $`q`$ 1, and 40$`\mathrm{°}i`$ 90$`\mathrm{°}`$. Results are shown in Figure 9 for $`\tau `$ = 0.5 (solid line), 1.0 (dot-dashed line), and 5.0 (dashed line) Gyr, and (from top to bottom) $`i`$ = 90$`\mathrm{°}`$, 80$`\mathrm{°}`$, 70$`\mathrm{°}`$, and 60$`\mathrm{°}`$. Only values of $`\mathrm{\Delta }`$Jecl for which $`q<`$ 0.63 and M2 $`>`$ 5 MJup are displayed, corresponding to the limits for sustained mass loss (Burgasser et al., 2000a). Significant drops in the observed flux are predicted by this model, particularly for $`i`$ $``$ 80$`\mathrm{°}`$, where fluxes at maximum eclipse over 0.7 mag fainter than the nominal brightness are possible over the full range of periods shown. These flux deviations are due primarily to the substantial brightness ratio between the secondary and primary, which is particularly important at high inclination angles where obscuration is most complete. At lower inclination angles, the limiting factor is the orbital separation, which is larger for lower mass secondaries. No eclipse variations are predicted to occur for $`i50\mathrm{°}`$ over the range of ages examined. Note that significant systematic errors may be present for the lowest mass ratios (i.e., low secondary Teff) due to the rough estimate of $`\beta `$ in Equation 15; however, our observations generally sample only the highest mass ratios due to the period cutoff.

The predicted $`\mathrm{\Delta }`$Jecl from these models are generally much larger than our observational limits (indicated by the hatched region in Figure 8). The constraints imposed by the non-detection of variability are summarized in Table 3 for $`i`$ $`>`$ 60$`\mathrm{°}`$; in this case, $`q`$ $`>`$ 0.70, 0.53, and 0.37, corresponding to secondary masses M2 $`>`$ 16, 16, and 20 MJup for $`\tau `$ = 0.5, 1.0, and 5.0 Gyr, respectively. Note that for the 0.5 Gyr case, our mass ratio limit exceeds the minimum $`q`$ required to sustain mass loss. The lowest inclination angle eclipses that could have been detected by our observations for $`\tau `$ = 1.0 and 5.0 Gyr are $`i_{min}`$ = 54$`\mathrm{°}`$ and 56$`\mathrm{°}`$, with corresponding minimum secondary masses of 20 and 34 MJup, respectively.

## 5 Discussion

Based on the results of the monitoring observations, we can confidently rule out an interacting binary system for 2MASS 1237+6526 for $`i`$ $`>`$ 60$`\mathrm{°}`$ and M2 $`>`$ 20 MJup. We cannot rule out a less inclined system, although such a geometry would fail to explain the apparent line shift seen in our 2000 March spectral data. A compromise solution would be a moderately inclined system with a massive primary and large mass ratio. Conclusive synthesis of the spectroscopic and photometric data, however, requires verification of H$`\alpha `$ line shift.

Other emission mechanisms should also be considered. 2MASS 1237+6526 could be a young, and hence very low mass, brown dwarf that is still accreting material from a circum(sub)stellar disk. If this object were in a weak-lined T-Tauri (WLTT) phase, its H$`\alpha `$ line could arise from a small, optically thin accretion boundary layer, consistent with our spectral and photometric observations. Assuming a typical WLTT age (30–100 Myr; Hartmann, Kenyon, & Hartigan 1993) and Teff $``$ 950 K, the evolutionary models of Burrows et al. (1997) predict a mass of 3–12 MJup, below the Deuterium-burning limit. Note that 2MASS 1237+6526 is not associated with any known star-forming region, although it could have been ejected from its nascent cluster (Reipurth & Clarke, 2001). The absence of the 1.25 $`\mathrm{\mu m}`$ K I doublet in this object (I. McLean, priv. comm.) marginally supports the possibility of it being a low-gravity source, although we cannot rule out a purely temperature effect (Burgasser et al., 2002a).

Our follow-up observations of 2MASS 1237+6526 have enabled us to rule out flaring as a viable means of emission, and have placed significant constraints on the geometry and membership of an interacting binary system. Parallax (to determine absolute brightness) and space motion measurements are clearly required to further characterize this object. Furthermore, the apparent velocity shift described in $`\mathrm{\S }`$2.4, which may prove to be a vital clue to the origin of emission in 2MASS 1237+6526, requires confirmation and more detailed, higher-resolution investigation. The nature of this enigmatic brown dwarf may hopefully be revealed with further observations.

A. J. B. would like to thank Telescope Operators Terry Stickel (Keck) and Barrett “Skip” Staples (Palomar) for their support during the observations presented here. A. J. B. also acknowledges useful discussions with C. A. L. Bailer-Jones, D. Koerner, M. Marley, I. McLean, and M. R. Zapatero Osorio during the preparation of the manuscript, and valuable comments from our anonymous referee. A. J. B., J. D. K., and J. E. G. acknowledge the support of the Jet Propulsion Laboratory, California Institute of Technology, which is operated under contract with the National Aeronautics and Space Administration. A. J. B. also acknowledges support provided by NASA through a Hubble fellowship grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. J. L. acknowledges support from NASA JPL grant 961040NSF. 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 University of Massachusetts and the Infrared Processing and Analysis Center, funded by the National Aeronautics and Space Administration and the National Science Foundation.