format: text

text:
# HST/NICMOS Observations of IZw18 – revealing a population of old AGB stars1footnote 11footnote 1Based on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.

## 1. Introduction

The blue compact galaxy (BCG) IZw18 (also known as Markarian 116) was first described by Zwicky in 1966. In a seminal paper Searle and Sargent (ss (1972)) showed that IZw18 had a very low oxygen abundance, and they concluded that the galaxy either must be young or form stars in short bursts intervened by long quiescent epochs (cf. Searle et al. ssb (1973)). This has been followed by many studies of the chemical properties of IZw18 and BCGs in general. Alloin et al. (alloin (1978)) derived an oxygen abundance of 12+log(O/H)=7.2 (2% solar) for IZw18 in agreement with more recent values, e.g. Izotov and Thuan (it (1999)). In spite of many searches for more metal-poor galaxies, IZw18 remains the most metal-poor one (when nebular oxygen abundances are considered) and is often regarded as the prototypical BCG (see Kunth and Östlin ko (2000) for a review).

If IZw18 was a genuinely young galaxy, that would imply that galaxies can form still at the present cosmic epoch and considerable effort has therefore been devoted to examine the ages of BCGs. Although it is now clear that most BCGs are old (see Kunth and Östlin ko (2000)), IZw18 has remained among a small number of galaxies which show no clear signs for the presence of an old stellar population. Kunth and Sargent (ks (1986)) argued that IZw18 could be young and that the low, but non-zero, metallicity could be due to self-pollution by massive stars. Recently an extensive study of chemical abundances of BCGs was published by Izotov and Thuan (it (1999)). They found that for very low oxygen abundances (less than 5% solar) the N/O and C/O ratios stay constant and have very little scatter. Their interpretation is that N and C are produced entirely by massive ($`M10M_{}`$) stars and that this implies that these galaxies must be younger than 40 Myr. While other interpretations are possible (Kunth and Östlin ko (2000)), this an interesting idea that deserves to be tested. With its record low oxygen abundance, IZw18 provides the critical test case.

The chemical abundances of IZw18 have been used to infer both a high and a low age (Kunth and Östlin ko (2000)), and it is clear that chemistry alone has limited value in constraining the age and star formation history given the many uncertain physical parameters involved (winds, yields, infall, initial mass function, etc.). Stronger constraints would be obtained by direct observations of the resolved stellar population. Hunter and Thronson (hunter:thronson (1995)) used deep HST/WFPC2 imaging to investigate the stellar population of IZw18, and found a population of young massive stars. No evidence for old stars was found, but neither any against. Similar results were obtained by Dufour et al. (dufour (1996)) from an independent WFPC2 data set. Lately, these two WFPC2 data sets have been carefully reanalyzed by Aloisi, Tosi and Greggio (aloisi (1999)), who discovered a population of red stars. Their results indicate that IZw18 is at least 0.5 Gyr old. This was however disputed by Izotov et al. (izotov (1999)) who claim that the distance to IZw18 has been underestimated, and that hence the age has been overestimated.

In this paper I report on the first results from a near infrared imaging study of IZw18 using NICMOS on board the HST. The goal was to search for an old stellar population, by identifying luminous red giants and asymptotic giant branch (AGB) stars. The data is however also valuable for studying the young massive stars, e.g. red super giants (RSGs) in the near infrared (NIR). A follow up program designed at investigating the presence of carbon stars by targeting molecular features were performed during Cycle 7-NICMOS. The results of this investigation will be presented in a forthcoming paper.

## 2. Observations and Reductions

I Zw 18 was observed for 3 orbits during Cycle 7 with NICMOS in the F110W and F160W passbands (roughly corresponding to J and H) under GO program 7461. The instrument and its performance are described in Thompson et al. (thompson (1997)) and the NICMOS instrument handbook (MacKenty et al. mackenty (1997)). The NIC2 camera was chosen to give good pixel sampling (0.075 $`\mathrm{}`$/pixel) and an acceptable field of view ($`19\mathrm{}\times 19\mathrm{}`$). All science exposures were obtained with the MULTIACCUM mode which employs consecutive non-destructive readouts of the detector. This enables e.g. an efficient cosmic-ray removal. For the F110W filter, two MIF1024 sequences were acquired, separated by an offset of 10 integer pixels in the Y direction. In addition, one MIF512 sequence was obtained with a 0.5 pixel offset in both the X and Y directions with respect to the first MIF1024 sequence. The use of a non-integer pixel offset was adopted since the NIC2 pixel scale slightly undersamples the point spread function (PSF) in F110W. For the F160W filter, two sets of images, consisting of one MIF512 and two MIF1024 sequences, were acquired at the same positions as the two MIF1024 sequences in F110W. The total exposure times are thus 2560 and 5120 seconds in the F110W and F160W filters respectively.

Inspection of the pipeline processed images revealed that these suffered from the so called “pedestal” effect, a residual bias not removed by the pipeline. To subtract the pedestal, the Pedestal Estimation and Quadrant Equalization Software developed by R.P. van der Marel was successfully utilized. After the pedestal had been removed, the images in each pass band were combined using the Drizzle package (Fruchter and Hook fruchter (1998)). The resulting F160W image is shown in Fig. 1.

Photometry of point sources was performed with the daophot package (Stetson et al. stetson (1990)). Unfortunately, the images contain few suitable stars for modeling the point spread function (PSF), since the bright stars are found in crowded regions. One bright source 4″ south-west (SW) of the main body turned out to be spatially extended and is probably a background galaxy. To assure that the wings of the PSF were not contaminated by neighboring stars, photometry was also performed using a synthetic PSF, which was constructed using the TinyTim software (version 4.4, Krist and Hook krist (1997)). The resulting photometry using real and synthetic PSF stars were nearly identical, except for a systematic difference of the order of 0.05 magnitudes in the F160W filter; the use of synthetic PSF stars yielded redder stars. Comparison with a control sample consisting of the eight best isolated stars showed that the synthetic PSF photometry did not suffer from any systematic offset compared to the aperture photometry, whereas the use of a natural PSF showed a small systematic offset but smaller scatter than the synthetic PSF photometry. Thus, the natural PSF photometry is probably more accurate but subject to a small systematic offset. No systematic difference was seen in the F110W data. The natural PSF photometry was finally adopted after correcting for the small offset in F160W. After the photometric reduction, stars that resulted too sharp or too extended in at least one filter were rejected. The photometries in the two filters were then combined by considering a matching radius of 0.5 pixels in X and Y for the stellar centers in the F110W and F160W images. All objects were carefully analyzed by eye and dubious sources rejected.

The DAOPHOT magnitudes refer to an aperture with radius 0.1875″. The aperture corrections in going from this to an 0.5″ aperture were determined using synthetic PSF and amount to 0.24 and 0.48 magnitudes for F110W and F160W respectively. To obtain total magnitudes an aperture correction of $`0.156`$ magnitudes for both filters were added to the the 0.5″ aperture magnitudes, and finally, the zeropoints of $`ZP_{F110W}=22.5436`$ and $`ZP_{F160W}=21.826`$ were added in order to have the final calibrated magnitudes in the HST Vegamag system (see NICMOS WWW pages). I will use the designations $`m_{110}`$ and $`m_{160}`$ to express the magnitudes in F110W and F160W respectively, and $`\mathrm{F110W}\mathrm{F160W}`$ for the color. The resulting color–magnitude diagrams (CMDs) and the magnitude–error relations are shown in Fig. 2. The detected stars span a similar range in color and brightness as those detected with HST/NICMOS in VIIZw403 (Schulte-Ladbeck et al. schulte (1999)).

The photometric completeness was investigated by using the addstar task. Artificial stars were added in half magnitude bins at random locations in two separate regions: “SE” covering the mildly crowded south-east region; and “NW” which is the severely crowded north-west region. Then, repeated tests were made of how many test stars within each half magnitude bin that were successfully recovered by daophot. More than 100 stars were added for each bin. The 75 and 50 % completeness levels are summarized in Table 1. In the NW region there are still bright unresolved sources left after the repeated daophot runs.

## 3. The distance of IZw18

Usually a distance of 10 Mpc is assumed for IZw18, based on a radial velocity of $``$750 km/s and a Hubble constant of $`H_0=75`$km/s/Mpc. Recently Izotov et al. (izotov (1999)) have claimed that the distance to IZw18 has been severely underestimated, and if Virgo-centric flow is taken into account a distance of 15-20 Mpc is obtained. They present also other circumstantial evidence (the presence of Wolf-Rayet stars and the ionization state of IZw18) for a distance $`D20`$ Mpc. It is true that Virgo-centric motion will slightly increase the distance, but the distance proposed by Izotov et al. (izotov (1999)) refer to old values of the distance of the Virgo cluster (21.7 Mpc). A lot of recent progress has been made on the distance to Virgo (e.g. Macri et al. macri (1999)) showing that its distance is $`D_{Virgo}16`$ Mpc. The Virgo-centric flow models express the distance of a galaxy as $`x`$, relative to the Virgo cluster, i.e. $`D=x\times D_{Virgo}`$ . Applying the linear model by Schechter (schechter (1980)) gives $`x=0.7`$ while the non-linear models by Kraan-Korteweg (kraan (1986)) give $`x=0.8`$, corresponding to 11 and 13 Mpc respectively. Hereafter I will adopt $`D=12.6`$ Mpc for IZw18, which corresponds to a distance modulus of $`(mM)=30.5`$. The distance should not be much larger (maximum 14.5 Mpc given the Virgocentric flow models and the uncertainty in $`D_{Virgo}`$) than this, but may be as low as $`10`$ Mpc. The distances suggested by Izotov et al. (izotov (1999)) would require IZw18 to be more distant than the Virgo cluster, in spite of a much lower radial velocity.

## 4. Analysis and Discussion

Internal extinction has only a moderate influence on the F110W-F160W colors and in general $`E(\mathrm{F110W}\mathrm{F160W})=0.5\times E(BV)`$, which has been estimated with the Synphot software. All optical investigations have found a low internal extinction in IZw18. I have adopted $`E(BV)=0.05`$ throughout the galaxy, and this correction has been applied to the data-points in Fig. 3.

