text:
# Electronic properties of the dimerized one-dimensional Hubbard model using lattice density-functional theory

## I Introduction

Hohenberg and Kohn (HK) replaced the wave function by the electronic density $`\rho (\stackrel{}{r})`$ as the fundamental variable of the many-body problem and thereby achieved a crucial breakthrough in the theoretical description of the electronic properties of matter.hk  Since then, density functional theory (DFT) has been the subject of remarkable developments. Formal improvements, extensions, and uncountable successful applications to a large variety of problems have made of this theory the most efficient, albeit not infallible, method of determining physical and chemical properties of matter from first principles.parr-book ; gross-book  Density functional (DF) calculations are usually based on the Kohn-Sham (KS) scheme which reduces the correlated $`N`$-particle problem to the solution of a set of self-consistent single-particle equations.ks  While this transformation is formally exact, the form of the interaction-energy functional $`W[\rho (\stackrel{}{r})]`$ involved in the KS equations is not known explicitly. Practical implementations of DFT always require approximations to $`W[\rho (\stackrel{}{r})]`$, or equivalently to the exchange and correlation (XC) functional $`E_{\mathrm{XC}}[\rho (\stackrel{}{r})]`$, on which the quality of the results depends crucially. Therefore, understanding the functional dependence of $`W[\rho (\stackrel{}{r})]`$ and exploring new ways of improving its approximations are central to the development of DF methods.

The most extensively used forms for $`W[\rho (\stackrel{}{r})]`$ are presently the local density approximation (LDA), ks  its spin-polarized version or local spin-density approximation (LSDA)hedin  and the gradient corrected extensions,gc  which were originally derived from exact results for the homogeneous electron gas. Despite an unparalleled success in the most diverse areas, the LDA-based approach fails systematically in accounting for phenomena that involve strong electron-correlation effects as observed, for example, in Mott insulators, heavy-fermion materials, or high-$`T_c`$ superconductors. These systems are usually described in the framework of parametrized lattice Hamiltonians such as Anderson,andmod  Hubbard,hub  Pariser-Parr-Pople,ppp  and related models which focus on the most relevant electron dynamics at low energies. However, even with simplified model interactions and a minimal number of orbitals per atom, a detailed understanding of the electronic properties in the strongly correlated limit remains a serious theoretical challenge. Exact results are rare or numerically very demanding and a variety elaborate many-body techniques have been specifically developed in order to study this problem.fulde  Being in principle an exact theory, the limitations of the DF approach have to be ascribed to the approximations used for exchange and correlation and not to the underlying formalism. It is therefore very interesting to extend the range of applicability of DFT to the many-body lattice models that describe the physics of strongly correlated systems. Moreover, the development of lattice density-functional theory (LDFT) constitutes an intrinsically inhomogeneous approach and provides a true alternative to the LSDA and related gradient-corrected methods. Thus, studies on simple models can open new insights into the properties of $`W`$ that should also be useful for future extensions to more realistic Hamiltonians and first principles calculations.

Several physical problems have been already investigated by applying the concepts of DFT to lattice models, for example, the band-gap problem in semiconductors,gunn  the role of off-diagonal elements of the density matrix and the non-interacting $`v`$ representability in strongly correlated systems,godby  or the development of energy functionals of the density matrix with applications to Hubbard and Anderson models.carl  In previous worksprbxcfun ; ldftdim  we have considered a density-matrix functional theory of many-body lattice models, that is analogous to Gilbert’s approach in the continuum,gilb  and applied it to the Hubbard Hamiltonian with uniform nearest-neighbor (NN) hopping integrals $`t_{ij}=t`$. The interaction energy $`W`$ of the Hubbard model has been calculated exactly as a function of the density matrix $`\gamma _{ij}`$ for various periodic lattices having $`\gamma _{ij}=\gamma _{12}`$ for all nearest neighbors $`i`$ and $`j`$. An analysis of the functional dependence of $`W(\gamma _{12})`$ for different band fillings and lattice structures revealed very interesting scaling properties.prbxcfun  On this basis, a simple general approximation to $`W(\gamma _{12})`$ has been derived which yields a remarkable agreement with available exact results in 1D systems and which predicts successfully the ground-state energy and charge-excitation gap of the 2D Hubbard model in the complete range of interaction strength.ldftdim  This shows that DFT with an appropriate approximation to $`W`$ is an efficient tool for determining the electronic properties of many-body lattice models.

The purpose of this paper is to extend the method by allowing for alternations of the density-matrix elements $`\gamma _{ij}`$ between nearest neighbors in order to study the dimerized 1D Hubbard model. This problem has motivated a considerable research activity in past years, particularly concerning the role of electron correlations in the dimerization of polymer chains like polyacetylene.teo-poly  In this context two qualitatively different regimes may be distinguished depending on the relative importance of the intra-atomic Coulomb repulsion $`U`$ and the NN hopping integral $`t`$. On the one side, for small $`U/t`$, the dimerization can be regarded as a bond-order wave that opens a gap at the Fermi surface of the one-dimensional (1D) single-particle band-structure (Peierls distortion). On the other side, for large $`U/t`$, local charge fluctuations are severely reduced and the low-energy properties are dominated by antiferromagnetic (AF) correlations between spin degrees of freedom. In this case the dimerization can be regarded as an alternation of the strength of AF correlations along the chain or spin-Peierls state. One of our aims is to analyze the differences between these two types of behaviors in the framework of LDFT. The properties of dimerized chains are also very interesting from a purely methodological point of view. They provide in fact a simple, physically motivated means of exploring the functional dependence of $`W[\gamma _{ij}]`$ by including additional degrees of freedom, thereby allowing for a larger flexibility. Moreover, several exact results are available to compare with (e.g., Bethe-Ansatz solution for the non-dimerized Hubbard chain, finite-ring Lanczos diagonalizations, or density-matrix renormalization-group calculations) which allow to quantify the accuracy of the final results.

The remainder of the paper is organized as follows. In Sec. II the main steps in the formulation of LDFT are briefly recalled. The properties of the interaction-energy functional $`W`$ of the dimerized Hubbard model are discussed in Sec. III. The domain of representability of $`\gamma _{ij}`$ and the scaling behavior of $`W`$ are investigated. A simple explicit approximation to $`W`$ is derived, which is appropriate for direct calculations. Sec. IV is concerned with applications. The ground-state energy and the charge-excitation gap of finite and infinite dimerized chains are determined as a function of Coulomb-repulsion strength $`U/t`$ and hopping-integral alternation $`\delta t`$. The LDFT results are contrasted with accurate numerical solutions in order to discuss goals and limitations of the present approach. Finally, Sec. V summarizes the main conclusions and points out some future perspectives.

## II Density functional theory on a lattice

We consider the many-body Hamiltonian

$$H=\underset{ij\sigma }{}t_{ij}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }+\frac{1}{2}\underset{\genfrac{}{}{0pt}{}{klmn}{\sigma \sigma ^{}}}{}V_{klmn}\widehat{c}_{k\sigma }^{}\widehat{c}_{m\sigma ^{}}^{}\widehat{c}_{n\sigma ^{}}\widehat{c}_{l\sigma },$$
(1)
where $`\widehat{c}_{i\sigma }^{}`$ ($`\widehat{c}_{i\sigma }`$) is the usual creation (annihilation) operator for an electron with spin $`\sigma `$ at site or orbital $`i`$. The hopping integrals $`t_{ij}`$ define the lattice (e.g., 1D chains, square or triangular 2D lattices) and the range of single-particle interactions (e.g., up to first or second neighbors). From the ab initio perspective $`t_{ij}`$ is given by the external potential $`V_{\mathrm{ext}}(\stackrel{}{r})`$ and by the choice of the basis. $`V_{klnm}`$ defines the type of many-body interactions which may be repulsive (Coulomb like) or attractive (in order to simulate electronic pairing) and which are usually approximated as short ranged (e.g., intra-atomic). Eq. (1) is mainly used in this section to present the general formulation which can then be applied to various specific models by simplifying the interactions. A particularly relevant example, to be considered in some detail in following sections, is the single-band Hubbard Hamiltonian with NN hoppings.hub  The non-dimerized form of this model is obtained from Eq. (1) by setting $`t_{ij}=t`$ for $`i`$ and $`j`$ NN’s, $`t_{ij}=0`$ otherwise, and $`V_{klnm}=U\delta _{kl}\delta _{nm}\delta _{kn}`$.

The hopping matrix $`t_{ij}`$ plays the role given in conventional DFT to the external potential $`V_{ext}(\stackrel{}{r})`$. Consequently, the single-particle density matrix $`\gamma _{ij}`$ between lattice sites replaces the continuum density $`\rho (\stackrel{}{r})`$ as basic variable. The situation is similar to the density-matrix functional theory proposed by Gilbert for the study of non-local pseudo-potentials $`V_{\mathrm{ext}}(\stackrel{}{r},\stackrel{}{r}^{})`$.gilb ; donnelly ; val  The ground-state energy $`E_{gs}`$ and density matrix $`\gamma _{ij}^{gs}`$ are determined by minimizing the energy functional

$$E[\gamma _{ij}]=E_K[\gamma _{ij}]+W[\gamma _{ij}]$$
(2)
with respect to $`\gamma _{ij}`$. $`E[\gamma _{ij}]`$ is defined for all density matrices that can be written as

$$\gamma _{ij}=\underset{\sigma }{}\gamma _{ij\sigma }=\underset{\sigma }{}\mathrm{\Psi }|\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }$$
(3)
for all $`i`$ and $`j`$, where $`|\mathrm{\Psi }`$ is an $`N`$-particle state. In other words, $`\gamma _{ij}`$ derives from a physical state and is said to be pure-state $`N`$-representable.foot\_ens  The first term in Eq. (2) is given by

$$E_K=\underset{ij}{}t_{ij}\gamma _{ij}.$$
(4)
It represents the kinetic energy associated with the electronic motion in the lattice and includes all single-particle contributions. Notice that Eq. (4) yields the exact kinetic energy and that no corrections on $`E_K`$ have to be included in other parts of the energy functional as in the KS approach. The second term in Eq. (2) is the interaction-energy functional given bylevy

$$W[\gamma _{ij}]=min\left[\frac{1}{2}\underset{\genfrac{}{}{0pt}{}{klmn}{\sigma \sigma ^{}}}{}V_{klmn}\mathrm{\Psi }[\gamma _{ij}]|\widehat{c}_{k\sigma }^{}\widehat{c}_{m\sigma ^{}}^{}\widehat{c}_{n\sigma ^{}}\widehat{c}_{l\sigma }|\mathrm{\Psi }[\gamma _{ij}]\right],$$
(5)
where the minimization implies a search over all $`N`$-particles states $`|\mathrm{\Psi }[\gamma _{ij}]`$ that satisfy

$$\mathrm{\Psi }[\gamma _{ij}]|\underset{\sigma }{}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }[\gamma _{ij}]=\gamma _{ij}$$
(6)
for all $`i`$ and $`j`$. $`W[\gamma _{ij}]`$ represents the minimum value of the interaction energy compatible with a given density matrix $`\gamma _{ij}`$. It is often expressed in terms of the Hartree-Fock energy

$$E_{\mathrm{HF}}[\gamma _{ij}]=\frac{1}{2}\underset{\genfrac{}{}{0pt}{}{ijkl}{\sigma \sigma ^{}}}{}V_{ijkl}\left(\gamma _{ij\sigma }\gamma _{kl\sigma ^{}}\delta _{\sigma \sigma ^{}}\gamma _{il\sigma }\gamma _{kj\sigma }\right)$$
(7)
and the correlation energy $`E_\mathrm{C}[\gamma _{ij}]`$ as