The obtained color magnitude diagram (CMD) clearly presents two distinct features. One population of red luminous stars and one with fainter ($`m_{110}>23`$) red stars. This suggests that stars with at least two different ages are present. To interpret the observed CMD I have mainly used the Padua stellar models for a metallicity $`Z=0.02Z_{}`$ (Fagotto et al. fagotto (1994), Bertelli et al. bertelli (1994)), i.e. identical to the gas metallicity of IZw18. In general, the Bertelli et al. (bertelli (1994)) isochrones are more useful than the Fagotto et al. (fagotto (1994)) tracks since the former contains later evolutionary stages, as the thermally pulsating AGB phase, which is very important for the present study.

Before interpreting the CMD, the effective temperature ($`T_{\mathrm{eff}}`$) and luminosity ($`L`$) of the stellar models have to be related to colors and absolute magnitudes in the observed system. To translate $`T_{\mathrm{eff}}`$ into F110W-F160W colors, and $`L`$ into $`M_{110}`$ for a given surface gravity, the BCP models by Origlia and Leitherer (ol (2000)) were used for a metallicity of $`[\mathrm{Fe}/\mathrm{H}]=2`$.

In figure 3 the Bertelli et al. (bertelli (1994)) isochrones are over-plotted on the CMD. For a distance modulus of $`(mM)=30.5`$ it is obvious that stars of ages 11-20 Myr are required to reproduce the luminous upper branch in the CMD ($`m_{110}<23.2,F110WF160W<1`$). These ages equals the lifetime of stars with initial masses of 20 to 13 $`M_{}`$. The upper branch is dominated with stars $`14`$ Myr old when the current star formation event appears to have peaked. For comparison, distances of 10 and 15 Mpc would give ages of 18 and 12 Myr respectively. Comparison with the Fagotto et al. (fagotto (1994)) stellar evolutionary tracks yield the same conclusion that the most luminous red stars are the evolved descendants of stars with initial masses 15-20 $`M_{}`$, and the same results are obtained when using “Geneva” tracks for the same metallicity (Schaerer 1999, private communication; cf. de Mello et al. demello (1998)).

For low mass stars (initial mass $`M5M_{}`$, age $`>100`$ Myr) only the isochrones (Bertelli et al. bertelli (1994)) have been considered, since the tracks do not include the thermally pulsating AGB phase. The most interesting result is that to reproduce the faint red part of the CMD ($`m_{110}>23.5,F110WF160W>1`$), an age of 1 Gyr is required at the distance of 12.6 Mpc. Even if the distance would be as high as 15 Mpc, stars older than 0.5 Gyr would be required, and $`D=10`$ Mpc would lead to a much higher age than 1 Gyr. Not even at 20 Mpc can the faint red stars be made younger than 300 Myr. For $`D=12.6`$ Mpc the red faint population can be nicely explained as the upper part of AGB stars with a range of ages from 0.3 to 5 Gyr. This corresponds to the lifetime of stars with initial masses 2.5 to 1.1 $`M_{}`$. A single 0.5 Gyr old AGB population does not suffice since that would not produce enough faint ($`m_{110}24`$) red stars, only when extending to ages in excess of 1 Gyr can the CMD of faint stars be satisfactory explained. This indicates that IZw18 has experienced more or less continuous star formation for at least 1 Gyr prior to the present burst.

The distribution of stars in the CMD is broadened by the observational error function, and it has to be assured that it cannot be reproduced by young stars plus photometric errors. If one assumes that no star is older than 100 Myr, what distribution would one expect? An isochrone of that age has a color $`F110WF160W=0.79`$ at $`m_{110}=24`$. The one sigma photometric errors at $`m_{110}=24,F110WF160W1`$ are $`\sigma _{F110WF160W}=0.125`$, but let’s be very conservative and adopt $`\sigma _{F110WF160W}=0.2`$. The one sigma range around $`F110WF160W=0.79`$, $`m_{110}24`$ then contain 14 stars. For a normal distribution, 68% of the stars should be contained within this range, and one would thus expect to find only $`3`$ stars (i.e. 16%) with $`F110WF160W1.0`$, wheras one finds more than a dozen. Thus the observed distribution of red stars cannot be reproduced by 0.1 Gyr old stars plus photometric errors. Even if rejecting the reddest stars ($`>1.5`$), an age of 1 Gyr is required to be consistent with the error distribution. Since this was a conservative estimate, we can conclude that IZw18 has an age of at least 1 Gyr. This result on the near-infrared CMD analysis is confirmed when using the synthetic CMD technique (cf. Tosi et al. (tosi (1991)); Aloisi et al. aloisi (1999)): stars older than 1 Gyr are required in order to reproduce the observed CMD, the best fit being obtained by including stars as old as 5 Gyr (Östlin 2000, in preparation).

An analysis of the spatial distribution of stars of different age show that young red super giants occupy preferentially the NW region, while the old AGB stars are more evenly distributed. Further constraints on the age of the stellar population may be obtained from the color of the underlying galaxy after the found stars have been subtracted. For a single stellar population, the integrated F110W-F160W color increases with age, from $`F110WF160W=0.1`$ at 10 Myr up to $`F110WF160W=0.6`$ at 125 Myr when the AGB stars start to dominate, and stays roughly constant thereafter. This has been calculated using the PEGASE spectral synthesis code (Fioc & Rocca-Volmerange 1999) and the NICMOS F110W and F160W response curves. The integrated colors in $`15\times 15`$ pixels boxes are typically $`F110WF160W=0.3`$ to $`0.6`$ in the NW and $`F110WF160W>0.6`$ in the SE region. This indicates that in the NW region, the light is dominated by an unresolved population of stars with ages 20 to 100 Myr, while in the SW region, the underlying population is dominated by stars older than 125 Myr.

The result of an age of 1 Gyr or larger for the underlying population in IZw18 is confirmed also by deep optical/ near-IR surface photometry (Östlin 2000, in preparation; Kunth & Östlin ko (2000)) and is in good agreement with the results of Aloisi et al. (aloisi (1999)) for the optical CMDs. We can thus conclude that all these different studies on independent sets of data in various spectral regions give the same answer: that IZw18 is an old galaxy.

## 5. Summary and Conclusions

I have analyzed the near infrared color-magnitude diagram from HST/NICMOS images of IZw18, the most metal-poor galaxy known. Assuming a distance of 12.6 Mpc, I find support for the presence of AGB stars at least 1 Gyr old. The results presented in this letter hold qualitatively even if the distance to IZw18 is significantly different. The latest star formation event appears to have peaked 14 Myr ago. My results are in good agreement with Aloisi et al. (aloisi (1999)) who from a careful reanalysis of HST/WFPC2 archival data concluded that IZw18 must be at least 0.5 Gyr old at their adopted distance of 10 Mpc. Thus, there is no need for galaxies with oxygen abundance less than 5 % of the solar value to be young as was claimed by Izotov and Thuan (it (1999)). This closes the question whether IZw18 is a genuinely young galaxy, still in the process of formation: the answer is – No!

Thanks to D. Calzetti and W. Freudling for practical help and suggestions, and to D. Schaerer for making revised evolutionary tracks available ahead of publication. N. Bergvall, D. Kunth, A. Lançon, M. Longhetti, M. Mouhcine and R. Schulte-Ladbeck are thanked for useful discussions. This work made use of the Pedestal Estimation and Quadrant Equalization Software developed by Roeland P. van der Marel. The referee A. Aloisi is thanked for many useful suggestions. This work was financially supported by the Swedish foundation for international cooperation in research and higher education (STINT), and the Swedish National Space-board. Many thanks to Anna and Snigel.


text:
# INTERSTELLAR HYDROGEN AND COSMIC BACKGROUND RADIATION

## Abstract

It is shown that a collection of photons with nearly the same frequency exhibits a ”condensation” type of phenomenon corresponding to a peak intensity. The observed cosmic background radiation can be explained from this standpoint in terms of the radiation due to fluctuations in interstellar Hydrogen.

00footnotetext: E-mail:birlasc@hd1.vsnl.net.in

In a previous communication it was suggested that the origin of the Cosmic Background Radiation is the random motion of interstellar Hydrogen. We will now deduce the same result from a completely different point of view.
We start with the formula for the average occupation number for photons of momentum $`\stackrel{}{k}`$ for all polarizations:

$$n_\stackrel{}{k}=\frac{2}{e^{\beta \mathrm{}\omega }1}$$
(1)
Let us specialize to a scenario in which all the photons have nearly the same energy so that we can write,

$$n_\stackrel{}{k}=n_\stackrel{}{k}^{}\delta (kk^{}),$$
(2)
where $`n_\stackrel{}{k}^{}`$ is given by (1), and $`k|\stackrel{}{k}|.`$ The total number of photons $`N`$, in the volume $`V`$ being considered, can be obtained in the usual way,

$$N=\frac{V[k]}{(2\pi )^3}_o^{\mathrm{}}𝑑k4\pi k^2n_k$$
(3)
where $`V`$ is large. Inserting (2) in (3) we get,

$$N=\frac{2V}{(2\pi )^3}4\pi k^2[ϵ^\theta 1]^1[k],\theta \beta \mathrm{}\omega ,$$
(4)
In the above, $`[k][L^1]`$ is a dimensionality constant, introduced to compensate the loss of a factor $`k`$ in the integral (3), owing to the $`\delta `$-function in (2): That is, a volume integral in $`\stackrel{}{k}`$ space is reduced to a surface integral on the sphere $`|\stackrel{}{k}|=k^{},`$ due to our constraint that all photons have nearly the same energy.
We observe that, $`\theta =\mathrm{}\omega /KT1,`$ since by (2), the photons have nearly the same energy $`\mathrm{}\omega `$. We also introduce,

$$\upsilon =\frac{V}{N},\lambda =\frac{2\pi c}{\omega }=\frac{2\pi }{k}\text{and}z=\frac{\lambda ^3}{\upsilon }$$
(5)
$`\lambda `$ being the wave length of the radiation. We now have from (4), using (5),

$$(e1)=\frac{\upsilon k^2}{\pi ^2}[k]=\frac{8\pi }{k^{}z}[k]$$

Using (5) we get:

$$z=\frac{8\pi }{k^{}(e1)}=\frac{4\lambda }{(e1)}[k]$$
(6)
From (6) we conclude that, when

$$\lambda =\frac{e1}{4}=0.4[L]$$
(7)
then,

$$z1$$
(8)
or conversely.
Though equation (7) is quite general, as it stands, we have to assign suitable units to it depending on the particular physical situation: We must get an additional input, in the order of magnitude sense, from the system under consideration to fix the units.
Let us now consider the case of radiations due to fluctuations of the cold interstellar Hydrogen as considered from an alternative view point in ref.. In this case it is infact known that $`\frac{1}{\upsilon }1`$ molecule per c.c.. On the other hand, the energy range for these cold molecules is small, so that the above considerations apply. So from (7), owing to the fact that $`\upsilon ^{1/3}\lambda cm.,`$ it follows that,

$$\lambda =0.4cm.$$
(9)
Remembering that from (5), $`\lambda `$ is the wave length and $`\upsilon `$ is the average volume per photon, the condition (8) implies that all the photons are very densely packed as in the case of Bose condensation. This means that from (9), we conclude that at the wave length $`0.4cm,`$ in the micro-wave region, the radiation has a peak intensity. It is remarkable that the cosmic background radiation has the maximum intensity exactly at the wave length given by (9). So, even without the Big Bang event it is possible to receive the observed cosmic background radiation due to fluctuations in interstellar Hydrogen (cf.ref. and for Hoyle and Wickramasinghe’s attempt to explain the background radiation in terms of Helium synthesis, without invoking the Big Bang.)
The same conslusion can be obtained from yet another argument also. Interestingly, El Naschie has studied the Cosmic Background Radiation from the point of view of fractal quantum space time.


text:
# Ferromagnetism and phase separation in one-dimensional d–p and periodic Anderson models

## I Introduction

The nature of the microscopic description of metallic ferromagnetism, as found in transition metals such as iron, nickel and cobalt, is a long-standing problem in strongly correlated electron systems. The Hubbard model was formulated in the 1960’s  in order to describe ferromagnetism in such materials. However, a ferromagnetic ground state in the single–band Hubbard model on nonfrustrated bipartite lattices has not been found at physical parameter values, and can be excluded in a large portion of the ground–state phase diagram. While mean–field theory yields a large region of stable ferromagnetism, fluctuations tend to destabilize the ferromagnetic phase, and numerical and variational calculations have narrowed the possible extent of a ferromagnetic state to a small region around the Nagoaka point. Therefore, it has become clear that additional features must be added to provide a description of metallic ferromagnetism.

Recently, a number of such possible extensions to the Hubbard model have been investigated. In particular, there are three classes of additions which can enhance ferromagnetism: (i) a change in the noninteracting density of states through the addition of frustrating hopping terms or the treatment of a geometrically frustrated lattice (ii) the inclusion of multiple orbitals per site or unit cell, and (iii) the addition of more general nearest–neighbor interactions.  In relation to (i), Vollhardt et al.  have emphasized that peaks or singularities in the occupied portion of the noninteracting density of states minimizes the kinetic energy loss due to polarization of the electrons. This picture is supported by numerical calculations in one  and two dimensions, and within the Dynamical Mean-Field Theory with a number of different forms of the noninteracting density of states.  In addition, Mielke  has shown rigorously that a Hubbard model with a less than half–filled flat band has a fully polarized ground state. This theorem has been extended to the case of nearly flat bands.  For (ii), a number of multiband models, such as the periodic Anderson model, the Kondo lattice model, and the multiband Hubbard models with Hund’s rule coupling have been found to have ferromagnetic ground states. In relation to (iii), Strack and Vollhardt  have studied a quite general form of the nearest–neighbor interaction and have found a fully polarized ground state at half–filling and one hole from half–filling (the Nagaoka state) under certain conditions.

In this paper, we consider a one–dimensional model which contains elements of points (i) and (ii) from above. The model which we will study is a Cu–O chain that contains one oxygen $`p`$–orbital and and copper $`d`$–orbital. Since one–dimensional models can usually be treated using well-controlled analytical or numerical techniques, one can obtain a complete picture of the ground–state phase diagram without resorting to approximations necessary in two or three dimensions. For one–dimensional systems with nearest–neighbor hopping and an arbitrary, but real and particle–symmetric interaction, a theorem by Lieb and Mattis  rules out a magnetized ground state. Because the Lieb–Mattis theorem no longer applies when electrons can pass around each other, the minimal addition to the one–dimensional Hubbard model necessary to obtain ferromagnetism is a next–nearest–neighbor hopping. The resulting Hubbard model on a “zigzag” ladder has a ferromagnetic ground state in a substantial parameter regime.  One can also formulate more general, related models composed of triangular elements. Two such models, studied by Tasaki  and Penc et al.  have been found to have ferromagnetic phases. The $`d`$$`p`$ model we study here is (aside from an added near–neighbor Coulomb repulsion) a particular case of the model studied in Ref. .

Another approach to ferromagnetism in one–dimensional models, involving point (ii) from above, are models with multiple non-degenerate orbitals per site. The simplest of these models are the one–dimensional Periodic Anderson Model (PAM) and the Kondo Lattice model (KLM). Both of these models have ferromagnetic ground states in a relatively large parameter regime.  These models can be related to one another in the regime of large $`f`$–repulsion for the PAM and small Kondo coupling for the KLM. The PAM is similar to the $`d`$$`p`$ model in that there are two linked non-degenerate levels per unit cell, and, as we will discuss in more detail in the following, exhibit similar behavior with the appropriate choice of parameters.

Another motivation for studying the $`d`$$`p`$ model is the relationship with models for the CuO2 planes of the high–$`T_c`$ superconductors. One dimensional $`d`$$`p`$ models have been thought to contain some of the essential ingredients needed to describe superconductivity in CuO2 planes.  Because of this, previous work on the one–dimensional $`d`$$`p`$ model has concentrated on determining the conditions under which superconducting correlations are dominant for models with  and without  direct hopping between the $`p`$–orbitals.

In this work, we will concentrate on the ferromagnetic phases in the ground state phase diagram and will discuss their physical origin. In Sec. II, we introduce the models which we will study and discuss choices of the models parameters. In Sec. III A, we compare the phase diagram of the Anderson lattice model from Ref.  with that of the $`d`$$`p`$ model. We find that the ferromagnetic region persists in the $`d`$$`p`$ model but it is pushed to higher values of the interaction and higher fillings. There is a region of phase separation in which ferromagnetic and antiferromagnetic states coexist. In Sec. III B, we set the parameters of the $`d`$$`p`$ system to realistic values for the cuprates. We then investigate the phase diagram as a function of the oxygen–oxygen hopping amplitude, $`t_{pp}`$, and the band filling, and find three disconnected ferromagnetic regions with three physically different mechanisms.

## II  Models

We consider a one-dimension $`d`$$`p`$ system with Hamiltonian

$`H=`$ $`t_{pp}`$ $`{\displaystyle \underset{[jj^{}]\sigma }{}}(p_{j\sigma }^{}p_{j^{}\sigma }^{}+p_{j^{}\sigma }^{}p_{j\sigma }^{})+\mathrm{\Delta }{\displaystyle \underset{j\sigma }{}}n_{j\sigma }^p+U_p{\displaystyle \underset{j}{}}n_j^pn_j^p`$ (1)
$`+`$ $`U_d{\displaystyle \underset{i}{}}n_i^dn_i^dt_{pd}{\displaystyle \underset{<ij>\sigma }{}}(d_{i\sigma }^{}p_{j\sigma }^{}+p_{i\sigma }^{}d_{j\sigma }^{})+V_{pd}{\displaystyle \underset{<ij>}{}}n_i^dn_j^p,`$ (2)
where $`[jj^{}]`$ denotes a sum over nearest–neighbor oxygen pairs, and $`ij`$ a sum over copper–oxygen nearest neighbors. In the $`d`$$`p`$ system we work in the hole representation so that $`d_{i,\sigma }^{}`$ $`(p_{j,\sigma }^{})`$ creates a hole on $`d`$ ($`p`$) site $`i`$ ($`j`$) with spin $`\sigma `$, and $`n_{i\sigma }^d=d_{i\sigma }^{}d_{i\sigma }^{}`$ and $`n_{j\sigma }^p`$ are the local hole densities on a copper and oxygen site, respectively. The parameter $`t_{pp}`$ is the strength of the direct hopping between the $`p`$-orbitals, $`\mathrm{\Delta }`$ is the difference in on–site energies, $`U_p`$ and $`U_d`$ are the on–site Coulomb repulsion on the $`p`$ and $`d`$ sites respectively, $`t_{pd}`$ is the hybridization and $`V_{pd}`$ the Coulomb repulsion between nearest–neighbor $`p`$ and $`d`$ orbitals. The lattice structure of the $`d`$$`p`$ model including a schematic representation of the Hamiltonian parameters is shown in Fig. 1(a). Unless otherwise stated, we will consider lattices consisting of $`N`$ copper sites and $`N+1`$ oxygen sites, with open boundary conditions at the ends of the lattice so that the ends consist of oxygen sites with connections in only one direction. When the system has $`N_h`$ total holes, we will discuss $`N_pN_hN`$, which corresponds to the excess number of holes in the oxygen band when $`U_d>\mathrm{\Delta }>0`$, and the corresponding $`p`$–band filling, $`n_pN_p/N`$.

This model has a gauge symmetry which can be used to permute the signs of the hopping matrix elements around each triangular element. For example, alternating the sign of $`t_{pd}`$ on succesive bonds and taking $`t_{pp}`$ positive is equivalent to having the same sign of $`t_{pd}`$ and taking $`t_{pp}`$ negative. Here we take the latter case, as shown in Eq. (2), restrict $`t_{pd}`$ to be positive, and allow the sign of $`t_{pp}`$ to vary.

In order to understand the physics of ferromagnetism in the $`d`$$`p`$ model, it is useful to make a comparison with a model that has a similar structure in an appropriately chosen parameter regime, the one–dimensional periodic Anderson model. The PAM has Hamiltonian

$$H=t\underset{i\sigma }{}(c_{i\sigma }^{}c_{i+1\sigma }+c_{i+1\sigma }^{}c_{i\sigma })+\epsilon _f\underset{i\sigma }{}n_{i\sigma }^f+U\underset{i}{}n_i^fn_i^f+V\underset{i\sigma }{}(c_{i\sigma }^{}f_{i\sigma }+f_{i\sigma }^{}c_{i\sigma }),$$
(3)
where $`c_{i\sigma }^{}`$ and $`c_{i\sigma }`$ create and annihilate conduction electrons with spin $`\sigma `$ at lattice site $`i`$, and $`f_{i\sigma }^{}`$ and $`f_{i\sigma }`$ create and annihilate local $`f`$–electrons. Here $`t`$ is the hopping matrix element for conduction electrons between neighboring sites, $`\epsilon _f`$ is the energy of the localized $`f`$–orbital, $`U`$ is the on–site Coulomb repulsion of the $`f`$–electrons, and $`V`$ is the on–site hybridization matrix element between electrons in the $`f`$–orbitals and the conduction band. We denote the number of electrons by $`N_e`$, and $`N`$ is the number of unit cells (each consisting of one $`f`$ site and one conduction site) in the lattice. Since there are two electronic orbitals in each site, the quarter–filled case corresponds to $`N_e=N`$ and the half–filled case has $`N_e=2N`$. In analogy to the $`d`$$`p`$ model, we will discuss $`N_c=N_eN`$, the excess number of electrons in the conduction band when $`U>\epsilon _f>0`$, and the corresponding filling, $`n_cN_c/N`$. Fig. 1(b) shows a schematic represention of the PAM.