$$W[\gamma _{ij}]=E_{\mathrm{HF}}[\gamma _{ij}]+E_\mathrm{C}[\gamma _{ij}].$$
(8)
$`W`$ and $`E_\mathrm{C}`$ are universal functionals of $`\gamma _{ij}`$ in the sense that they are independent of $`t_{ij}`$, i.e., of the system under study. They depend on the considered interactions or model, as defined by $`V_{klmn}`$, on the number of electrons $`N_e`$, and on the structure of the many-body Hilbert space, as given by $`N_e`$ and the number of orbitals or sites $`N_a`$.

$`E[\gamma ]`$ is minimized by expressing $`\gamma _{ij}=\gamma _{ij}+\gamma _{ij}`$ in terms of the eigenvalues $`\eta _{k\sigma }`$ (occupation numbers) and eigenvectors $`u_{ik\sigma }`$ (natural orbitals) as

$$\gamma _{ij\sigma }=\underset{k}{}u_{ik\sigma }\eta _{k\sigma }u_{jk\sigma }^{}.$$
(9)
Lagrange multipliers $`\mu `$ and $`\lambda _{k\sigma }`$ ($`\epsilon _{k\sigma }=\lambda _{k\sigma }/\eta _{k\sigma }`$) are introduced in order to impose the usual constraints $`_{k\sigma }\eta _{k\sigma }=N_e`$ and $`_i|u_{ik\sigma }|^2=1`$. Derivation with respect to $`u_{jk\sigma }^{}`$ and $`\eta _{k\sigma }`$ ($`0\eta _{k\sigma }1`$) yields the eigenvalue equationsldftdim ; gilb

$$\eta _{k\sigma }\underset{i}{}\left(t_{ij}+\frac{W}{\gamma _{ij}}\right)u_{ik\sigma }=\epsilon _{k\sigma }u_{jk\sigma }$$
(10)
with the following conditions relating $`\eta _{k\sigma }`$ and $`\epsilon _{k\sigma }`$:

$`\epsilon _{k\sigma }`$ $`<`$ $`\mu \mathrm{if}\eta _{k\sigma }=1,`$ (11)
$`\epsilon _{k\sigma }`$ $`=`$ $`\mu \mathrm{if}0<\eta _{k\sigma }<1,`$ (12)
and

$`\epsilon _{k\sigma }`$ $`>`$ $`\mu \mathrm{if}\eta _{k\sigma }=0.`$ (13)
Self-consistency is implied by the dependence of $`W/\gamma _{ij}`$ on $`\eta _{k\sigma }`$ and $`u_{ik\sigma }`$. Eqs. (10)–(13) hold exactly in all interaction regimes. They are analogous to well-known results of density-matrix functional theory in the continuum.gilb  However, notice the difference with the KS-like approach considered in Ref. godby , which assumes non-interacting $`v`$-representability, and where only integer occupations are allowed.

The importance of fractional natural-orbital occupations has already been stressed in previous density-matrix functional studies in the continuum.gilb  In fact, in the case of models one observes that $`0<\eta _{k\sigma }<1`$ for all $`k\sigma `$ except in very special situations like the uncorrelated limit ($`V_{klmn}=0`$) or the fully-polarized ferromagnetic state in the Hubbard model. This can be understood from perturbation-theory arguments —none of the $`\eta _{k\sigma }`$ should be a good quantum number for $`V_{klmn}0`$— and has been explicitly demonstrated in exact solutions for finite systems or for the 1D Hubbard chain.lieb-wu  Therefore, the case (12) is the only relevant one in general. All $`\epsilon _{k\sigma }`$ in Eq. (10) must be degenerate or equivalently

$$t_{ij}+W/\gamma _{ij\sigma }=\delta _{ij}\mu .$$
(14)
Clearly, approximations of $`W`$ in terms of the diagonal $`\gamma _{ii}`$ alone can never yield such a behavior. Given a self-consistent scheme that implements the variational principle, the challenge remains to find good approximations to $`W[\gamma _{ij}]`$ that are simple enough to be applied in practical calculations.

## III Interaction-energy functional

In order to determine $`W[\gamma _{ij}]`$ from Eq. (5) we look for the extremes of