One can set the parameters of $`d`$$`p`$ model so that it is the same as the PAM except that the $`d`$-orbital hybridizes with the two nearest–neighbor copper atoms, while the PAM hybridizes only on–site. In order to compare the two models, we set $`V_{pd}=U_p=0`$, $`t_{pp}=t`$, $`t_{pd}=V`$, $`\mathrm{\Delta }=ϵ_f`$ and $`U_d=U`$.

## III  Calculation and Results

Since itinerant ferromagnetism is an intrinsically strong–coupling phenomenon that occurs due to a subtle competition between kinetic and potential energy,  it is generally quite hard to treat with analytic methods. While exact statements can be made in some special cases, such as the Nagaoka state  and the case of flat bands, in general one must resort to mean–field theory or variational techniques, both of which give phase diagrams which can be at best qualitatively accurate. Therefore, sufficiently accurate numerical methods can be useful to determine the properties of models with ferromagnetic phases. Here we use the Density Matrix Renormalization Group to calculate the ground–state properties of the $`d`$$`p`$ model and the PAM. The DMRG is a variational numerical method closely related to exact diagonalization, but which can be used to treat much larger systems. Being able to treat sufficiently large systems is important since the models we study here have two fermionic sites, i.e. sixteen degrees of freedom per unit cell, which would severely limit the maximum size available for exact diagonalization. Quantum Monte Carlo methods would suffer from the fermion sign problem at all band fillings on the $`d`$$`p`$ lattice.

We use the finite–system version of the ground–state DMRG algorithm  to accurately calculate the energy expectation values of equal–time operators in the ground state including local spin, density, and various correlation functions. We keep up to 800 states in the system block and treat lattices of up to $`N=32`$ unit cells. The maximum sum of discarded density matrix eigenvalues is approximately $`5\times 10^6`$.

In order to search for ferromagnetic phases, it is essential to be able to determine the total spin of the ground state. Since our version of the DMRG algorithm does not allow direct control of the total spin, we use a combination of methods to do this. First, we can calculate the expectation value $`\psi _0|S^2|\psi _0`$ directly, where $`\psi _0`$ is the variational DMRG ground state calculated in a particular $`S_z`$ sector, and $`S^2`$ is the explicit operator for the square of the total spin. Since the DMRG ground state is variational, $`\psi _0|S^2|\psi _0`$ does not always take on the quantized values due to mixing between states of different $`S_z`$. Second, we can add a term of the form $`\lambda S^2`$ to the Hamiltonian, raising the energy of higher spin states within a particular $`S_z`$ sector. We can use this to calculate the ground–state energy as a function of $`S^2`$ by examining the energy of the minimum $`S^2`$ ground state in a particular $`S_z`$ sector. Finally, we can examine the ground–state energies as a function of $`S_z`$ and use the degeneracy in $`S_z`$ to determine $`S^2`$. For application of these methods to the ferromagnetism in the PAM, see Ref.  and to the Hubbard chain with next–nearest–neighbor hopping, see Ref. .

### A  Comparison with the periodic Anderson model

In the first part of this work, we present a comparison between the $`d`$$`p`$ system and the PAM. The qualitative picture developed for the ferromagnetic phase of the PAM  should also apply to the $`d`$$`p`$ model. In order to make the correspondence to the PAM, we choose a regime in which the $`d`$-orbital lies below the bottom of the $`p`$-band. At exactly quarter filling, all the electrons are in the $`d`$-sites and all the spin states are degenerate if there is no hybridization. When the hybridization is turned on, the super-exchange favors antiferromagnetic interactions. The exchange coupling constant can be calculated from perturbation theory giving

$$J_A=\frac{4t_{pd}^4}{(\mathrm{\Delta }+V_{pd})^2}\left(\frac{1}{U_d}+\frac{2}{2\mathrm{\Delta }+U_p}\right),$$
(4)
or when $`U_d=2\mathrm{\Delta }`$ and $`U_p=V_{pd}=0`$, $`J_A=6t_{pd}^4/\mathrm{\Delta }^3`$. When a hole is added to the quarter–filled system, it will prefer to go to the $`p`$-sites since $`U_d>|\mathrm{\Delta }|`$. This extra hole will tend to form singlets with the mostly localized spins in the d-sites. The hole can delocalize and thus lower its energy if the spins of the $`d`$–holes are oriented in the same direction, i.e. are ferromagnetically ordered. The binding energy for these singlets can be calculated within perturbation theory in a manner similar to that applied to the two-dimensional case by Zhang and Rice,  yielding

$$J_S=2t_{pd}^2\left(\frac{1}{\mathrm{\Delta }V_{pd}}+\frac{1}{\mathrm{\Delta }V_{pd}+U_p}+\frac{2}{U_d\mathrm{\Delta }V_{pd}}\right).$$
(5)
For $`U_d=2\mathrm{\Delta }`$ and $`U_p=V_{pd}=0`$, $`J_S=8t_{pd}^2/\mathrm{\Delta }`$.

In general, when there are $`N_p=N_hN`$ holes in the $`p`$-band, this effect will favor a ferromagnetic ground state with total spin $`S=(NN_p)/2`$. We will denote this value of $`S`$ for the ground state as complete ferromagnetism, while a value of $`S`$ smaller but still greater than the minimum, will be labeled incomplete. Note that a complete state has a lower magnetization than a saturated one. It represents a state where all the uncompensated d-spins are aligned.

Thus, one would expect a competition between ferromagnetism and antiferromagnetism near quarter filling similar to that found in the PAM. For the PAM there is an exact result that proves the ground state is ferromagnetic with complete magnetization for $`N_c=1`$ ($`N_e=N+1`$) and $`U_d=\mathrm{}`$. There is no such result, however, for the $`d`$$`p`$ model. Furthermore, while the antiferromagnetic exchange is a sixth–order process in the Anderson lattice model, it is fourth order (and therefore stronger) in the $`d`$$`p`$ system. Therefore, it is not clear that the ferromagnetism will still be present in the $`d`$$`p`$ system. Yanagisawa  studied a $`2\times 2`$ CuO cluster with exact diagonalization and found a ferromagnetic ground state. However, it is difficult to draw general conclusions about the behavior in the thermodynamic limit from such a small cluster.

In our DMRG calculations, we consider chains with $`N=16`$ and open boundary conditions and set $`t_{pp}=0.5`$, $`t_{pd}=0.375`$, $`\mathrm{\Delta }=U_d/2`$ and $`U_p=V_{pd}=0`$. This choice allows us to compare the phase diagram with the one of the Anderson lattice model from Ref.  where $`t=0.5`$, $`V=0.375`$, $`ϵ_f=U/2`$. In Fig. 2, we present both phase diagrams in the plane of band filling and the local Coulomb repulsion, $`U`$ or $`U_d`$. Here $`n_c=0`$ ($`n_p=0`$ for the $`d`$$`p`$ model) represents the quarter-filled case, which in the large Coulomb repulsion regime means that the localized orbitals are singly occupied and the extended orbitals are empty in the zero hybridization limit. In the $`d`$$`p`$ system, Fig. 2(b), the region of complete ferromagnetism is pushed towards higher values of $`U_d`$ and higher densities and there is a wide region of incomplete states at low densities. This is consistent with the fact that antiferromagnetic correlations at quarter filling in this system are of lower order (and therefore stronger) than in the Anderson lattice model.

In Fig. 3 we show the Cu–Cu correlation function, $`S_d^+(r)S_d^{}(0)`$, as a function of distance, $`r`$ (measured in units of the lattice constant), for the $`d`$$`p`$ model with the parameters described above and $`U_d=5`$ for $`N_p=0`$, 2 and 4. Here $`S_d^+(r)=d_{r,}^{}d_{r,}`$ is the $`d`$–spin raising operator on site $`r`$. For $`N_p=0`$ (quarter filling), the correlations are clearly antiferromagnetic with the amplitude decaying slowly with distance, while for $`N_p=2`$ and $`N_p=4`$, the correlations are ferromagnetic. This illustrates that the system goes from an antiferromagnetic state at quarter filling to a ferromagnetic state as soon as additional holes are added for these parameter values.

The competition between the tendency to ferromagnetism and the antiferromagnetic exchange gives rise to the region of incomplete ferromagnetism at low doping. In this region, we find that the system phase separates into domains of complete ferromagnetism in which the $`p`$–holes are localized and antiferromagnetic domains in which the $`p`$–band is empty. Indeed, the first symptom of this phase separated state can be seen in Fig. 3, in which there is a jump in the correlation function for $`r=12`$ for the case of 18 holes. This peculiar behavior is due to the tendency of the holes to localize and form a region of complete ferromagnetism near the center of the chain, leaving the p-orbitals empty near the ends.

To illustrate this phase–separated regime, we have calculated the chemical potential $`\mu =E_o(N_h+1)E_o(N_h)`$, where $`E_o(N_h)`$ is the ground state energy with $`N_h`$ holes, for chains of $`N=8`$, 16 and 32 unit cells. We show $`\mu `$ as a function of the filling for $`U_d=5`$ in Fig. 4. For 8 unit cells, $`\mu `$ increases monotonically with $`n_p`$, but for 16 unit cells there is a shallow minimum at $`n_p=0.0625`$. For 32 unit cells, there is a flat region which is shown in detail in the inset, with a minimum at $`n_p=0.09375`$ which is indicative of phase separation. We find a similar behavior of the chemical potential is found for $`U_d=6`$ and $`U_d=4`$.

To illustrate the localization of the holes in the $`p`$–band, we display the density in the $`p`$–sites as a function of position in the lattice for $`U_d=4`$ in an $`N=32`$ unit–cell chain with $`N_p=0,1,2,3`$ in Fig. 5. At quarter filling, the $`p`$–levels are almost empty. Upon doping, they begin to fill with holes, but instead of spreading out over the lattice, the additional holes localize in a region which is smaller than the lattice size. This behavior is also seen for $`U_d=5`$ and $`U_d=6`$.

In Fig. 6 we show the local spin density $`S_i^z`$ as a function of position in the lattice for an $`N=32`$ chain with $`U_d=5`$, $`N_p=2`$ and $`S_z=S=5`$ (the ground state) for both the $`p`$–band and the $`d`$–band. There are two distinct domains in the lattice: one in which the $`d`$–holes are completely polarized and another one in which they tend to order antiferromagnetically. The coupling between $`p`$ and $`d`$ holes is clearly antiferromagnetic. The asymmetry of the distribution is a numerical artifact due to the near–degeneracy of states with the phase–separated region centered at different points.

To gain some physical insight into the phase separation mechanism, one can use energetic arguments to estimate the stability of competing phases. The energy scale in the antiferromagnetic state is set by $`J_A`$, the antiferromagnetic exchange energy, Eq. (4). The energy of the uniform ferromagnetic phase can be estimated by the binding energy of the Zhang–Rice–like singlets, $`J_SN_p`$. For a rough estimate of the magnetic energy scale, consider an antiferromagnetically ordered state versus one that has complete ferromagnetism. The Bethe Ansatz energy of a one–dimensional Heisenberg antiferromagnet  consisting of $`N1`$ bonds (measured relative to the ferromagnetic state and neglecting end effects) is $`J_A^{\mathrm{eff}}(N1)`$, where $`J_A^{\mathrm{eff}}=\mathrm{ln}2J_A`$. In order for complete ferromagnetism to occur, the condensation energy of the singlets, $`J_SN_p`$ must be larger than this, leading to the condition

$$n_p\mathrm{ln}2J_A/J_Sn_{pc1}.$$
(6)
This estimate is only valid for very low densities of $`p`$–electrons, since it neglects their kinetic energy, and for strong coupling where the magnetic picture is valid. According to this estimate, ferromagnetism is favored over antiferromagnetism for $`U_d3`$ at all densities above quarter-filling on the lattice sizes that we have studied. However, it yields a smaller $`n_{pc1}`$ than we find in the numerical calculations. For example, for $`U_d=4`$ and the parameter set used in this section, $`n_{pc1}=0.018`$.

Now let us examine the possibility of phase separation between a region of undoped antiferromagnetism and a region of complete ferromagnetism. In order for such a phase to exist, the energy gained by forming an undoped antiferromagnetic domain must compensate the loss in kinetic energy of the Zhang–Rice singlets, which have to localize in a region smaller than the lattice size $`N`$. We estimate the kinetic energy of the Zhang–Rice singlets by taking them to be noninteracting hard–core bosons confined to a lattice of length $`X`$ with open boundary conditions. (We study lattices with open boundary conditions here, so that this form is valid up to the number of $`p`$–sites, $`X=N+1`$.) The single–particle energy levels are given by

$$\epsilon _j=2t_{ZR}\mathrm{cos}\frac{\pi j}{X+1},$$
(7)
where $`t_{ZR}`$ is the effective hopping of the singlets and $`j=1,\mathrm{},N`$. Since hard–core bosons are equivalent to spinless fermions in one dimension, the total kinetic energy is given by successively occupying single–particle states,

$`T(N_p,X)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N_p}{}}}\epsilon _j`$ (8)
$``$ $`{\displaystyle _0^{N_p}}𝑑j\epsilon _j2t_{ZR}{\displaystyle \frac{X}{\pi }}\mathrm{sin}{\displaystyle \frac{N_p\pi }{X}},`$ (9)
where the approximation on the right–hand side is valid in the continuum limit: $`X1`$. The energy gained by constricting the ferromagnetic domain to a length $`X`$ is

$`E_{PS}(N_p,X)=`$ $`J_A^{\mathrm{eff}}(N1X)`$ (10)
$`+`$ $`T(N_p,N+1)T(N_p,X).`$ (11)
The phase–separated state will be favored over the homogeneous ferromagnetic state if $`E_{PS}`$ is positive. In that case, the size of the ferromagnetic region is given by the value of $`X`$ that maximizes the energy gain, determined by

$$0=\frac{E_{PS}}{X}J_A^{\mathrm{eff}}+\frac{2\pi ^2t_{ZR}N_p^3}{3X^3},$$
(12)
where we have taken $`N_p\pi 1`$ in order to expand the sine and have used the continuum expression from the right–hand side of Eq. (9). Eq. (11) has one maximum for $`X>0`$ which occurs at

$$X_{\mathrm{max}}=N_p\left(\frac{2\pi ^2t_{ZR}}{3J_A^{\mathrm{eff}}}\right)^{1/3}N_px_o.$$
(13)
It is easy to verify that $`E_{PS}(N_p,X_{\mathrm{max}})`$ is always positive as long as $`0<X_{\mathrm{max}}<N`$. Therefore $`X_{\mathrm{max}}`$ is the length of the ferromagnetic domain when phase separation is present. Note that $`X_{\mathrm{max}}N_p`$, in agreement with the behavior seen in Figs. 5 and 6. In order to numerically estimate $`x_0`$ for the $`d`$$`p`$ system, we take $`t_{ZR}=t_{pp}`$ and $`J_A^{\mathrm{eff}}=\mathrm{ln}2J_A`$, where $`J_A`$ is given by Eq. (4) for the $`d`$$`p`$ model. (One could calculate $`t_{ZR}`$ more accurately in perturbation theory, but only its order of magnitude is important, since it appears within the cube root.) For the parameters used in this section, this gives $`x_0=7`$, 8.5 and 10 for $`U_d=4`$, 5 and 6, respectively. This is in reasonable agreement with our numerical results (See Fig. 5), which also show an increase in the localization length with $`U_d`$.

The boundary between the uniform ferromagnetic phase and the phase–separated phase is determined by the condition $`X_{max}=N`$, leading to a critical density

$$n_{pc2}=\left(\frac{3J_A^{\mathrm{eff}}}{2\pi ^2t_{ZR}}\right)^{1/3}.$$
(14)
The densities $`n_{pc1}`$ and $`n_{pc2}`$ calculated from Eqs. (6) and (14) are represented in the phase diagram of Fig. 2(b) as dotted and dashed lines, respectively. Notice that the numerical calculation yields a somewhat wider region of phase separation than given by these estimates. In order to estimate the finite–size effects, we have also examined the spin of the ground state on $`N=8`$ and $`N=32`$ unit–cell chains. We find that the points of complete polarization in Fig. 2(b) behave consistently with system size, i.e. are also completely polarized for $`N=8`$ and $`N=32`$. However, the scaling behavior of the incomplete states is more complicated. In particular, some states which are incompletely polarized in the $`N=16`$ chain become completely polarized for $`N=32`$. These cases are circled with a dashed line in Fig. 2(b). Therefore, the phase–separated region tends to become narrower in the thermodynamic limit. There is, however, evidence that some of the incompletely polarized states remain in the thermodynamic limit. For $`U_d=5`$ and 6, the two densities closest to quarter filling on the $`N=16`$ chain have total spin $`S=5/2`$ and $`5`$, and $`S=7/2`$ and $`6`$, respectively. On the $`N=32`$ chain, the ground–state $`S`$ at the same density doubles, suggesting that $`S/N`$ will remain constant in the thermodynamic limit. At $`U_d=4`$ and $`U_d=3`$, none of the incomplete states show such a simple scaling and we believe that larger system sizes are needed to determine the nature of the phase boundaries. For $`U_d=3`$, we find no evidence of phase separation at any finite system size.

Let us now examine the relevance of the phase separation arguments to the PAM. Our energetic arguments should also apply with the appropriate mapping of $`J_A`$ and $`J_S`$ to the model parameters. For the symmetric case, $`\epsilon _f=U/2`$, $`J_S=8V^2/U`$ and $`J_A=8V^4t^2/\epsilon _f^5`$. For the parameter sets we have compared on the two models, $`J_A`$ is smaller in the PAM than in the $`d`$$`p`$ system. This implies that phase separation occurs over much larger localization lengths in the PAM. Indeed, for the conduction electron density $`n_c=1/32`$ we obtain $`x_o28`$, $`21`$ and $`14`$ for $`U=6`$, $`5`$ and $`4`$, respectively. Because these lengths are larger than those for the $`dp`$ system, we could not find definite numerical evidence of a phase–separated state on the lattice sizes that we were able to study; larger systems sizes would be required. In order to confirm our picture of the origin of the phase separation, we can instead reduce $`x_0`$ by increasing the hybridization $`V`$; we double $`V`$ to $`V=0.75`$, decreasing the localization length by a factor of two to $`x_o14`$, $`10`$ and $`7`$ for $`U=6`$, $`5`$ and $`4`$, respectively. Using the DMRG, we then examine $`N=32`$ systems with $`t=0.5`$, $`ϵ_f=U/2`$ and $`V=0.75`$ for $`U=4`$, 5, and 6. The results, presented in Fig 7, show numerical evidence for phase separation for all three $`U`$–values. The chemical potential, Fig. 7(a), has a minimum as a function of the filling, and the conduction electron density, Fig. 7(b), tends to localize in a region smaller than the size of the system. Therefore, we find that phase separation can also occur in the PAM under the right conditions.

### B  Phase Diagram for “Realistic” Parameters

In the second part of our work, we think of the $`d`$$`p`$ system as the one–dimensional analog of the copper–oxide planes in high temperature superconductors. We set the parameters $`t_{pd}=1`$, $`U_d=8`$, $`U_p=4`$, $`\mathrm{\Delta }=3`$ and $`V_{pd}=0.5`$, values reasonable for the cuprates, and then vary $`t_{pp}`$ and the hole filling. The resulting phase diagram is shown in Fig. 8 for an $`N=16`$ chain and both possible signs of $`t_{pp}`$. The cases with complete polarization, $`S=(NN_p)/2`$, are circled. For $`t_{pp}>0`$ the phase diagram has one ferromagnetic region consisting almost exclusively of complete polarization for $`t_{pp}>1`$ and from $`n_p0.4`$ to $`n_p0.75`$. The arguments of the first part of our work should also apply to this region, i.e. the ferromagnetic ordering of the $`f`$–spins is caused by the condensation of Zhang–Rice singlets. As can be seen from the lack of regions of incomplete polarization, we find no phase separation at these parameter values. It is clear that Eq. (6) cannot properly predict the boundary of the ferromagnetic phase since it is invalid at the relatively large $`N_p`$ values at which we find the phase. From Eq. (13), the charactistic size of the localization of the singlets is given by $`x_04.5t_{pp}^{1/3}`$ for these parameter values. Our energetic arguments therefore predict no phase separation at the densities at which we do find a ferromagnetic phase.