$`F`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{klmn}{\sigma \sigma ^{}}}{}}\left[V_{klmn}\mathrm{\Psi }|\widehat{c}_{k\sigma }^{}\widehat{c}_{m\sigma ^{}}^{}\widehat{c}_{n\sigma }\widehat{c}_{l\sigma }|\mathrm{\Psi }\right]+\epsilon \left(1\mathrm{\Psi }|\mathrm{\Psi }\right)`$ (15)
$`+`$ $`{\displaystyle \underset{i,j}{}}\lambda _{ij}\left(\mathrm{\Psi }|{\displaystyle \underset{\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }\gamma _{ij}\right)`$
with respect to $`|\mathrm{\Psi }`$. Lagrange multipliers $`\epsilon `$ and $`\lambda _{ij}`$ have been introduced to enforce the normalization of $`|\mathrm{\Psi }`$ and the representability of $`\gamma _{ij}`$ as required by Eq. (6). Derivation with respect to $`\mathrm{\Psi }|`$, $`\epsilon `$ and $`\lambda _{ij}`$ yields the eigenvalue equationsprbxcfun

$$\underset{ij\sigma }{}\lambda _{ij}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }+\underset{\genfrac{}{}{0pt}{}{klmn}{\sigma \sigma ^{}}}{}V_{klmn}\widehat{c}_{k\sigma }^{}\widehat{c}_{m\sigma ^{}}^{}\widehat{c}_{n\sigma }\widehat{c}_{l\sigma }|\mathrm{\Psi }=\epsilon |\mathrm{\Psi },$$
(16)
and the auxiliary conditions $`\mathrm{\Psi }|\mathrm{\Psi }=1`$ and $`\gamma _{ij}=\mathrm{\Psi }|_\sigma \widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }`$. The Lagrange multipliers $`\lambda _{ij}`$ play the role of hopping integrals to be chosen in order that $`|\mathrm{\Psi }`$ yields the given $`\gamma _{ij}`$. The pure-state representability of $`\gamma _{ij}`$ ensures that there is always a solution. The subset of $`\gamma _{ij}`$ that can be represented by a ground-state of Eq. (16) for some $`\lambda _{ij}`$ is the physically relevant one, since it necessarily includes the absolute minimum $`\gamma _{ij}^{gs}`$ of $`E[\gamma _{ij}]`$. Nevertheless, it should be noted that pure-state representable $`\gamma _{ij}`$ may be considered that can only be described by excited states or by linear combinations of eigenstates of Eq. (16).prbxcfun

The general functional $`W[\gamma _{ij}]`$, valid for all lattice structures and for all types of hybridizations, can be simplified at the expense of universality if the hopping integrals are short ranged. For example, if only NN hoppings are considered, the kinetic energy $`E_K`$ is independent of the density-matrix elements between sites that are not NN’s. Therefore, the constrained search in Eq. (5) may be restricted to the $`|\mathrm{\Psi }[\gamma _{ij}]`$ that satisfy Eq. (6) only for $`i=j`$ and for NN $`ij`$. This reduces significantly the number of variables in $`W[\gamma _{ij}]`$ and renders the determination and interpretation of the functional dependence far simpler.

In Sec. III.1 we present and discuss exact results for the interaction energy $`W[\gamma _{ij}]`$ of the dimerized Hubbard model on representative finite and infinite chains. These are obtained by solving Eq. (16) using accurate numerical methods. The dependence of the interaction energy on the alternating NN density-matrix elements $`\gamma _{12}`$ and $`\gamma _{23}`$ is analyzed. Scaling properties are identified within the domain of representability of $`\gamma _{ij}`$. On the basis of these results we propose, in Sec. III.2, a simple general approximation to $`W(\gamma _{12},\gamma _{23})`$ which is useful for practical applications. A first test on the accuracy of this approximation is also provided by comparison with available exact solutions.

### III.1 Exact calculated $`W[\gamma _{ij}]`$ of the dimerized Hubbard model

In the following we consider the dimerized 1D Hubbard model which in the usual notation is given byhub

$$H=\underset{ij\sigma }{}t_{ij}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }+U\underset{i}{}\widehat{n}_i\widehat{n}_i.$$
(17)
The NN hopping integrals $`t_{ij}`$ take two alternating values: $`t_{i,i+1}=t_{12}=t+\delta t`$ for $`i`$ odd and $`t_{i,i+1}=t_{23}=t\delta t`$ for $`i`$ even. The corresponding interaction-energy functional reads

$$W[\gamma _{ij}]=min\left[U\underset{l}{}\mathrm{\Psi }[\gamma _{ij}]|\widehat{n}_l\widehat{n}_l|\mathrm{\Psi }[\gamma _{ij}]\right],$$
(18)
where the minimization is performed with respect to all $`N`$-particle states $`|\mathrm{\Psi }[\gamma _{ij}]`$ satisfying $`\mathrm{\Psi }[\gamma _{ij}]|_\sigma \widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }[\gamma _{ij}]=\gamma _{ij}`$ for NN $`ij`$. For repulsive interactions $`W[\gamma _{ij}]`$ represents the minimum average number of double occupations corresponding to a given degree of electron delocalization, i.e., to a given $`\gamma _{ij}`$. Eq. (16) then reduces to

$$\underset{\genfrac{}{}{0pt}{}{ij}{\sigma }}{}\lambda _{ij}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }+U\underset{i}{}\widehat{n}_i\widehat{n}_i|\mathrm{\Psi }=\epsilon |\mathrm{\Psi }.$$
(19)
This eigenvalue problem can be solved numerically for finite systems with various boundary conditions. To this aim we expand $`|\mathrm{\Psi }[\gamma _{ij}]`$ in a complete set of basis states $`|\mathrm{\Phi }_m`$ which have definite occupation numbers $`\nu _{i\sigma }^m`$ at all orbitals $`i\sigma `$: $`\widehat{n}_{i\sigma }|\mathrm{\Phi }_m=\nu _{i\sigma }^m|\mathrm{\Phi }_m`$ with $`\nu _{i\sigma }^m=0`$ or $`1`$. The values of $`\nu _{i\sigma }^m`$ satisfy the usual conservation of the number of electrons $`N_e=N_e+N_e`$ and of the $`z`$ component of the total spin $`S_z=(N_eN_e)/2`$, where $`N_{e\sigma }=_i\nu _{i\sigma }^m`$. For not too large clusters, the state $`|\mathrm{\Psi }_0[\gamma _{ij}]`$ corresponding to the minimum in Eq. (18) —the ground state of Eq. (19)— can be determined by sparse-matrix diagonalization procedures such as the Lanczos iterative method.lanczos  For large chains, the properties of $`|\mathrm{\Psi }_0[\gamma _{ij}]`$ can be calculated using the density-matrix renormalization-group (DMRG) methodwhite  which allows reliable extrapolations to the infinite-length limit. Finally, in the absence of dimerization ($`\delta t=0`$) translational symmetry implies that all NN $`\gamma _{ij}`$ are the same, and therefore one may set $`\lambda _{ij}=\lambda `$ for all NN $`ij`$. The lowest eigenvalue of Eq. (19) can then be determined from Lieb and Wu’s exact solution of the 1D Hubbard model following the work by Shiba.lieb-wu

In Fig. 1 the interaction energy $`W`$ of dimerized Hubbard chains is shown in the form of constant-energy curves given by $`W(\gamma _{12},\gamma _{23})=\lambda E_{\mathrm{HF}}`$, where $`\gamma _{12}=\gamma _{i,i+1}`$ for $`i`$ odd and $`\gamma _{23}=\gamma _{i,i+1}`$ for $`i`$ even are the density-matrix elements or bond orders between NN’s. $`E_{\mathrm{HF}}=U/4`$ stands for the Hartree-Fock energy, and $`\lambda `$ is a constant ($`0\lambda 1`$). Results are presented for the $`N_a=12`$ site ring and for the infinite 1D chain which were obtained from Eq. (19) by using Lanczos-diagonalization and DMRG methods, respectively.lanczos ; white  Only positive $`\gamma _{12}`$ and $`\gamma _{23}`$ are considered since this is the relevant domain when all the hopping integrals have the same sign. In bipartite lattices, like open chains or rings with even $`N_a`$, the sign of the NN bond orders can be changed without altering $`W`$ by changing the sign of the local orbitals at one of the sublattices. Thus, $`W(\gamma _{12},\gamma _{23})=W(\gamma _{12},\gamma _{23})`$. Moreover, $`W(\gamma _{12},\gamma _{23})=W(\gamma _{23},\gamma _{12})`$ as even an odd sites can be interchanged by a simple translation ($`N_a`$ even for rings).

The domain of definition of $`W`$ is restricted by the pure-state representability of $`\gamma _{ij}`$. The axes $`\gamma _{12}=0`$ and $`\gamma _{23}=0`$ in Fig. 1 represent a collection of disconnected dimers or fully dimerized states, while $`\gamma _{12}=\gamma _{23}`$ corresponds to non-dimerized states. In between, the degree of dimerization can be measured by the angle $`\varphi =\mathrm{arctan}(\gamma _{12}/\gamma _{23})`$. The degree of electron delocalization for each $`\varphi `$ is characterized by $`\gamma =\sqrt{\gamma _{12}^2+\gamma _{23}^2}`$, which is bounded by $`\gamma ^{\mathrm{}}(\varphi )\gamma \gamma ^0(\varphi )`$ in order that $`\gamma _{ij}`$ remains pure-state representable. The density-matrix elements along the curve $`\gamma =\gamma ^0(\varphi )`$ are the largest bond orders that can be achieved on a given lattice and for given $`N_a`$ and $`N_e`$ \[$`(\gamma _{12}^0,\gamma _{23}^0)=\gamma ^0(\varphi )(\mathrm{cos}\varphi ,\mathrm{sin}\varphi )`$\]. They represent the maximum electron delocalization for each $`\varphi `$ and yield the extremes of the kinetic energy $`E_K=t_{ij}\gamma _{ij}`$, with different $`\varphi `$ corresponding to different $`t_{12}/t_{23}`$. Thus, for $`\gamma =\gamma ^0(\varphi )`$ the density matrix can be represented by the ground state of the uncorrelated Hubbard model for some $`t_{12}/t_{23}`$ ($`U=0`$). In the absence of degeneracy the underlying electronic state $`|\mathrm{\Psi }_0`$ is a single Slater determinant and $`W(\gamma _{12}^0,\gamma _{23}^0)=E_{\mathrm{HF}}`$. Consequently, the upper bound for $`\gamma `$ coincides with the $`\lambda =1`$ curve in Fig. 1. The correlation energy $`E_C=WE_{\mathrm{HF}}`$ vanishes as expected in the fully delocalized limit. For $`U=0`$ the minimization of the energy $`E=E_K`$ as a function of $`\gamma _{ij}`$ can be stated in terms of the representability of $`\gamma _{ij}`$ alone. In this case the equilibrium condition yielding $`\gamma _{ij}^{gs}`$ is achieved at the borders of the domain of representability, more precisely, when the normal to the curve $`\gamma =\gamma ^0(\varphi )`$ is parallel to $`\stackrel{}{}E_K=(t_{12},t_{23})`$.

Concerning the lower bound $`\gamma ^{\mathrm{}}(\varphi )`$, one should first note that as $`\gamma `$ decreases, $`\gamma <\gamma ^0(\varphi )`$, it is possible to construct correlated states $`|\mathrm{\Psi }[\gamma _{ij}]`$ having increasingly localized electrons. Charge fluctuations can then be reduced more efficiently for smaller $`\gamma `$, and therefore the Coulomb interaction energy decreases with decreasing $`\gamma `$ \[see Eq. (18) and Fig. 1\]. $`W`$ reaches its minimum value $`W_{\mathrm{}}=Umax\{0,N_eN_a\}`$ in the strongly correlated limit where $`\gamma =\gamma ^{\mathrm{}}(\varphi )`$. For half-band filling this corresponds to a fully localized state having $`\gamma ^{\mathrm{}}(\varphi )=0`$ and $`W_{\mathrm{}}=0`$. However, note that for $`N_eN_a`$, $`W`$ reaches $`W_{\mathrm{}}`$ already for $`\gamma ^{\mathrm{}}(\varphi )>0`$ since partially delocalized states can be found having minimal Coulomb repulsion. This is the case for example in a fully-polarized ferromagnetic state.

Fig. 1 also provides a qualitative picture of the functional dependence of $`W`$ for dimerized chains. On the one side, for strongly dimerized states ($`\varphi 0`$ or $`\varphi \pi /2`$) the constant-$`W`$ curves resemble circumference arcs, the gradient $`\stackrel{}{}W`$ being approximately radial. This type of behavior is most clearly seen for weak or moderate correlations ($`\lambda 0.3`$), while in the localized regime ($`\lambda 0.1`$) it holds only for a very limited range of $`\varphi `$ around $`\varphi =0`$ or $`\varphi =\pi /2`$. On the other side, for weakly to moderately dimerized states ($`\varphi \pi /8`$$`\pi /4`$) the level curves can be regarded in first approximation as straight lines parallel to $`\gamma _{12}=\gamma _{23}`$. The very weak dependence of $`\stackrel{}{}W`$ on $`\varphi `$ implies that for $`\varphi \pi /4`$ the ground-state values of $`\gamma _{12}^{gs}`$ and $`\gamma _{23}^{gs}`$, which result from the minimization of $`E=E_K+W`$, are very sensitive to the hopping alternation $`\delta t`$. In fact, significant variations of $`\varphi `$ are necessary until $`\stackrel{}{}W=\stackrel{}{}E_K(1+\delta t/t,1\delta t/t)`$ even for $`\delta t/t1`$. This is particularly notable for weak correlations since the $`W=E_{\mathrm{HF}}`$ curve is strictly linear for $`|\varphi \pi /4|<0.05`$. Therefore, a discontinuous change from $`\gamma _{12}^{gs}/\gamma _{23}^{gs}=1`$ to $`\gamma _{12}^{gs}/\gamma _{23}^{gs}=0.91`$ is found at $`U=0`$ and arbitrary small $`\delta t`$. For $`U>0`$, $`\gamma _{12}^{gs}`$ and $`\gamma _{23}^{gs}`$ are continuous functions of $`\delta t`$, although the dependence on $`\delta t`$ remains very strong for small $`\delta t`$, as can be inferred from the level curves in the figure. Comparing subfigures (a) and (b) one observes that the results for $`N_a=12`$ and $`N_a=\mathrm{}`$ are quite similar. The rather rapid convergence with chain length suggests that $`W(\gamma _{12},\gamma _{23})`$ is not very sensitive to the details of the considered system, even if the minimization constraints in Eq. (18) apply only to NN bond orders. This is of interest for practical applications, as it will be discussed below.

In Fig. 2 the interaction energy $`W`$ is shown as a function of $`\gamma `$ for representative values of $`\varphi `$, including in particular the fully-dimerized ($`\varphi =0`$) and non-dimerized ($`\varphi =\pi /4`$) cases. Despite the quantitative differences among the various $`\varphi `$, several qualitative properties are shared by all the curves: (i) As already discussed, the domain of representability of $`\gamma `$ is bound for each $`\varphi `$ by the bond order $`\gamma ^0(\varphi )`$ in the uncorrelated limit. $`\gamma ^0`$ decreases monotonously with increasing $`\varphi `$ for $`0\varphi \pi /4`$ showing that a compromise between $`\gamma _{12}`$ and $`\gamma _{23}`$ is made when the two bonds are active. This is an important contribution to the $`\varphi `$-dependence of $`W`$. (ii) In the delocalized limit, $`W(\gamma ^0,\varphi )=E_{\mathrm{HF}}=U/4`$ for all $`\varphi `$, since the electronic state yielding the largest $`\gamma `$ is a single Slater determinant. Moreover, one observes that $`W/\gamma `$ diverges at $`\gamma =\gamma ^0`$ for all $`\varphi `$. This is a necessary condition in order that the ground-state density matrix satisfies $`\gamma ^{gs}<\gamma ^0`$ for arbitrary small $`U>0`$, as expected from perturbation theory. (iii) Starting from $`\gamma =\gamma ^0`$, $`W`$ decreases with decreasing $`\gamma `$ reaching its lowest possible value $`W=0`$ for $`\gamma =0`$ ($`N_e=N_a`$). The decrease of $`W`$ with decreasing $`\gamma `$ means that the reduction of the Coulomb energy due to correlations is done at the expense of kinetic energy or electron delocalization. (iv) In the limit of small $`\gamma `$ one observes that $`W\gamma ^2`$. Therefore, for $`U/t1`$, $`\gamma ^{gs}t/U`$ and $`E_{gs}t^2/U`$, a well known result in the Heisenberg limit of the half-filled Hubbard model.fulde

### III.2 Scaling approximation to $`W[\gamma _{ij}]`$

In order to compare the $`\gamma `$-dependence of $`W`$ for different $`\varphi `$ and to analyze its scaling behavior it is useful to bring the domains of representability for different $`\varphi `$ to a common range by considering $`W(\gamma ,\varphi )`$ as a function of $`\gamma /\gamma ^0(\varphi )`$, as displayed in Fig. 3. In this form the results for different $`\varphi `$ appear as remarkably similar, showing that the largest part of the dependence of $`W`$ on the ratio $`\gamma _{12}/\gamma _{23}`$ comes from the domain of representability of $`\gamma _{ij}`$ given by its upper bound $`\gamma ^0(\varphi )`$ \[$`\gamma ^{\mathrm{}}(\varphi )=0`$ for half-band filling\]. An analogous scaling behavior has been found in previous numerical studies of $`W(\gamma _{12})`$ of non-dimerized Hubbard models, where $`\gamma _{12}`$ refers to the NN density-matrix element.prbxcfun  In this case one observes that $`W(\gamma _{12})`$ depends weakly on system size $`N_a`$, band-filling $`n=N_e/N_a`$, and lattice structure, if $`W`$ is measured in units of the Hartree-Fock energy $`E_{\mathrm{HF}}`$ and if $`\gamma _{12}`$ is scaled within the relevant domain of representability $`[\gamma _{12}^{\mathrm{}},\gamma _{12}^0]`$. In the present context, Fig. 3 implies that the change in $`W`$ associated to a given change in the degree of delocalization $`\gamma /\gamma ^0(\varphi )`$ can be regarded as nearly independent of $`\varphi `$ and system size. A good general approximation to $`W(\gamma ,\varphi )`$ can then be obtained by applying such a scaling to the functional dependence extracted from a simple reference system. An appropriate choice is provided by the fully-dimerized chain corresponding to $`\varphi =0`$, which can be worked out analytically. In this case the system consists of a collection of dimers and the exact interaction energy reads

$$W(\gamma ,\varphi =0)=\frac{UN_a}{4}\left(1\sqrt{1\gamma ^2}\right).$$
(20)
Scaling the functional dependence of the dimer interaction energy to the $`\varphi `$-dependent domain of representability one obtains

$$W_0(\gamma ,\varphi )=\frac{UN_a}{4}\left(1\sqrt{1\left[\frac{\gamma }{\gamma ^0(\varphi )}\right]^2}\right),$$
(21)
which we propose as approximation to $`W`$ for dimerized systems. Notice that $`W_0(\gamma ,\varphi )`$ preserves the previous general properties (i)–(iv) and that it is of course exact for $`\varphi =0`$ \[$`\gamma ^0(\varphi =0)=1`$\]. In practice, the system specific function $`\gamma ^0(\varphi )`$ can be easily obtained by integration of the single-particle spectrum.

It is important to remark that the density matrices $`\gamma _{ij}`$ involved in the approximate functional $`W_0`$ are pure-state $`N`$-representable. Eq. (21) applies to the $`\gamma _{ij}`$ obtained by scaling the off-diagonal elements of the density matrices $`\gamma _{ij}^0`$ that derive from uncorrelated states $`|\mathrm{\Psi }_0`$ having $`N_e`$ electrons on $`N_a`$ sites, and a uniform density distribution $`\mathrm{\Psi }_0|_\sigma \widehat{n}_{i\sigma }|\mathrm{\Psi }_0=N_e/N_a=1`$. In other terms, $`\gamma _{ij}`$ has the form $`\gamma _{ij}=\lambda \gamma _{ij}^0`$ with $`0\lambda 1`$ for all $`ij`$, and $`\gamma _{ii}=\gamma _{ii}^0=1`$ for all $`i`$. In order to show the pure-state representability of $`\gamma _{ij}`$, we consider two normalized $`N`$-particles states $`|\mathrm{\Psi }_a`$ and $`|\mathrm{\Psi }_b`$ satisfying $`\mathrm{\Psi }_a|_\sigma \widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }_b=0`$ for all $`ij`$. This condition is fulfilled, for example, by states $`a`$ and $`b`$ having different defined total spins $`S`$ or $`S_z`$, or by superpositions of pure-$`S`$ or pure-$`S_z`$ states sharing no common eigenvalues. The density-matrix represented by $`|\mathrm{\Psi }=\alpha |\mathrm{\Psi }_a+\beta |\mathrm{\Psi }_b`$ with $`\alpha ^2+\beta ^2=1`$ is then given by $`\gamma _{ij}=\mathrm{\Psi }|_\sigma \widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }=\alpha ^2\gamma _{ij}^a+\beta ^2\gamma _{ij}^b`$, where $`\gamma _{ij}^a`$ and $`\gamma _{ij}^b`$ are the density matrices corresponding to $`|\mathrm{\Psi }_a`$ and $`|\mathrm{\Psi }_b`$. Therefore, all the density matrices in the segment defined by $`\gamma _{ij}^a`$ and $`\gamma _{ij}^b`$ are pure-state $`N`$-representable. The representability of a scaled uncorrelated $`\gamma _{ij}^0`$ at half-band filling follows from the previous lemma by taking $`|\mathrm{\Psi }_a=|\mathrm{\Psi }_0`$, which has $`S=0`$ or $`1/2`$, and $`|\mathrm{\Psi }_b`$ equal to the fully localized state with one electron per site and maximal $`S=N_e/2`$, for which $`\gamma _{ij}^b=0`$ for all $`ij`$, and $`\gamma _{ii}^b=1`$ for all $`i`$. Consequently, the $`\gamma _{ij}`$ in the domain of definition of $`W_0`$ and the ground-state density matrices $`\gamma _{ij}^{gs}`$ derived from it are all pure-state $`N`$-representable.

Fig. 4 compares Eq. (21) with the exact $`W(\gamma ,\varphi )`$ for a $`12`$-site Hubbard ring and for the infinite chain. One observes that the proposed approximation follows rather closely the exact results for all $`\gamma `$ and $`\varphi `$. The largest discrepancies are found for vanishing or moderate dimerization (e.g., $`\varphi =3\pi /16`$ or $`\varphi =\pi /4`$) and relatively large $`\gamma `$ ($`\gamma 0.8`$). In all cases the quantitative differences remain small ($`|W_0W|/U0.047`$ for $`\varphi =\pi /4`$ and $`|W_0W|/U0.045`$ for $`\varphi =3\pi /8`$) which is quite remarkable taking into account the simplicity of the approximation. In the following, Eq. (21) is applied in the framework of LDFT to determine several properties of the dimerized 1D Hubbard model. Comparison is made with exact results whenever possible in order to assess the performance of the method.

## IV Dimerized Hubbard chains

In Figs. 5 and 6 the ground-state energy $`E_{gs}`$, kinetic energy $`E_K`$, and Coulomb energy $`E_C`$ of the 1D Hubbard model are given as a function of the Coulomb repulsion strength $`U/t`$ for different hopping alternations $`\delta t`$. Accurate numerical results are also shown which were obtained using the Lanczos diagonalization methodlanczos  for $`N_a=12`$ or the DMRG methodwhite  for the infinite chain. In the case of the non-dimerized infinite chain, Lieb and Wu’s exact solutionlieb-wu  is taken as reference. The results for $`N_a=12`$ and $`N_a=\mathrm{}`$ are qualitative very similar. $`E_{gs}`$ increases monotonically with $`U/t`$ since $`E_{gs}/U=\widehat{n}_i\widehat{n}_i>0`$, vanishing in the limit of $`U/t=\mathrm{}`$. For $`U/t<4`$ this is essentially a consequence of the increase of $`E_CU`$, as $`E_K`$ and $`\gamma _{ij}`$ remain very much like in the uncorrelated $`U=0`$ state. In contrast, for $`U/t>4`$ the electrons become increasingly localized, and the increase of $`E_{gs}`$ results form the increase of $`E_K`$ which approaches zero as $`|\gamma _{ij}|`$ decreases. At the same time $`E_C`$ tends to zero as charge fluctuations are suppressed (see Figs. 5 and 6).

Comparison between LDFT and the exact results shows a very good agreement. This concerns not only $`E_{gs}`$ but also the separate kinetic and Coulomb contributions indicating that electron localization and intra-atomic correlations are correctly described for all $`U/t`$. Moreover, this also shows that the results obtained for the ground-state energy are not the consequence of strong compensations of errors. Concerning the accuracy of $`E_K`$ and $`E_C`$ one generally observes that a somewhat higher precision is achieved for $`E_K`$, which functional dependence is known exactly, as compared to $`E_C`$, which derives from an approximation to $`W`$ \[Eq. (21)\]. For $`\delta t/t0.1`$ the LDFT calculations are nearly indistinguishable from the exact ones (e.g., $`|E_{gs}E_{gs}^{ex}|/t0.03`$ for $`\delta t/t=0.1`$). Even the largest quantitative discrepancies, found for the non-dimerized chain at intermediate $`U/t`$, are pretty small. For instance, for $`\delta t=0`$ and $`U/t=4`$ we obtain $`|E_{gs}E_{gs}^{ex}|/t=0.020`$ for the 12-site ring and $`|E_{gs}E_{gs}^{ex}|/t=0.044`$ for the infinite chain. Comparing Figs. 5 and 6 one observes that the performance of the method is sometimes higher for the 12-site ring than for the infinite chain. For example, Fig. 6 shows that $`E_K`$ ($`E_C`$) is slightly overestimated (underestimated) for $`\delta t=0`$ and $`U/(U+4t)=0.7`$$`0.8`$, whereas for $`N_a=12`$ a much better agreement with the exact result is found (see Fig. 5). In any case it is important to recall that no artificial symmetry breaking is required to describe correlation-induced localization correctly, as it often occurs in other approaches (e.g., mean-field spin-density-wave state). Moreover, the present calculations remain simple and numerically not demanding, since the minimization of $`E[\gamma _{ij}]`$ is performed using analytical expressions for $`E_K`$ and $`W`$ \[see Eqs. (4) and (21)\]. One concludes that LDFT, combined with Eq. (21) as approximation to the interaction energy functional, provides an efficient and correct description of the ground-state properties of the 1D Hubbard model in the complete range of interaction strength and dimerization.

The charge-excitation or band gap

$$\mathrm{\Delta }E_c=E_{gs}(N_e+1)+E_{gs}(N_e1)2E_{gs}(N_e)$$
(22)
is a property of considerable interest in strongly correlated systems which can can be related to the discontinuity in the derivative of the kinetic and correlation energies per site with respect to the electron density $`n`$. The determination of $`\mathrm{\Delta }E_c`$ constitutes a much more serious challenge than the calculation of ground-state properties like $`E_{gs}`$, $`E_K`$, and $`E_C`$ particularly in the framework of a density-functional formalism. Results for $`\mathrm{\Delta }E_c`$ of the 1D Hubbard model are given in Figs. 7 and 8 as a function of the Coulomb repulsion strength $`U/t`$ for different values of the hopping alternation $`\delta t`$ ($`n=1`$). $`\mathrm{\Delta }E_c`$ vanishes for $`\delta t=0`$ and $`U/t=0`$, and increases with increasing $`U/t`$ or $`\delta t`$. Comparison between LDFT and Lanczos exact diagonalizations ($`N_a=12`$) or the Bethe-Ansatz solutionlieb-wu  ($`N_a=\mathrm{}`$ and $`\delta t=0`$) shows fairly small quantitative discrepancies. In the most difficult non-dimerized case we find $`|\mathrm{\Delta }E_c\mathrm{\Delta }E_c^{ex}|<0.18t`$ for $`N_a=12`$, and $`|\mathrm{\Delta }E_c\mathrm{\Delta }E_c^{ex}|<0.34t`$ for $`N_a=\mathrm{}`$. For small $`U/t`$ and $`\delta t=0`$, $`\mathrm{\Delta }E_c`$ is somewhat underestimated for $`N_a=12`$ and overestimated for $`N_a=\mathrm{}`$. The latter is mainly due to the fact that Eq. (21) fails to reproduce the exponential decrease of $`\mathrm{\Delta }E_c`$ for $`U/t0`$ ($`N_a=\mathrm{}`$ and $`\delta t=0`$).lieb-wu  As in previous properties the accuracy improves with increasing $`\delta t`$. Fig. 7 shows that the LDFT results for non-vanishing dimerization and $`N_a=12`$ are very close to the exact ones ($`|\mathrm{\Delta }E_c\mathrm{\Delta }E_c^{ex}|/t<0.011`$ already for $`\delta t/t=0.1`$). Therefore, one expects that the predictions for $`N_a=\mathrm{}`$ and $`\delta t>0`$ should be reliable. Finally, one may note that in the limit of large $`U/t`$, the hopping alternation $`\delta t`$ has little effect on the charge gap. As the electrons tend to localize for $`U/t\mathrm{}`$, $`\mathrm{\Delta }E_cU+E_b`$ where $`E_b=4t`$ is the energy of the bottom of the single-particle band. The present lattice-density-functional scheme describes correctly the crossover from a band insulator to a Mott insulator which occurs in dimerized chains as $`U/t`$ is varied from the weak-interaction to the strong-interaction regime.

## V Summary and Outlook

A novel density-functional approach to lattice-fermion models has been applied to the dimerized 1D Hubbard Hamiltonian. In this framework the basic variable is the single-particle density matrix $`\gamma _{ij}`$ and the key unknown is the interaction-energy functional $`W[\gamma _{ij}]`$. In the present paper we have first investigated the functional dependence of $`W`$ on the density-matrix elements $`\gamma _{12}`$ and $`\gamma _{23}`$ between nearest neighbors in dimerized chains ($`\gamma _{i,i+1}=\gamma _{12}`$ for $`i`$ odd and $`\gamma _{i,i+1}=\gamma _{23}`$ for $`i`$ even). Rigorous results for $`W(\gamma _{12},\gamma _{23})`$ were derived from finite-ring Lanczos diagonalizations and from DMRG calculations for the infinite chain. An analysis of these exact results shows that $`W`$ can be appropriately scaled as a function $`\gamma /\gamma ^0(\varphi )`$, where $`\gamma =\sqrt{\gamma _{12}^2+\gamma _{23}^2}`$ and $`\gamma ^0(\varphi )`$ is the largest representable $`\gamma `$ for a given $`\varphi =\mathrm{arctan}(\gamma _{12}/\gamma _{23})`$. A simple general approximation to $`W`$ was then proposed which takes advantage of this scaling behavior and which provides with a unified description of correlations from weak to strong coupling regimes. Finally, using this approximation, several ground-state properties and the charge-excitation gap of dimerized chains have been determined successfully as a function of Coulomb repulsion strength $`U/t`$ and hopping alternation $`\delta t`$.

The accuracy of the results encourages more or less straightforward applications of the present approach to related problems such as multi-leg ladders or the two-dimensional square lattice with first and second NN hoppings ($`t`$-$`t^{}`$ Hubbard model). Moreover, the possibility of generalizing the present scaling approximation to an arbitrary number of independent variables $`\gamma _{ij}`$ deserves to be investigated in detail, since it would open the way to applications in very low symmetry situations including metal clusters and disordered systems.

###### Acknowledgements.

Helpful discussions with Dr. Ph. Maurel are gratefully acknowledged. This work has been supported by CONACyT (Mexico) through the project W-8001 (Millennium initiative). Computer resources were provided by IDRIS (CNRS, France).


format: text

text:
# The Type Ia Supernova 2001V in NGC 398711footnote 1 Based on observations obtained at German-Spanish Astronomical Centre, Calar Alto, operated by the Max-Planck-Institute for Astronomy, Heidelberg, jointly with the Spanish National Commission for Astronomy

## 1 Introduction

SN 2001V was discovered spectroscopically by P. Berlind at the F. L. Whipple Observatory on Feb. 19, 2001 (Jha et al., (2001)). The supernova was immediately classified as of Type Ia based on the Si ii absorption trough at 6150 Å. The blue continuum suggested that this SN was discovered before maximum light. This was also strengthened by the high expansion velocity (14 000 kms-1) derived from the center of the Si ii line. Because the light curves of Type Ia SNe are now frequently used to infer distances to the host galaxies, such SNe that are identified before maximum are very important, because i) the light curve can be sampled more effectively, ii) the inferred distances are better constrained by the light variation around maximum than at later phases.

The host galaxy, NGC 3987 belongs to the small group centered on NGC 4005 (Nilson, (1973), Gregory & Thompson, (1978)). This group is located within the medium-sized cluster Zw 127-10 (Zwicky, et al., (1961)). NGC 3987 is an edge-on, Sbc-type galaxy with remarkable central dust lane. The IRAS point source 11547+2528 is located close to the optical center of this galaxy, which was also detected in radio (Condon & Broderick, (1986)).

SN 2001V is the first supernova discovered in this galaxy. The SN occured in the outskirts, almost at the visible edge of the host (Fig.1). The position of the SN as well as its blue colour at maximum may suggest that the effect of interstellar extinction is probably not very high for SN 2001V. Indeed, the galactic component of the reddening in this direction is only $`E(BV)=0.02`$ mag, according to Schlegel et al., (1998).

There are two recent distance estimates for NGC 3987: Mould et al, (1993) lists $`\mu _0=33.82`$ mag (58.1 Mpc) based on near-IR Tully-Fisher relation, while van Driel et al., (2001) gives $`v=4361`$ kms-1 as a group-averaged radial velocity corrected for Virgo-infall that results in $`d=67.1`$ Mpc (adopting $`H_0=65`$ kms-1Mpc-1). Based on the average of these distance estimates, the expected maximum brightness of a Type Ia SN is about 14.6 mag. According to the NASA/IPAC Extragalactic Database, the redshift of NGC 3987 is $`z=0.015`$. Thus, SN 2001V is a relatively nearby, low-redshift SN.

In the followings, details of the observations and data reductions are given, then the results of the photometric analysis are presented and discussed.

## 2 Observations and data reduction

The light variation of SN 2001V was followed on 18 nights, starting from $`t=8`$ days (with respect to $`B`$maximum), extending up to $`t=+62`$ days. The applied telescopes and detectors are listed in Table 1. All data were collected through standard Johnson-Cousins $`BVRI`$ filters.

The CCD-frames have been reduced in IRAF111IRAF is distributed by NOAO which is operated by the Association of Universities for Research in Astronomy (AURA) Inc. under cooperative agreement with the National Science Foundation. First, the instrumental magnitudes of the SN and the selected comparison stars (Fig.1) were derived with aperture photometry using the task digiphot/apphot. The radius of the aperture was 6 pixels (about 2 x FWHM, see Table 1), while the sky level was determined in a 5 pixel-wide annulus with inner radius of 10 pixels. Because the SN is located close to the “tip” of NGC 3987, most of the pixels in the annulus were not affected significantly by the light from the host galaxy. The sky level was determined by calculating the modal average (3 x median $``$ 2 x mean) of the intensities in the annulus. The calculations were done interactively, and the results were plotted on the screen and inspected visually in order to detect any obvious systematic errors. The average sky level around the SN was always very close to the background around the comparison stars, no clear systematic effect could be found.

Second, the whole dataset have been re-reduced using PSF-photometry (digiphot/daophot), as also advised by the referee of this paper. PSF-photometry is a superior photometric method if the background strongly varies around the objects, and its removal is complicated. Since SN 2001V is somewhat contaminated by its host, the use of PSF-photometry may be useful to separate the light of the SN from that of the galaxy. This method requires bright stars with high S/N to construct a good PSF. Because there are only a few of such stars in the field of NGC 3987, and also the field of view of most of the telescopes used in this study was quite small, a lot of frames contained only 4 - 5 stars around the SN. Nevertheless, the PSF of each frame was determined interactively. The analytic component of the PSF was approximated by the built-in penny2 function in daophot, but in most cases it was based on only 2-3 stars. The frames made with the Schmidt-telescope (Table 1) contained much more field objects, but the PSF on these frames showed significant positional dependence. Therefore, a second-order variable PSF model has been constructed for these pictures, while constant PSF was determined for all other frames. Individual sky levels were calculated for all objects, and the background was subtracted iteratively during the fitting of the PSF in daophot/allstar. Then, the residuals were examined visually on the subtracted frames. The stars as well as the SN were adequately removed from most frames, but slight residuals were present at the position of the brightest stars on some frames. This was probably caused by the uncertainty of the PSF due to the small number of PSF-stars.

Differential magnitudes of SN 2001V have been computed using the comparison stars labeled in Fig.1 (see also Table 2 below). After transforming them to the standard system (Vinkó et al., (2001)), the brightness of the SN was calculated from the calibrated BVRI magnitudes of each comparison star (see below). Then, the SN magnitudes belonging to the same frame were averaged.

In order to search for any systematic effect caused by the reduction procedure, we have compared the SN magnitudes from the aperture- and PSF-photometry. Fig.2 presents the aperture minus PSF-magnitudes in all filters as a function of time. It is apparent that most of the differences are within $`\pm 0.1`$ mag. Naturally, the differences are higher at later phases, when the SN became fainter, but there is no visible systematic trend in the data. We conclude that both the aperture- and PSF-photometry of SN 2001V presented here is affected by approximately the same amount of random errors at the $`\pm 0.1`$ mag level, mainly due to the technical limitations of the applied instruments (lower S/N, small field of view, very few PSF-stars). In order to reduce the uncertainties introduced by the reduction method, the SN magnitudes resulted from both aperture- and PSF-photometry were averaged, and these magnitudes were accepted as the final result.

At first, the magnitudes of the comparison stars were calibrated via Landolt standards, observed at Calar Alto Observatory with the 1.2 m Cassegrain, under photometric conditions on Aug.11, 2001. The reliability of this dataset have been checked by using the calibrated photometry of some of the field stars made by the CfA Supernova Group at the F. L. Whipple Observatory. Because the comparison stars used in this study, and those adopted by the CfA-group were different, we have re-computed the magnitudes of our comparison stars using the CfA-dataset as secondary standards on a set of $`BVRI`$ frames collected with the Schmidt telescope of Konkoly Observatory on March 14, 2001. After performing aperture photometry (with the same parameters described above) and standard transformation, the “Calar Alto” calibration and the “CfA/Konkoly” calibration was found to be fully consistent. The differences between the two datasets were within $`\pm 0.020.03`$ mag. Finally, the results from the CfA/Konkoly calibration was adopted, in order to keep maximal consistency with data of SN 2001V from other observatories. The final standard magnitudes of the comparison stars are collected in Table 2. The uncertainties reflect the statistical errors given by daophot and also the magnitude differences between the two calibrating datasets.

The calibrated standard magnitudes of SN 2001V are listed in Table 3. The estimated errors of each point are given in parentheses. These uncertainties were calculated as $`\sigma ^2=\sigma _p^2+\sigma _m^2+\sigma _t^2`$, where $`\sigma _p`$ is the statistical error given by daophot, $`\sigma _m`$ is the error due to the photometric method (assumed to be equal to the difference between the aperture- and PSF-photometry, plotted in Fig.2) and $`\sigma _t`$ is the error introduced by the standard transformation (adopted as 0.03 mag for each data).

The observations listed in Table 3 were compared with the $`R`$band CCD photometric data of K.Hornoch (Hornoch, (2001)) and found reasonable agreement within $`\pm 0.1`$ mag. More recently Mandel et al., (2002) published $`B_{\mathrm{max}}=14.64\pm 0.03`$ mag based on a more extensive photometric dataset. This agrees very well with our template fitting (see next section) resulting in $`B_{\mathrm{max}}=14.60\pm 0.16`$ (the larger uncertainty is due to the lack of data around maximum in our photometry). This agreement may give some support to our belief that the data in Table 3 probably do not suffer from large systematic errors.

## 3 Results

The BVRI light curves were analyzed by the Multi-Colour Light Curve Shape (MLCS) method of Riess et al., (1998) (the SN template vectors were kindly provided by A. Riess). The timescale of the observations was divided by $`1+z`$ = 1.015 to remove the effect of time dilation. The value of $`K`$-correction in the $`V`$ band at this redshift varies between $`0.01`$ and $`0.03`$ mag (Hamuy et al., (1993)), which is definitely smaller than the uncertainty of our measurements, thus, it was not taken into account. The maximum absolute magnitude of the fiducial light curve was adopted as $`M_V^{\mathrm{max}}=19.46`$ mag (Richardson et al., (2002)). The fitting was computed simultaneously to all light curves. Details of the procedure are described in Vinkó et al., (2001). The minimum of the $`\chi ^2`$ was found at the following parameters: $`T_0(B)`$ = JD 2451973.0 $`\pm `$ 0.2 (fiducial $`B`$-maximum), $`E(BV)`$ = $`0.05\pm 0.02`$ mag, $`\mu _0`$ = $`34.36\pm 0.14`$ mag and $`\mathrm{\Delta }`$ = $`0.47\pm 0.04`$ mag. The given uncertainties are formal errors corresponding to the geometry of the $`\chi ^2`$ function around minimum. The fitted light curves are shown in Fig.3. Note, that the $`T_0(B)`$ parameter given here corresponds to the maximum of the fiducial light curve (with $`\mathrm{\Delta }=0.0`$), the actual maximum of SN 2001V in $`B`$ occured one day later, at JD 2451974.

The parameters above imply some interesting properties of SN 2001V. Comparing the $`E(BV)`$ colour excess with the reddening due to Milky Way dust ($`E(BV)_{\mathrm{gal}}=0.02`$, Schlegel et al., (1998)), it is visible that the ISM in the host galaxy contributes to the dust absorption. Since NGC 3987 is an edge-on galaxy with strong central dust lane, it is not an unexpected result. The host-galaxy reddening could have been even more severe if the SN had occured closer to the central plane of NGC 3987, like e.g. in the case of SN 1986G (Phillips et al., (1987)), or, more recently SN 2002cv in NGC 3190 (Meikle & Mattila, (2002)).

The reddening was also estimated using the “Lira-Phillips law” (Phillips et al., (1999)) as an independent check of the value given by the MLCS method. This method uses the apparent homogeneity of the unreddened $`(BV)_0`$ colours of SNe Ia in the $`30<t_V<90`$ phase interval ($`t_V`$ is the phase measured from $`V`$maximum). Unfortunately, there is only one observed $`(BV)`$ data point in Table 3 that can be used for this purpose ($`(BV)=0.98\pm 0.10`$ mag). The moment of $`V`$maximum was estimated from the template light curves fitted by the MLCS method. It occured at $`T_0(V)`$ = 2451976 $`\pm 1`$ JD. After correcting for time dilation, the phase of the observed $`(BV)`$ point was found to be $`t_V=46.3\pm 1`$ days. The relation

$$(BV)_0=0.7250.0118(t_V60)$$

given by Phillips et al., (1999) resulted in $`(BV)_0=0.89\pm 0.02`$. This implies $`E(BV)=0.09`$ mag, which is not too far from the value given by MLCS. However, its uncertainty is probably higher, because this estimate is based on a single observation, and the “Lira-Phillips law” itself may have an intrinsic uncertainty at the $`0.1`$ mag level (Li et al., 2001a ). On the other hand, the similarity of the reddening values given by two independent methods may mean that the estimated $`E(BV)`$ is close to the real value.

Adopting $`E(BV)=0.05`$ mag, the total absorption in the $`B`$, $`V`$, $`R`$ and $`I`$ bands is $`A_B=0.21`$, $`A_V=0.16`$, $`A_R=0.13`$, $`A_I=0.09`$ mag, respectively, using the galactic reddening law given by Schlegel et al., (1998).

The luminosity parameter $`\mathrm{\Delta }=0.47`$ implies that SN 2001V is a super-luminous Type Ia SN, relative to the majority of such SNe. Super-luminous SNe are often referred as SN 1991T-, or SN 1999aa-like events (see Li et al., 2001b ). The distinction between different subtypes of SNe Ia is usually based on spectroscopic properties, e.g. the strength of the Si ii trough at 6150 Å. Overluminous SNe Ia have weaker Si ii 6150 Å line, and the Fe iii lines are strong. In addition, SN 1999aa-like events usually show strong Ca ii H & K lines, which are weak in SN 1991T-like SNe. These spectroscopic characteristics correlate with photometric properties: the light curves of SN 1991T-like SNe declines more slowly, and they are bluer in maximum than normal SNe Ia. Thus, overluminous SNe may be recognizable from photometry, e.g. via the MLCS-method.

We do not have spectroscopic data of SN 2001V at our disposal, but two spectra of SN 2001V, obtained around maximum, are available in graphical form from the CfA-website 222http://cfa-www.harvard.edu/cfa/oir/Research/supernova. As the referee of this paper, Dr. Weidong Li, pointed out, these spectra are very similar to that of SN 1999aa: broad, but less prominent Si ii 6150 Å line, strong Fe iii and Ca ii H & K lines. These spectral features fully support our conclusion (based on photometry), that SN 2001V is an overluminous SN Ia.

Note, however, that it is not clear whether the higher luminosity is necessarily connected with the spectroscopic peculiarities. There are examples, such as SN 1992bc (Riess et al., (1998)), that were super-luminous, but otherwise showed normal spectra around maximum. Thus, the overluminosity of some Type Ia SNe might be due to e.g. statistical fluctuation of the ejected 56Ni mass that powers the light curve, rather than a difference in the explosion mechanism or the physical state of the progenitor. Recently Richardson et al., (2002) found that the maximum brightness $`M_B`$ of SNe Ia has a Gaussian distribution with $`\sigma =0.56`$ mag (corrected for the effect of extinction).

From yet unpublished photometry, recently Mandel et al., (2002) concluded that SN 2001V was a “normal” Type Ia event, with light and colour curves very similar to SN 1990N. On the other hand, they also determined the initial decline rate, and found $`\mathrm{\Delta }m_{15}(B)=0.99\pm 0.05`$. This decline rate may suggest an overluminous, rather than a “normal” SN Ia (e.g. Hamuy et al., (1996)). Although our photometry has less accuracy and phase coverage, we also attempted to estimate this parameter and found $`\mathrm{\Delta }m_{15}(B)=0.9\pm 0.1`$. According to the correlation between $`\mathrm{\Delta }m_{15}(B)`$ and $`\mathrm{\Delta }`$ (e.g. Riess et al., (1998)), this is consistent with the $`\mathrm{\Delta }`$ parameter given by the MLCS-method. Furthermore, Fig.4 shows the essential similarity of the light curves of SN 2001V with that of SN 1991T ($`\mathrm{\Delta }m_{15}(B)=0.94`$, Hamuy et al., (1996)).

Thus, all available data consistently support the conclusion that SN 2001V was an intrinsically bright SN Ia. At maximum, it was brighter than the fiducial Type Ia SN by about 0.5 mag. Adopting the parameters given by the MLCS-method, the inferred absolute magnitudes of SN 2001V at maximum are

$`M_B^{\mathrm{max}}`$ = $`19.97\pm 0.15`$ mag

$`M_V^{\mathrm{max}}`$ = $`19.96\pm 0.15`$ mag

$`M_R^{\mathrm{max}}`$ = $`19.96\pm 0.15`$ mag

$`M_I^{\mathrm{max}}`$ = $`19.63\pm 0.15`$ mag

The true distance modulus, $`\mu _0=34.36\pm 0.14`$ mag, corresponds to $`74.5\pm 5`$ Mpc. The uncertainty given here is the formal error of the fitting, and does not include possible systematic effects, such as the zero-point of the SN distance scale. The MLCS distance marginally agrees with the kinematic distance of NGC 3987 corrected for Virgo-infall (67 Mpc, see Sect.1), the difference is about $`1.5\sigma `$. The Tully-Fisher distance (58 Mpc, Sect.1) is about 20 percent shorter. The TF distance modulus ($`\mu _0(TF)=33.82`$, Mould et al, (1993)) can be brought into agreement with the MLCS distance by adding 0.54 mag to the previous one, in accord with the suggestion by Shanks (Shanks, (1997)) for matching the TF and SN distance scales.

## 4 Summary

The results of this paper can be summarized as follows.

1. We presented BVRI photometry of SN 2001V starting from 8 days before maximum light and extending up to 60 days past maximum.
2. The reddening was estimated by applying the MLCS template fitting method for all light curves. This resulted in $`E(BV)=0.05\pm 0.02`$ mag, indicating that there was significant extinction in the host galaxy affecting the light from SN 2001V. The empirical “Lira-Phillips law” gave similar, but slightly higher colour excess ($`E(BV)=0.09`$ mag). Its uncertainty is probably higher, because it is based only a single observation. Therefore, the result from the MLCS-method was accepted as final.
3. The light curve analysis indicates that SN 2001V was an over-luminous SN by about 0.5 mag relative to the fiducial SN Type Ia. This is consistently supported by the overall similarity of the light curves and the $`\mathrm{\Delta }m_{15}(B)`$ decline rate with those of SN 1991T, and also the publicly available spectra of SN 2001V. The spectral features suggest that SN 2001V may be a SN 1999aa-like object (Li et al., 2001b ).
4. The SN distance was inferred from the MLCS method and found to be 74.5 $`\pm 5`$ Mpc. This marginally agrees with the previous kinematic distance, but about 20 percent longer than the Tully-Fisher distance of NGC 3987.

###### Acknowledgements.

This work was supported by Hungarian OTKA Grants No. T032258, T034615, the “Bolyai János” Research Scholarship to JV, the OM FKFP 0010/2001 grant from Hungarian Ministry of Education, Pro Renovanda Cultura Hungariae Foundation and MTA-CSIC Joint Project No.15/1998 from Hungarian and Spanish Academy of Sciences. We acknowledge the permission by Prof. A. Riess to use the MLCS template vectors. We thank the referee, Dr. Weidong Li, for a lot of critical comments and suggestions that helped us to improve the paper. We are grateful to Dr. T. Matheson and the Supernova Group led by Prof. R. Kirshner at the Harvard-Smithsonian Center for Astrophysics for providing their calibration data prior to publication. Thanks are also due to Dr. L.L.Kiss for his kind help during some observing sessions carried out at Szeged Observatory. The NASA Astrophysics Data System, the SIMBAD and NED databases and the Canadian Astronomy Data Centre were frequently used to access data and references. The availability of these services are gratefully acknowledged.


text:
# Pulsation at the tip of the first giant branch?

## 1 Introduction

After leaving the main sequence, stars of low to intermediate mass (0.8$`M_{}M8M_{}`$) become red giants and climb the first giant branch (RGB). At the tip of the RGB (TRGB), they experience the so-called helium flash, followed by a sudden decrease in luminosity. Then, the stars become red giants second time and ascend the asymptotic giant branch (AGB), increasing the scale of mass-loss and pulsation. On the upper AGB they become optically obscured because of the circumstellar dust shell formed by heavy mass-loss, making them accessible only in the infrared (Iben & Renzini 1983).

Recently, large and homogeneous NIR photometry data for the Large Magellanic Cloud (LMC) has become available owing to large-scale NIR surveys such as 2MASS (Skrutskie et al. 1997) and DENIS (Epchtein et al. 1997). However, as they are single-epoch surveys, they do not provide variability information. Variability information is necessary to get accurate mean magnitudes that are also indispensable for a proper understanding of the red giant evolution.

Many optically-visible variables were found in the LMC, as the by-products of the gravitational lensing experiments like MACHO (Wood 2000), MOA (Noda et al. 2002) and OGLE (Udalski et al. 1999). However, most, if not all, of high mass-losing variables might have been missed by these projects, because they are not readily detectable in the optical wave bands due to their surrounding circumstellar dust shells. Typical examples of such dust-enshrouded variables are OH/IR stars and dusty carbon stars (e.g., Habing 1996). As these dust-enshrouded variables are thought to be at the latest stage of AGB evolution (e.g., Habing 1996), the understanding of final evolution of low to intermediate mass star is greatly helped if their accurate quantities such as their number density and other physical parameters (mean magnitudes and periods etc.) would be obtained through deeper and larger survey.

Many low-amplitude variables were newly discovered by the aforementioned optical surveys. The vast majority are at the luminosities below the TRGB. Alves et al. (1998) suggested that they are AGB stars in the early evolutionary phase prior to entering the thermal-pulsing phase of the AGB. However, it is not clear yet if they really are or not.

In this letter, we present some first results of our NIR monitoring survey which was started in December 2000, using the InfraRed Survey Facility (IRSF) at Sutherland, South African Astronomical Observatory. The limiting magnitudes (10$`\sigma `$) of this survey are 17.4, 17.0 and 15.8 at $`J`$, $`H`$ and $`K`$, respectively. At these limits, we can cover all of the AGB stars except extremely red ones (i.e., $`JK5.0`$, such as those found in NGC419 and NGC1978 by Tanabé et al. 1998) and partway down the RGB populations in the LMC. We have been monitoring the total area of 3 square degrees along the LMC bar (see Figure. 1), sufficiently large to do statistical analysis and to make complete catalog of variable red giants.

## 2 Observations and reductions

The IRSF consists of a specially constructed 1.4m alt-azimuth telescope to which is attached a large-format 3-channel infrared camera, SIRIUS (Simultaneous three-colour InfraRed Imager for Unbiased Surveys). The SIRIUS can observe the sky in the three wave bands $`J`$(1.25$`\mu `$m), $`H`$(1.63$`\mu `$m) and $`K_s`$(2.14$`\mu `$m) simultaneously and has a field of view of 7.8 arcmin square with a scale of 0.45 arcsec/pixel. Details of the instrument are found in Nagashima et al. (1999) and Nagayama et al. (2002).

### 2.1 Photometry

Standard image reduction (i.e., dark subtraction, flat-field correction, sky subtraction) was made by using the SIRIUS pipeline (Nakajima, private communication). All images were reduced in the same way and a single image comprises 10 dithered 5 sec exposures. Photometry was performed on the dithered images with DoPHOT in the fixed position mode (Schechter, Mateo & Saha 1993).

At first, the instrumental magnitudes were transformed to the LCO system ones based on observations of a few dozen stars from Persson et al. (Persson et al. 1998). Then, the LCO system magnitudes were further transformed to the CIT/CTIO system ones through the equations in Persson et al. We estimate the uncertainties in the series of these system transformations should not exceed 0.05 mag. The photometry has a signal-to-noise ratio (S/N) of about 20 at $`J`$, $`H`$ and $`K`$ magnitudes of at 16.1, 16.1 and 15.0, respectively. Typical photometric errors are summarized in Table 1.

In the present study, the data are not dereddened. Instead, we show the direction and magnitude of the reddening vector in each diagram, based on the relations in Koornneef (1982).

### 2.2 Detections of variable stars

We used Image Subtraction Method (ISM) package ISIS.V2.1 (Alard & Lupton 1998; Alard 2000) to detect variable stars. This method can find variable stars even in very crowded fields such as LMC bar centre, and its high efficiency in detecting variables is well known from the previous/ongoing surveys (MACHO, e.g., Alcock et al. 2001; MOA, e.g., Bond et al. 2001; OGLE, e.g., Alard 1999; etc.). In this study, variable stars with light variations of 0.03 mag or more would all have been detected. So far, nearly ten thousand variables have been detected in our surveying area, and their $`JHK`$ light curves have been accumulated.

## 3 Results

In the survey area, 717,239 sources were detected by DoPHOT. Because of the scan overlaps, there should be many multiple entries of a single source in this original sample. We eliminated such multiple entries based on (1) spatial proximity ($`|\mathrm{\Delta }r|1.5^{^{\prime \prime }}`$) and (2) difference in photometry ($`|\mathrm{\Delta }m_{J,HandK}|0.1`$ mag). It leaves 562,429 sources. Further, requiring detections at all wave bands with S/N$``$20 in $`K`$ leaves 226,793 sources. Among the 226,793 sources, DoPHOT classified 184,677 sources as ‘perfect’ stars in all three wave bands (i.e., object type $`=`$ 1 in $`J`$, $`H`$ and $`K`$). Hereafter, we refer these 184,677 ‘perfect’ stars as the main sample.

Among the main sample, 5,188 stars were found to be variables. The full discussions of their pulsation properties, such as pulsation period, amplitude, variability type and mean magnitude etc. will be made in the future paper. In this letter, we will concentrate only on the information that a star is variable or not.

## 4 Discussion

In Figure 2 left panel, we show the colour-magnitude diagram ($`JK,K`$) and colour-colour diagram ($`HK,JH`$) of the main sample. In each diagram, numbers of stars in bins of area 0.025$`\times `$0.025 mag2 are computed and the fiducial colour is applied according to the number density of stars in each position (see the annotated colour wedge). Because the plotted magnitudes are single-epoch ones, these are affected by the scatter due to the light variation of each source. As one can see, the lower right area and the upper part ($`K9.5`$) of the colour-magnitude diagram are sparsely populated. These are attributed to the $`J`$ limiting magnitude and $`K`$ band saturation of the SIRIUS, respectively.

Nikolaev & Weinberg (2000) analyzed the 2MASS LMC colour-magnitude diagram and identified the major stellar populations. Compared to their colour-magnitude diagram, ours has a narrow and sharp shape, suggesting a smaller Galactic foreground content in our sample. As this indicates, we can reasonably ignore the contribution of Galactic foreground stars here, because our survey region is centred on the LMC bar (see Figure. 1), where foreground contribution is negligibly small compared to the number of the intrinsic LMC stars (Cioni et al. 2000a).

From the result of our variable star survey, now we know which star is a variable. Since neither 2MASS nor DENIS provide variability information, it will be interesting if we add it to the colour-magnitude and the colour-colour diagram. Figure 2 right panel is the same diagram as the left one, but it only contains the 5,188 variables extracted from the main sample. By comparing left and right panels, we can tell which part of the diagram the variable stars reside in.

### 4.1 Obscured variables

Stars in the dust-enshrouded phase are particularly interesting because of their importance in many areas, such as in the supply source of interstellar matter and as a key to answer long standing problems, like the carbon star mystery (Iben 1981). Because of their redness, most, if not all, of these red variables have been possibly missed by the previous optical surveys, such as MACHO, MOA and OGLE. Previous observations of dust-enshrouded variables in the LMC (e.g., Zijlstra et al. 1996; Reid 1991; Wood et al. 1992; van Loon et al. 1998) were biased toward bright ones such as IRAS sources, and were carried out only for a few tens of stars. Therefore, it will be invaluable if their pulsation period, amplitude, and mean magnitudes would be known by our survey. We have to defer the full discussion of the properties of the obscured variables pending further variability observations, but it is of interest to show some preliminary results.

The Midcourse Space Experiment (MSX) satellite surveyed the whole LMC in the mid-infrared and many dust-enshrouded stars were detected (Egan, van Dyk & Price 2001; Wood & Cohen 2001). These red stars are thought to be obscured by the dusty circumstellar envelope, which was developed by their mass-loss (Tanabé et al. 1997). In the main sample, there are 254 very red ($`JK>2.4`$) stars. Most of them might be included in the MSX sources analysed by Wood & Cohen (2001). By comparing the left and right panel of Figure. 2, one can easily notice the fact that most of such mass-losing stars are variables (194 stars out of 254 stars, 76.4%). This observational result corroborates the idea that the mass-loss and variability are connected.

Some stars from Hughes & Wood (1990) have also very red colour and most of them are Carbon-rich. However, their sample only extends up to $`JK3.0`$, and there is still a possibility that dusty Oxygen-rich variables like OH/IR stars (e.g., Wood et al. 1992; van Loon et al. 1998) may have redder colours. Many obscured variables (up to $`JK5.0`$) were found by our survey, but for most of them, the spectral type remains unknown. Infrared spectroscopic follow-up work would be interesting.

### 4.2 RGB variables?

In spite of the difficulty of distinguishing between the AGB and the RGB stars on the colour-magnitude diagram, variability is usually regarded as a sign of AGB stars. However, Origlia et al. (2002) recently found a significant scale of mass-loss among the first giant branch stars near the RGB tip. Since a strong mass-loss is believed to be associated with pulsation, it is natural to search for the light variation in the RGB stars.

In the upper panel of Figure 3, we show the differential stellar $`K`$ magnitude distribution $`N(m)`$ ($`dK=`$ 0.1 mag) of the main sample. The discontinuity in the $`N(m)`$ is clearly seen around $`K12.2`$, which corresponds to the TRGB. Most of the stars lying above the TRGB are intermediate age AGB stars, while both RGB and AGB stars are present just under the TRGB with the ascendancy of the RGBs. Renzini (1992) suggested that in a composite population the $`N(m)`$ must exhibit a drop by more than a factor of 4 at the luminosity corresponding to the TRGB, and the ratio of AGB to RGB stars near the tip of the RGB is about 1/3 for the intermediate age stars. Our result confirms the predominancy of the intermediate population in the LMC.

We show the $`K`$ magnitude distribution for the 5,769 variables in the lower panel of Figure 3. As is clearly seen in it, there are two peaks in the luminosity of variables; one consists of genuine AGB variables and the other consists of stars just below the TRGB discontinuity. The position of the dip between the two peaks coincides with the above mentioned discontinuity for the main sample.

Alves et al. (1998) and Wood et al. (1999) discussed the MACHO red variables that lie below the tip of the first giant branch. They suggested that all red giant variables are most likely AGB stars. According to Origlia et al. (2002), about 15% of RGB stars near the TRGB have dust shells. Assuming all AGB stars vary and all of these 15% of RGB stars near the TRGB also vary, and using Renzini’s AGB/RGB ratio of 1/3, the ratio of variable RGB stars to AGB stars is at least 0.45. Although the second peak might be made of AGB variables with much fainter luminosities than the typical AGB variables, the coincidence of the dip luminosity for the variables with the discontinuity luminosity for the main sample can be more naturally explained by assuming that the second peak is due to a contribution from RGB variables since there is no reason to assume, according to models, that early AGB variables accumulate at the TRGB (Cioni et al. 2000b).

## 5 Summary

In this letter, we presented the first results of our ongoing NIR survey of variable red giants in the LMC. Many obscured stars ($`JK>2.4`$) were detected and we found that most of them are variables, indicating that the mass-loss and variability are connected. The nature of the obscured variables are still unknown and further observational investigations are necessary and in progress. We also found that luminosity distribution of variable stars has two peaks. One is well above the TRGB and the other is just around it. Because there is no obvious reason to assume that Early AGB variables ”pile-up” around the TRGB, we attributed the second peak to a contribution from RGB variables. From this point of view, we propose the evolutionary scenario of red giants as: (1) pulsate at the tip of the first giant branch; (2) pulsation fades away with the helium flash; (3) pulsation recommences as they ascend the asymptotic giant branch. However, the discrimination between RGB and AGB is still a challenging subject, and more direct evidence for their classification is needed.

## Acknowledgments

The authors would like to thank Michael Feast and John Menzies for a critical reading of the manuscript. We also thank Motonori Kamiya for kindly offering us a picture used in Figure 1. This research is supported in part by the Grant-in-Aid for Scientific Research (C) No. 12640234 and Grant-in-Aids for Scientific on Priority Area (A) No. 12021202 and 13011202 from the Ministry of Education, Science, Sports and Culture of Japan. The IRSF/SIRIUS project was initiated and supported by Nagoya University, National Astronomical Observatory of Japan and University of Tokyo in collaboration with South African Astronomical Observatory under a financial support of Grant-in-Aid for Scientific Research on Priority Area (A) No. 10147207 of the Ministry of Education, Culture, Sports, Science, and Technology of Japan. Y. Nakajima, C. Nagashima and T. Nagayama are financially supported by Japan Society for the Promotion of Science. Finally, we thank the referee for kind and useful comments that helped to improve the paper.


text:
# High 𝑝_𝑇 Measurements from PHENIX

## 1 Introduction

In Run II at RHIC ($`\sqrt{s_{NN}}=200`$ GeV), PHENIX has been able to measure particle yields up to large transverse momenta, $`p_T10`$ GeV/c, where the cross section is expected to be dominated by hard processes. This high statistics data set allows us to study the systematics of high $`p_T`$ phenomena in detail and at much higher $`p_T`$ than possible in Run I ($`\sqrt{s_{NN}}=130`$ GeV) data.

In p+p collisions, it is known that hard scattering and the fragmentation of the scattered partons dominates the production of hadrons above $`p_T2`$ GeV/c . In Au+Au collisions, these hard scatterings are of particular interest because they occur early in the collision, leaving the hard-scattered partons sensitive to the properties of the collision medium. Further interactions of these partons in the dense medium might cause the high transverse momentum tail of the hadron spectrum, where the hadrons are likely to be the leading particles of jets, to be suppressed . Due to the high multiplicity environment, jets cannot be directly observed in a heavy ion collision. We can however determine the contribution of jet fragmentation to hadron yields at high $`p_T`$ via correlation measurements and use this knowledge to interpret the measured hadron $`p_T`$ spectra.

First, we give a brief overview of the PHENIX detector. In the following sections, we proceed assuming that the hadrons at high $`p_T`$ emanate from jets, with the objective of measuring the effect of the dense medium on the hard-scattered partons through modification of the particle spectra at high $`p_T`$. We define the nuclear modification factor and show such a quantity for both neutral pions and charged hadrons. With the measurement of spectra for both charged hadrons and neutral pions, we also investigate the particle composition at high $`p_T`$. Finally, we use two-particle correlations to test the assumption that the high $`p_T`$ hadrons are due to the fragmentation of hard-scattered partons.

## 2 Overview of the PHENIX Detector

PHENIX is a versatile detector, with tracking chambers, ring-imaging Cerenkov (RICH) detectors, a time-of-flight (TOF) wall, and electromagnetic calorimeters in the central arms, which cover $`\mathrm{\Delta }\eta =\pm 0.35`$ and $`180^o`$ in azimuth. Along the beamline there are beam-beam counters and zero-degree calorimeters, which are used for event triggering and centrality classification. Neutral pions are measured via their $`\gamma \gamma `$decay using the electromagnetic calorimeters (EMCal). There are 6 sectors of Lead-Scintillator (PbSc) and 2 sectors of Lead-Glass (PbGl) calorimeter, providing two independent measurements of the $`\pi ^0`$ $`p_T`$ spectrum. The charged hadron spectrum is measured using the tracking detectors: drift chambers and pad chambers. At $`p_T`$ above $`5`$ GeV/c the RICH detector, together with an energy cut in the PbSc, is used to identify charged pions. With various detection methods, PHENIX can make multiple overlapping measurements at high $`p_T`$, providing valuable consistency checks. Details of the detector are described elsewhere .

## 3 High $`p_T`$ Particle Spectra

In Run II, PHENIX recorded more than 30 million minimum bias Au+Au events. This event sample enables us to reach 10 GeV/c in the transverse momentum spectra of both identified neutral pions and charged hadrons. In addition, we recorded 140 million p+p events in which we measure neutral pions up to 13 GeV/c. Shown in Fig. 1 is the neutral pion production cross section as a function of $`p_T`$, measured with the PbSc, for p+p collisions .

The measured p+p spectrum is compared to a fit to UA1 data  over the range $`p_T<6`$ GeV/c, extrapolated to higher $`p_T`$. Although our measurement agrees with the UA1 data, it disagrees with this extrapolation at very high $`p_T`$. This new measurement provides an important reference spectrum for comparison with the Au+Au data.

For Au+Au collisions, the events are binned into centrality selections, and the spectra are shown for the most central bin and a peripheral bin in Figs. 3 and 3. Centrality is expressed as a percentage of the total inelastic cross section. In the Run II data, the charged hadron yields are measured up to $`p_T=10`$ GeV/c for central events and $`p_T=5`$ GeV/c for peripheral events, while the neutral pion yields are measured up to 8 and 6.5 GeV/c for central and peripheral events respectively.

## 4 $`R_{AA}`$ for Neutral Pions

The nuclear modification factor quantifies the effect of A+A compared to p+p collisions on particle yields for point-like processes. It is defined as the ratio of the particle yield in a A+A collision to the yield in a p+p collision scaled by the mean number of binary (nucleon+nucleon) collisions $`N_{coll}`$ in the A+A event sample. Because hard processes are generally believed to scale with $`N_{coll}`$ (“binary-scaling”), this ratio is expected to be one at high $`p_T`$ in the absence of any nuclear effects.

$$R_{AA}(p_T)=\frac{(\mathrm{Yield}\mathrm{per}\mathrm{A}+\mathrm{A}\mathrm{collision})}{N_{coll}(\mathrm{Yield}\mathrm{per}\mathrm{p}+\mathrm{p}\mathrm{collision})}=\frac{d^2N^{A+A}/dp_Td\eta }{N_{coll}(d^2\sigma ^{p+p}/dp_Td\eta )/\sigma _{inelastic}^{p+p}}.$$
(1)
Shown in Fig. 5 is $`R_{AA}`$ as a function of $`p_T`$ for central Pb+Pb collisions at $`\sqrt{s_{NN}}=17`$ GeV from WA98  and central Au+Au collisions at $`\sqrt{s_{NN}}=130`$ GeV from PHENIX Run I . Already in this comparison, the suppression relative to binary-scaling at 130 GeV is striking. It is quite different from the enhancement observed at 17 GeV. This enhancement is known as the “Cronin effect”  and has been attributed to $`p_T`$ broadening due to initial state scatterings . Figure 5 shows $`R_{AA}`$ for central and peripheral Au+Au collisions at $`\sqrt{s_{NN}}=200`$ GeV.

The central data again show a suppression as seen at 130 GeV. With a measurement that extends to much larger $`p_T`$, the suppression is shown to persist up to $`p_T8`$ GeV/c and is as much as a factor of 5-6 at the largest $`p_T`$. The measurement of the p+p reference spectrum with the same detector greatly reduces the systematic error in $`R_{AA}`$ for central events. The $`R_{AA}`$ for peripheral events is consistent with binary-scaling (p+p yields scaled by the number of binary collisions in the peripheral event sample). The error in this ratio is dominated by the uncertainty in $`N_{coll}`$ for peripheral collisions. Within the rather large errors, the spectrum for peripheral collisions does not show effects of the nuclear medium. If there is a Cronin effect at RHIC, it is not strongly evident in peripheral collisions.

There are a number of model predictions, shown in Fig. 6, that were made prior to the release of the Run II data. In the case of one model which assumes energy loss that is constant with respect to the parton’s energy , $`R_{AA}`$ increases by a factor of $`2`$ from $`p_T3.5`$ to $`p_T8`$ GeV/c. Although the magnitude of the suppression qualitatively describes the data at $`p_T4`$ GeV/c, the increasing $`R_{AA}`$ with $`p_T`$ is contrary to what is observed in central Au+Au collisions. Two other models which use a common formalism for including energy-dependent energy loss differ in the point where the calculation is begun (thus resulting in somewhat different predictions). One such model predicts an $`R_{AA}`$ that increases only slightly with increasing $`p_T`$ , while another model predicts a nearly constant $`R_{AA}`$ from $`p_T4`$ to at least $`p_T10`$ GeV/c .

The energy-dependent energy loss models give reasonable agreement with the data. This is also seen in another calculation made subsequent to the release of the Run II data .

## 5 Comparison of Central Collisions with Respect to Peripheral Collisions

An alternate measure of the nuclear modification is the central-to-peripheral ratio scaled by the number of binary collisions,

$$\mathrm{Binary}\mathrm{Scaled}\mathrm{Central}/\mathrm{Peripheral}=\frac{(\mathrm{Yield}\mathrm{per}\mathrm{central}\mathrm{collision})/\mathrm{N}_{\mathrm{coll}}^{\mathrm{central}}}{(\mathrm{Yield}\mathrm{per}\mathrm{peripheral}\mathrm{collision})/\mathrm{N}_{\mathrm{coll}}^{\mathrm{peripheral}}}.$$
(2)
Since the $`\pi ^0`$ spectrum in peripheral Au+Au collisions is consistent with the binary-scaled p+p $`\pi ^0`$ spectrum (Fig. 5), this ratio is similar to $`R_{AA}`$. An advantage of this ratio is that many of the systematic uncertainties in the measurement cancel, particularly for charged hadrons. A disadvantage is that it is sensitive to the centrality dependence of the Cronin effect, which is not known. Figure 7 shows this ratio as a function of $`p_T`$ for neutral pions, charged hadrons, and charged pions.

The suppression is observed in all three measurements. In the charged hadrons, the suppression factor reaches 3-4 at high $`p_T`$. The difference in this ratio between the identified pions and charged hadrons seems to be due to the particle composition. In particular, we find that $`p/\pi ^01`$ in central collisions for $`p_T=24`$ GeV/c while $`p/\pi ^00.4`$ in peripheral collisions . The error in the normalization, denoted by a shaded area for the charged hadrons and outlines for the neutral pions, is dominated by the uncertainty in $`N_{coll}`$ for peripheral collisions.

## 6 Centrality Dependence of $`R_{AA}`$

As shown in Fig. 5, the 10% most central events have a suppression factor reaching 5-6, while the 70-80% centrality selection is consistent with binary-scaled p+p collisions. We now address the evolution of the nuclear modification factor from the most peripheral to the most central collisions, in terms of the mean number of participating nucleons $`N_{part}`$ in a centrality selection. For the 10% most central collisions, $`N_{part}=327`$, and for the 70-80% central collisions, $`N_{part}=14`$. Figure 9 shows $`R_{AA}`$ vs. $`N_{part}`$ for neutral pions with $`4<p_T<6`$ GeV/c. The suppression gradually increases from peripheral to central events. The same trend is seen in the charged hadron yields. Figure 9 shows the integrated yield scaled by $`N_{coll}`$ vs. $`N_{part}`$ for charged hadrons. This quantity is similar to $`R_{AA}`$, differing only in normalization by the p+p reference yields.

Again the binary-scaled yield decreases, or the suppression increases, gradually with increasing $`N_{part}`$.

## 7 Particle Composition at High $`p_T`$

PHENIX can measure the $`p/\pi `$ ratio up to almost 4 GeV/c in $`p_T`$ , at which point protons can no longer be distinguished from other particles via their time of flight. Since neutral pions are identified to much higher $`p_T`$, we can look at the ratio of $`\pi ^0`$ to nonidentified charged hadrons $`(h^++h^{})/2`$ for $`p_T>4`$ GeV/c. This ratio is shown in Fig. 10 for minimum bias events. The surprising feature that this ratio does not increase for transverse momenta greater than 4 GeV/c, but remains nearly constant at a value around 0.5, indicates that approximately half of the charged hadrons at high $`p_T`$ are protons and/or kaons, assuming $`\pi ^0=(\pi ^++\pi ^{})/2`$. This result is rather different from hard processes as measured in e+e- collisions . The interpretation for this new result is not yet clear. There may be some other production mechanism for protons and/or kaons in Au+Au collisions at these large transverse momenta. Alternatively, the large proton and/or kaon content at high $`p_T`$ could be due to a difference in the suppression of pions relative to protons/kaons.

## 8 Evidence for Jets

Since jets cannot be directly observed in a Au+Au collisions, a correlation analysis is required to detect their presence. We take a high $`p_T`$ neutral particle, a cluster in the EMCal with energy greater than 2.5 GeV, as the trigger and correlate it with all charged particles in the event within a given range in $`p_T`$. Comparing the correlations that we measure in p+p collisions at $`\sqrt{s_{NN}}=200`$ GeV to those from the PYTHIA event generator , we establish that PYTHIA reproduces the behavior of jets known to be present in p+p collisions. PYTHIA is then compared to Au+Au collisions at $`\sqrt{s_{NN}}=200`$ GeV, which demonstrates that the correlations we observe in Au+Au collisions behave as expected for jets, both in the azimuthal and polar angles. Figure 11 shows the background-subtracted correlation in the azimuthal angle $`\mathrm{\Delta }\varphi `$ of the trigger particle with charged particles having $`p_T`$ between 2 and 4 GeV/c.

Also shown is a fit to the data of the correlation produced by PYTHIA with an added $`cos(2\varphi )`$ term. The $`cos(2\varphi )`$ term accounts for the contribution of flow to the correlation, while the coefficient of the PYTHIA term represents the contribution due to jets. Further details of this analysis are given elsewhere . From the figure, one can see that much of the correlation can be attributed to jets (PYTHIA), particularly the near-side correlation around $`\mathrm{\Delta }\varphi =0`$. The trigger particle is typically a photon from a $`\pi ^0`$ decay. This is strong evidence that when we trigger on a high $`p_T`$ $`\pi ^0`$, we are indeed looking at a leading particle from a jet.

## 9 Conclusions

We have presented $`p_T`$ spectra for charged hadrons in Au+Au collisions and neutral pions in Au+Au and p+p collisions and have shown that the high $`p_T`$ suppression observed in Run I RHIC data at $`\sqrt{s_{NN}}=130`$ GeV  persists up to $`p_T8`$ GeV/c at $`\sqrt{s_{NN}}=200`$ GeV. The neutral pions are suppressed by a factor of 5-6 and the charged hadrons by 3-4 for $`p_T>4`$ GeV/c in the most central collisions. The suppression increases gradually with increasing centrality, or $`N_{part}`$. We have also shown that when triggering on a high $`p_T`$ neutral particle, which is predominantly a photon from a $`\pi ^0`$ decay, the correlation between the neutral and charged particles in the event is jet-like. This is strong evidence that the neutral pions that we measure at high $`p_T`$ indeed have significant contributions from jet fragmentation. Finally, we have shown that approximately half of the charged hadrons at high $`p_T`$ are not pions, but protons and/or kaons, up to $`p_T9`$ GeV/c.