When $`t_{pp}<0`$ the phase diagram is quite different. There is one point at quarter filling and $`t_{pp}=2`$ for which the ground state is fully polarized ($`S=8`$ for $`N=16`$), surrounded by a region of smaller polarization. Additionally, there is a region with polarization greater than complete, but less than fully polarized for a narrow range of $`t_{pp}`$ but a wide range of hole occupation, $`N_p`$. More precisely, the extent of this region is $`0.5t_{pp}0.7`$ and $`7N_p14`$.

We first discuss the fully polarized phase at quarter filling and $`t_{pp}2.0`$. One can gain insight into the origin of this ferromagnetic phase by examining the noninteracting band structure. Diagonalization of the Hamiltonian, Eq. (2), with $`U_d=U_p=V_{pd}=0`$ and periodic boundary conditions leads to the hybridized bands

$$\epsilon _\pm (k)=\mathrm{\Delta }/2t_{pp}\mathrm{cos}k\pm \sqrt{(\mathrm{\Delta }/2t_{pp}\mathrm{cos}k)^2+2t_{pd}^2(1+\mathrm{cos}k)}.$$
(15)
These bands, plotted for three $`t_{pp}`$ values near the fully polarized point, are shown in Fig. 9. As can be seen, the lower band does not overlap with the upper band and has a dispersion that is quite flat; it becomes completely flat at $`t_{pp}=1.78`$. Therefore, all the holes will go into the half–filled flat band at this point. A theorem by Mielke,  subsequently extended to nearly flat bands by Mielke and Tasaksi  shows rigouously that for the Hubbard model with a low–lying completely flat band at half–filling or less will have a fully polarized ground state. Therefore, this region is due to flat–band ferromagnetism. This is substantiated by the strong dependence of the polarization on $`t_{pp}`$ and $`n_p`$; there is only a small region of partial polarization surrounding the fully polarized point.

We next treat the narrow strip of ferromagnetism at $`t_{pp}0.5`$ and intermediate $`N_p`$. The first important issue is whether the partial polarization in this region persists in the thermodynamic limit. In Fig. 10, we plot $`S/N`$ as a function of $`1/N`$ for three different fillings. As can be seen, $`S/N`$ increase linearly with $`1/N`$ and extrapolates to a finite value as $`N\mathrm{}`$.

One can gain some insight into the relationship between the sign of $`t_{pp}`$ and the mechanism for ferromagnetism by examining the behavior of small clusters. In particular, we examine the behavior of one triangular element of the lattice, i.e. a cluster consisting of one copper and two oxygen sites, occupied by two holes ($`N=1`$, $`N_p=1`$). This CuO2 cluster can be treated analytically if spin and reflection symmetry are taken into account: The total spin can either be zero or one, and states can be symmetric and antisymmetric with respect to exchange of the oxygens. The antisymmetric $`S=0`$ and the symmetric $`S=0`$ sectors can be represented by $`2\times 2`$ matrices, the symmetric $`S=0`$ sector is $`4\times 4`$, and the antisymmetric $`S=1`$ sector is $`1\times 1`$. All of these matrices can be diagonalized analytically (although it is easier, in practice, to treat the $`4\times 4`$ matrix numerically since the general solution is rather complicated). Taking $`t_{pd}=1`$, $`U_d=8`$, $`U_p=4`$, $`\mathrm{\Delta }=3`$ and $`V_{pd}=0.5`$, and varying $`t_{pp}`$, we find that for $`t_{pp}>0`$, the ground state is always a singlet with positive reflection symmetry. For $`t_{pp}<0`$, there is a competition between singlet formation, favored by direct antiferromagnetic exchange, and triplet formation, favored by third order exchange around the plaquette, which is ferromagnetic.

In Fig. 11, we show the dependence of the energy of the three lowest lying states of the CuO2 cluster on $`t_{pp}`$. As can be seen, for $`0<t_{pp}0.46`$ and $`t_{pp}4.5`$, the ground state is a singlet with positive reflection symmetry. In the intermediate range, $`0.46t_{pp}4.5`$, the ground state is an antisymmetric triplet, which is close in energy to the first excited state, an antisymmetric singlet. This illustrates the effect of the sign of $`t_{pp}`$ on the basic triangular element of the lattice: a localized triplet is only possible for negative $`t_{pp}`$ and can occur only for a certain range of $`t_{pp}`$, when there are two holes on the plaquette.

While the triplet ground state of the CuO2 cluster motivates the possibility of a ferromagnetic state in this region, and provides an explanation for the minimum value of $`t_{pp}`$ required, it cannot provide a complete picture of this phase. In particular, the range of $`t_{pp}`$ with a triplet ground state is much wider in the CuO2 cluster than the range of $`t_{pp}`$ in the ferromagnetic strip in Fig. 8. When more triangles are connected together, interaction between holes in the $`p`$–band can become important. The point of maximum polarization occurs when $`n_p0.5`$, i.e. there is one hole for every second $`p`$–site. At this filling, ferromagnetism based on exchange around the triangular elements would be somewhat frustrated since each $`p`$–site is shared by two triangular elements. These mechanisms would tend to suppress ferromagnetism. Also, as seen in Fig. 11, the triplet ground state is close in energy to a singlet excited state even on a single triangular element.

We can gain additional insight into these effects from diagonalization on clusters larger than three sites. For clusters of $`N=2`$, 3, and 4 unit cells (5, 7, and 9 total sites) with open boundary conditions, a ferromagnetic ground state appears at all fillings in the range $`0<N_pN`$ for $`t_{pp}t_{pp}^{\mathrm{min}}`$, with $`t_{pp}^{\mathrm{min}}`$ ranging from -0.52 to -0.48. However, the extent of this region in $`t_{pp}`$ and the value of the polarization depends on $`N_p`$. For $`N_p=1`$, one hole in the $`p`$–band, and all three sizes, we obtain a ground–state spin at or close to the maximum possible value, up to $`t_{pp}^{\mathrm{max}}`$ ranging between -2.9 and -4.2. This behavior is like that in the $`N=1`$ CuO2 cluster, since all the $`d`$–holes are ferromagnetically polarized. It is also reminiscent of the ferromagnet phase in the PAM, except that it is a triplet state which delocalizes rather than a singlet. Similar ferromagnetic states have been found in the one-dimensional Kondo lattice model with ferromagnetic exchange.

For larger $`N_p`$, we obtain a ferromagnetic ground state up to $`t_{pp}^{\mathrm{max}}`$ ranging from -0.62 to -0.70, a value in good agreement with that found in larger systems. The ground–state spin varies and depends on $`N_p`$. One can understand the nature of this phase by examining the distribution of spin polarization and the spin–spin correlation function. We display the local spin polarization, $`S_i^z`$, for parameters that yield the maximum polarization on a $`N=16`$ cluster in Fig. 12, taking $`S_z=5`$, the maximal value. One can see that the holes on the $`d`$–sites are antiferromagnetically aligned, which we have confirmed by examining the spin–spin correlation function. The $`p`$–sites, however, all have positive polarization, and are ferromagnetically correlated with the positively polarized $`d`$–sites, i.e. with every other $`d`$–site. The correlations between $`p`$–sites are generally ferromagnetic. The picture is then that there is an alternation between triangular plaquettes with a ferromagnetic moment and isolated $`d`$–holes antiferromagnetically aligned with the neighboring plaquettes. This phase is reminiscent of a ferrimagnet with alternating $`S=1/2`$ and $`S=1`$. Changing the band filling can frustrate this order, changing the amount of polarization. One domain wall is visible in Fig. 12. Since this phase is due to a subtle competition between ferromagnetic and antiferromagnetic exchange processes, it is not surprising that it is quite sensitive to the value of $`t_{pp}`$. Certainly the picture of localized plaquettes breaks down as the hopping between the $`p`$–sites increases.

## IV  Conclusions

We have studied ferromagnetism in the phase diagram of the one dimensional $`d`$$`p`$ model. In the first part of this work, we have compared the phase diagram near quarter filling of a $`d`$-$`p`$ system to the one-dimensional periodic Anderson model. In both cases, we have found a ferromagnetic region near quarter filling, but, for equivalent values of the parameters, the ferromagnetic region in the $`d`$$`p`$ system is pushed towards higher values of the Coulomb interaction and higher densities.

Between the antiferromagnetic state at quarter filling and the homogeneous ferromagnetic region, there is a phase-separated regime in which ferromagnetic domains are segregated from the antiferromagnetic background. The phase–separated state arises because the loss of kinetic energy of the holes due to their localization in a small region is compensated by the gain in antiferromagnetic bonds. This regime is more evident in the $`d`$$`p`$ system because the antiferromagnetic coupling is stronger than in the Anderson model and therefore the characteristic length for the ferromagnetic domains is smaller. Nevertheless, we have presented numerical evidence for the phase–separated state for both models.

Recently, similar phases with a hole-rich ferromagnetic domains segregating from hole-undoped antiferromagnetic regions have been reported for models such as the one-orbital ferromagnetic Kondo model and other models studied in relation with manganites. (See Ref.  and references therein.) In our case, it must be emphasized that the ferromagnetism is not saturated as in the case of these models. In the case of the ferromagnetic Kondo model, the phase–separated state was found near half-filling. However, it was pointed out that such a phase–separated state can arise at low densities if an antiferromagnetic Heisenberg coupling is added to the model. Several studies in different numbers of dimensions using different methods suggest that the phase–separated state is very robust. It is also mentioned that the inclusion of longer–range Coulomb repulsion can turn the phase–separated state into more exotic phases, such as striped phases with ferromagnetic domains. While we have only studied a one-dimensional system, the arguments presented for the physical mechanism of phase separation should still be valid in higher dimensionsal systems.

Also, it is important to notice that the Anderson model can be mapped to the antiferromagnetic Kondo model for large values of the Coulomb repulsion. It is well known that a ferromagnetic phase is also present in the antiferromagnetic Kondo lattice model. However, there is no phase–separated state because the antiferromagnetic coupling is higher order than the mapping. We believe that an additional antiferromagnetic Heisenberg coupling of appropriate strength between the localized spins would induce a phase–separated state in the Kondo lattice model.

In the second part of our work, we have studied the $`d`$$`p`$ model the parameters to values relevant to the cuprates. We have found that the sign of the oxygen hopping $`t_{pp}`$ is crucial in determining the physics of the ground state. Three disconnected ferromagnetic phases were found in the phase diagram. One phase occurs at positive $`t_{pp}`$ and is similar to the ferromagnetic phase found in the PAM and discussed in the first part of this work. A second phase appears for negative $`t_{pp}`$ at and near quarter filling and occurs at a point at which the occupied noninteracting band has a flat dispersion. This “flat-band” mechanism is in agreement with a picture based on exact theorems. A third ferromagnetic phase appears for small negative $`t_{pp}`$ and it is related to ferromagnetic electron exchange processes on plaquettes with a positive product of hopping integrals. This phase has elements of ferrimagnetism due to an ordering of triangular plaquettes with a ferromagnetic moment and isolated, antiferromagnetically aligned $`d`$–holes.

## Acknowledgments

We thank S. Trugman and S. Daul for helpful discussions. M. Guerrero acknowledges the support of the U. S. Department of Energy. R.M.N acknowledges an allocation of computer time from the Centro Svizzero di Calcolo Scientifico in Manno and support from the Swiss National Foundation under Grant Nos. 20–46918.96 and 20–53800.98.


text:
# Can the Universe escape eternal acceleration?

Recent astronomical observations of distant supernovae light-curves, -, suggest that the expansion of the universe has recently begun to accelerate. Acceleration is created by an anti-gravitational repulsive stress, like that produced by a positive cosmological constant, or universal vacuum energy. It creates a rather bleak eschatological picture. An ever-expanding universe’s future appears to be increasingly dominated by its constant vacuum energy. A universe doomed to accelerate forever will produce a state of growing uniformity and cosmic loneliness. Structures participating in the cosmological expansion will ultimately leave each others’ horizons and information-processing must eventually die out. Here, we examine whether this picture is the only interpretation of the observations. We find that in many well-motivated scenarios the observed spell of vacuum domination is only a transient phenomenon. Soon after acceleration starts, the vacuum energy’s anti-gravitational properties are reversed, and a matter-dominated decelerating cosmic expansion resumes. Thus, contrary to general expectations, once an accelerating universe does not mean always an accelerating universe.

The observed cosmic acceleration, if due to vacuum repulsive stresses, is very puzzling. In contrast to conventional forms of matter, like pressureless dust or radiation, the cosmological energy densities contributed by these repulsive stresses are not diluted by expansion: they remain constant. Consequently, the vacuum energy will eventually exceed the densities of all other forms of matter in an ever-expanding universe, causing its expansion to accelerate. Yet, in order to explain why the vacuum energy comes to dominate other matter only close to the present time, $`15Gyr`$ after the expansion began, we must assume it starts out with a density that is $`10^{120}`$ times smaller than that of other forms of matter, for no known reason. A possible explanation is the ’quintessence’ picture, in which a scalar matter field is present in the universe and can display repulsive gravitational behaviour late in the universe’s history, despite remaining innocuous for most of its life . In the most extreme case, its influence can become identical to the presence of a positive cosmological constant with constant energy density but the possible time-variation permitted in the evolution of its density allows acceleration to appear at late times following more natural initial conditions. We shall show that in well-motivated scenarios of this type the observed spell of vacuum domination is only a transient phenomenon.

We consider a homogeneous and isotropic universe with zero spatial curvature that contains two dominant forms of matter: a perfect fluid with pressure $`p`$ and density $`\rho `$ linked by an equation of state $`p=(\gamma 1)\rho `$, with $`\gamma `$ constant, together with a hypothetical scalar ’quintessence’ field $`\varphi `$ defined by its self-interaction potential $`V(\varphi )`$. All variables can depend only on the cosmic co-moving proper time $`t`$. The cosmological expansion scale factor is $`a(t)`$ and the Hubble expansion rate is defined by $`H\dot{a}/a.`$ Separate conservation of mass-energy for the two matter fields as the universe expands requires $`\rho a^{3\gamma }`$, while the scalar field, $`\varphi `$, obeys

$$\ddot{\varphi }+3H\dot{\varphi }=\frac{V}{\varphi },$$
(1)
where overdots represent time derivatives. The densities of the perfect fluid, $`\rho `$, and the $`\varphi `$ field, $`\rho _\varphi `$, drive the expansion, so

$$3H^2=\rho +\rho _\varphi ,$$
(2)
where the scalar-field’s energy density is the sum of its kinetic and potential parts, $`\rho _\varphi \frac{1}{2}\dot{\varphi }^2+V.`$ We use Planck units defined by $`\mathrm{}=c=(8\pi G)^{1/2}=1.`$

If we choose a monotonic potential of the form $`V(\varphi )=Ae^{\lambda \varphi }`$, with $`A`$ and $`\lambda >\sqrt{2}`$ constant parameters, then $`V`$ need not increasingly dominate the right-hand side of eq. (2) at large $`a`$. Instead, $`\rho _\varphi `$ evolves in direct proportion to $`\rho `$ $`t^2`$ during an era of radiation ($`\gamma =4/3`$) or matter ($`\gamma =1`$) domination. Hence $`H=2/3\gamma t`$ and the fraction of the total density in the form of quintessence during this ’scaling’ regime depends only on $`\lambda :`$

$$\mathrm{\Omega }_\varphi \frac{\rho _\varphi }{\rho _\varphi +\rho }=\frac{\rho _\varphi }{3H^2}=\frac{3\gamma }{\lambda ^2}.$$

This unusual type of behaviour shows how vacuum domination can be postponed for a long period since the universe’s expansion began. However, eventually, the expansion departs from this special scaling behaviour and the potential of the $`\varphi `$ field comes to dominate eq. (2), contributing an almost constant value of $`\rho _\varphi `$ and accelerating expansion. Eventually, the potential energy $`V`$ supersedes the kinetic energy,$`\frac{1}{2}\dot{\varphi }^2,`$ making quintessence behave like an effective cosmological constant. In the standard scenarios the early scaling behaviour is a robust feature of quintessence models, but subsequent late-time acceleration can only be achieved with a large degree of fine tuning of the defining parameters - either in the initial conditions or in the parameters of the potential.

Albrecht and Skordis  have proposed a particularly attractive model of quintessence. It is driven by a potential which introduces a small minimum to the exponential potential:

$$V(\varphi )=e^{\lambda \varphi }\left(A+(\varphi \varphi _0)^2\right).$$
(3)
Unlike previous quintessence models, late-time acceleration is achieved without fine tuning of the initial conditions. The authors argue that such potentials arise naturally in the low-energy limit of $`M`$-theory; the constant parameters, $`A`$ and $`\varphi _0`$, in the potential take values of order $`1`$ in Planck units, so there is also no fine tuning of the potential. They show that, regardless of the initial conditions, $`\rho _\varphi `$ scales, with $`\rho \rho _\varphi t^2`$ during the radiation and matter eras, but leads to permanent vacuum domination and accelerated expansion after a time which can be close to the present. Acceleration begins when the field gets trapped in the local minimum of the potential at $`\varphi =\varphi _0+(1\pm \sqrt{1\lambda ^2A})/\lambda `$, which is created by the quadratic factor in eq. (3) when $`1\lambda ^2A`$. Once the field gets stuck in the false vacuum its kinetic energy disappears ($`\varphi `$ constant), and the ensuing dominance of $`\rho +\rho _\varphi `$ by an almost constant value of the potential value triggers a period of accelerated expansion that never ends. The probability for quantum tunnelling through the barrier is negligible .

We have found that this type of behaviour is by no means generic. In Fig. 1 we plot the evolution of $`\mathrm{\Omega }_\varphi `$ and a measure of the effective equation of state of the total matter content of the universe $`w_{tot}=p_{tot}/\rho _{tot}`$, where the total density is $`\rho _{tot}=`$ $`\rho +\rho _\varphi `$ and the total pressure is $`p_{tot}=(\gamma 1)\rho +\frac{1}{2}\dot{\varphi }^2V`$. We show $`\mathrm{\Omega }_\varphi `$ and $`w_{tot}`$ in two qualitatively distinct cases, created by equally plausible values for potential parameters. Vacuum domination occurs when $`\mathrm{\Omega }_\varphi >1/2`$, and accelerated expansion when $`w_{tot}<1/3`$. The dashed line traces the case found by Albrecht and Skordis, but the solid line shows the general behaviour. In both cases there is scaling of the densities early in the expansion history: with $`w_\varphi =1/3`$ and $`\mathrm{\Omega }_\varphi 4/\lambda ^2`$ in the radiation era, followed by $`w_\varphi =0`$ and $`\mathrm{\Omega }_\varphi 3/\lambda ^2`$ in the matter era . The expansion decelerates in both scaling eras. But, then, in both cases, as $`\varphi \varphi _0`$, vacuum domination and accelerated expansion is triggered. Whereas in the first case this phenomenon is permanent, in the latter it is ephemeral. In the latter case, $`\varphi `$ soon continues rolling down the potential, and the universe resumes scaling evolution with $`\rho `$ $`\rho _\varphi `$ $`t^2`$ and another matter-dominated era ensues. The spectre of never-ending vacuum domination has been lifted.

Transient vacuum domination arises in two ways. When $`A\lambda ^2<1,`$ the $`\varphi `$ field arrives at the local minimum with enough kinetic energy to roll over the barrier and resume descending the exponential part of the potential where $`\varphi >>\varphi _0`$. This kinetic energy is determined by the scaling regime, and so by parameters of the potential and not initial conditions. Another instance of transient vacuum domination is the whole region $`A\lambda ^2>1`$. As $`A`$ increases towards $`\lambda ^2`$, the potential loses its local minimum, and flattens out into a point of inflexion. This is sufficient to trigger accelerated expansion temporarily, but the field never stops rolling down the potential, and matter-dominated scaling evolution with $`a(t)t^{2/3}`$ is soon resumed.

Further scrutiny of the theory reveals a critical dividing line separating permanent and transient vacuum domination. If we impose on the parameter space $`\{A,\varphi _0,\lambda \}`$ the condition that accelerated expansion be occurring by the recent past, as observations imply, then we require $`\lambda \varphi _0280`$. This leaves us with two degrees of freedom, which we parametrise by $`\lambda `$ and $`A\lambda ^2`$. The behaviour of the Universe in this parameter space is shown in Fig. 2. The unshaded area represents the regions of permanent vacuum domination. Outside this area we have plotted contours giving the number of e-foldings of accelerated expansion, defined by $`N=log(a_f/a_i)`$, where $`a_f`$ and $`a_i`$ are the values of the expansion scale factor at the end and start of acceleration, respectively.

Our investigations have revealed a new type of cosmological evolution. It is possible for the universe to exit from a period of accelerated expansion and resume decelerated expansion. Moreover, for the well-motivated family of Albrecht-Skordis potentials this is the most likely form of evolution, rather than a state of continuing acceleration. Models of this type create a different theoretical framework in which to interpret the observation results of - and change their consequences for theories of cosmic structure formation. Other escape routes from acceleration, like introducing a quintessence field that ultimately decays into matter and radiation after producing acceleration, are at present entirely ad hoc. We can also imagine a richer cosmological structure if we admit the possibility of a random variation of $`V(\varphi )`$ around the universe, in the spirit of chaotic inflationary universes. The condition for transient accelerated expansion will arise quasi-randomly around the universe as a set of parameters like $`\{A,\varphi _0,\lambda \}`$ varies stochastically. This will give rise to some environments displaying transient acceleration at different times in their history, in which galaxy formation and clustering can continue to develop for longer than in the rarer, perpetually-accelerating, regions. The ability of expanding regions to recover from the onset of accelerated vacuum-dominated expansion may be an important factor in the global evolution of regions of an infinite universe containing long-lived stars and galaxies.


text:
# Reply to Boglione and Pennington

## Abstract

In a recent article we presented an argument which, we believe, shows to be incorrect an estimate, by Boglione and Pennington, of corrections to the valence (quenched) approximation predictions for properties of the lightest scalar glueball. Boglione and Pennington’s reply to our article, it appears to us, fails to address the specific technical issues we raised.

In Ref.  a formula is assumed for the errors, in valence (quenched) approximation scalar glueball predictions, which arise from the valence approximation’s omission of glueball-quarkonium mixing. No derivation from QCD itself is given to suggest why this particular formula for valence approximation errors might be expected to hold. The errors are then evaluated in a model. In Ref.  we raised two main objections to the work of Ref.  and, as a consequence of these objections, concluded that the claims of Ref.  to have shown the valence approximation to be unreliable for properties of the lightest scalar glueball are not correct. Boglione and Pennington’s  reply to our article mentions only one of our two objections and, it appears to us, fails to address the key technical issue on which this objection rests. To the second problem discussed in Ref. , Boglione and Pennington’s reply offers no response.

In the present comment we will give a brief qualitative summary of the issues raised in Ref. , reply to Boglione and Pennington’s  response to one of these issues, and, finally, consider an additional question concerning the work of Ref.  which arises from the numerical lattice calculation of glueball-quarkonium mixing in Ref. .

One of the objections to Ref.  discussed in Ref.  is a by-product of a systematic expansion derived in Ref.  for full QCD vacuum expectation values. The first term in this expansion is the valence approximation. Higher terms then give the valence approximation’s error. Comparing these terms with the formula assumed for valence approximation errors in Ref. , we show in Ref.  that the error formula of Ref.  appears not to be a correct expression for the error in the valence approximation as it is generally applied.

The difficulty with the error formula of Ref.  has a simple origin. In lattice calculations for the valence approximation to QCD and for full QCD, with a common choice of lattice spacing and a common choice of renormalization conditions, the valence approximation QCD coupling constant $`g_{val}`$ and the full QCD coupling constant $`g`$ are not equal. The valence approximation $`g_{val}`$ may be thought of as $`g`$ divided by a chromoelectric analog of a dielectric constant taking into account the screening of $`g`$ by dynamical quark-antiquark pairs present in full QCD but absent in the valence approximation. In the expansion for full QCD of Ref. , the difference between $`g`$ of full QCD and $`g_{val}`$ of the expansion’s first term is compensated by counter-terms subtracted from the expansion’s higher terms which contain closed quark loops. By choosing $`g_{val}`$ entering the leading term so that full QCD and the valence approximation obey the same renormalization conditions, the valence approximation’s accuracy is maximized and, by means of the subtracted counter-terms, the correction terms are minimized. Phrased differently, an optimally chosen $`g_{val}`$ shifts as large a contribution as possible out of the expansion’s higher order correction terms and into the expansion’s first term.

In fact, however, the expansion remains logically correct for any choice of $`g_{val}`$. In particular, $`g_{val}`$ could be set equal to the full QCD $`g`$. With this choice of $`g_{val}`$, the valence approximation becomes less accurate and the expansion’s correction terms become larger. In particular, the counter-terms subtracted from quark loops in the error terms of the expansion entirely vanish. In the relation assumed in Rev.  for valence approximation errors, no counter-terms accompany quark loops in the error terms. As a consequence of this absence, the version of the valence approximation to which this formula applies must have $`g_{val}`$ equal to $`g`$. The resulting valence approximation will be unnecessarily inaccurate and will differ from the valence approximation as generally applied in lattice QCD. Correspondingly, the accompanying error terms will be unnecessarily large and will not, as claimed in Ref. , be correct estimates of the errors in the valence approximation with an optimal choice of $`g_{val}`$. For example using data of Ref.  for inverse lattice spacing of about 1 GeV, with $`g_{val}`$ forced to $`g`$, light hadron masses and meson decay constants differ from experiment by as much as 45%, rather than by less than 10% or less than 20%, respectively, with an optimal choice of $`g_{val}`$.

Further evidence mentioned in Ref.  that the approximation to QCD which Ref.  identifies as the valence approximation uses a poorly chosen $`g_{val}`$ is provided by the sign predicted in Ref.  for the error in valence approximation mass predictions. The valence approximation with correctly chosen $`g_{val}`$, as discussed in Ref. , will generally underestimate excited state masses and meson decay constants. This expectation is supported by the data of Refs. . On the other hand, consider full QCD expanded according to Ref.  but with the valence approximation first term using full QCD $`g`$ in place of the optimal screened $`g_{val}`$. As expected for the effect of quark-antiquark color screening and easily confirmed by an inspection of lattice data, $`g_{val}`$ is less than $`g`$. Thus the mass predictions of this first term will be the valence approximation’s predictions if, with no change in the value of lattice spacing, $`g_{val}`$ is pushed up to $`g`$. A simple agrument based on asymptotic freedom or another direct inspection of lattice data show that raising $`g_{val}`$ to $`g`$ increases mass predictions. Thus while the valence approximation with optimal $`g_{val}`$ underestimates excited state masses, the valence approximation with $`g`$ in place of $`g_{val}`$ overestimates masses. This second alternative is the error found in the mass predictions of the proposed valence approximation of Ref. .

Boglione and Pennington  offer, in effect, a two part response to the preceding discussion. First, they argue that there is a complex relation between their expansion, relying on an effective theory of interacting hadrons, and QCD itself, as an interaction of quarks and gluons. As a consequence, formulating the counter-terms of Ref.  in the language of Ref.  and identifying where they might occur in the arithmetic of Ref.  is difficult. Second, they argue that, by the nature of its construction as an effective field theory, their theory will somewhere include the counter-terms of Ref.  in so far as these are actually components of QCD. In reply to the first of these comments, however, we offer the observation of Ref.  that the counter-term to the one-loop quark diagram is simply a derivative of the leading valence approximation term with respect to $`g_{val}`$. Thus by expressing the parameters of the effective theory as functions of $`g_{val}`$ it should, in principle, be possible to express the one-loop counter-term in an effective field theory even though quark and gluon fields themselves are submerged. In reply to the second comment, meanwhile, we point out that the question is not whether the counter-terms are implicitly present but rather exactly where they are present. With a correctly chosen $`g_{val}`$ the counter-terms appear subtracted from valence approximation corrections and shift part of the would-be corrections forward into the valence approximation. With a $`g_{val}`$ set equal to $`g`$, however, the counter-terms are missing from the quark-loop corrections and have been hidden, implicitly, in the valence approximation itself. The valence approximation then becomes less accurate and the corrections become larger. It is this second state of affairs which we believe we have shown occurs in the work of Ref. . Whatever errors Ref.  may find in the valence approximation as formulated by them do not bear on the valence approximation as generally applied.

The other objection raised in Ref.  can be recast as the claim that the effective field theory on which the calculation of Ref.  is based, according to Boglione and Pennington’s reply, is missing an important set of terms coupling glueballs and quarkonium. Boglione and Pennington’s reply proposes that the lightest scalar quarkonium and glueball fields entering the Lagrangian of their effective field theory correspond to orthogonal states, respectively a pure, ideally mixed quarkonium nonet and a pure unmixed glueball. Mixing among the glueball and quarkonium states then occurs only through decay couplings linear in the scalar fields and second order in pseudoscalar quarkonium fields. Missing from this picture, however, are bilinear terms in the effective Lagrangian directly coupling the scalar glueball field to scalar quarkonium fields. Such terms are the effective field theory’s remnants of the QCD quark-antiquark annihilation process for which a typical Feynman diagram is shown in Figure 1. The coefficients of these terms can not be determined from an effective field theory and must be taken as additional inputs from some other source. These terms are calculated using lattice QCD in Ref.  and shown to make significant contributions to quarkonium-glueball mixing. The mixing arising from these terms explains several otherwise puzzling features of observed scalar glueball and quarkonium data. It is shown further in Ref.  that mixing through these terms is quite probably more important than mixing through pseudoscalar pairs considered in Ref.  if mixing through pairs is calculated according to the systematic expansion of Ref. . To the objection that Ref.  omits the dominant mixing process, Boglione and Pennington’s reply offers no response.

A third issue concerning Ref.  arises from the calculation in Ref.  of the consequences of the mixing process of Fig. 1 for the glueball decay couplings calculated in Ref.  and taken as input by Ref. . It is explained in Ref.  that the coupling calculation of Ref.  includes the first order effect of quarkonium-glueball mixing by the process of Fig. 1. If this first order effect is removed, Ref.  shows, it appears possible that the remaining coupling of the pure pseudoscalar glueball to pseudoscalar pairs may be zero. In this case, glueball-quarkonium mixing by the mechanism considered in Ref.  would not only be small in comparison to the contribution of Fig. 1, it would be zero.